王小慶， 金先龍， 王建煒

(1.上海交通大學(xué) 機(jī)械系統(tǒng)與振動(dòng)國(guó)家重點(diǎn)實(shí)驗(yàn)室，上海 200240； 2.上海交通大學(xué) 機(jī)械與動(dòng)力工程學(xué)院，上海 200240)

土體-結(jié)構(gòu)非線性耦合系統(tǒng)動(dòng)力響應(yīng)并行計(jì)算方法研究

王小慶1,2， 金先龍1,2， 王建煒2

(1.上海交通大學(xué) 機(jī)械系統(tǒng)與振動(dòng)國(guó)家重點(diǎn)實(shí)驗(yàn)室，上海 200240； 2.上海交通大學(xué) 機(jī)械與動(dòng)力工程學(xué)院，上海 200240)

針對(duì)土體-結(jié)構(gòu)非線性耦合(Soil-Structure Interaction，SSI)系統(tǒng)動(dòng)力響應(yīng)數(shù)值模擬帶來的大規(guī)模計(jì)算量問題，提出基于SSI負(fù)載均衡及對(duì)偶圖理論兩種區(qū)域分解算法的并行計(jì)算方法。結(jié)合傳統(tǒng)的貪婪法及遞歸坐標(biāo)對(duì)分方法，對(duì)這四種方法的并行性能進(jìn)行研究。SSI采用基于對(duì)稱罰函數(shù)的方法處理，系統(tǒng)方程采用顯式中心差分有限元方法求解。對(duì)典型的SSI工程問題動(dòng)力響應(yīng)進(jìn)行并行數(shù)值模擬，并對(duì)這四種方法的可擴(kuò)展性進(jìn)行分析。結(jié)果表明：基于SSI負(fù)載均衡的并行計(jì)算方法，充分考慮土體和結(jié)構(gòu)耦合負(fù)載的均衡，并行效率最優(yōu)，基于對(duì)偶圖理論區(qū)域分解的方法和遞歸坐標(biāo)對(duì)分方法效率次之，貪婪法并行效率最低；隨核數(shù)增加，并行效率下降，需根據(jù)實(shí)際模型規(guī)模合理選擇并行計(jì)算核數(shù)，獲得最優(yōu)的并行計(jì)算效益；基于罰函數(shù)的顯式有限元方法能夠較好的解決SSI動(dòng)力響應(yīng)問題。

結(jié)構(gòu)-土體非線性耦合；區(qū)域分解；并行計(jì)算；動(dòng)力響應(yīng)；顯式有限元；對(duì)稱罰函數(shù)法

土體-結(jié)構(gòu)耦合(Soil-Structure Interaction，SSI)問題廣泛存在于各種工程系統(tǒng)中，如土體-隧道、列車-隧道-土體、土體-高層建筑、土體-橋梁、土體-核島等耦合系統(tǒng)，因而得到學(xué)者的廣泛研究[1-4]。隨著問題規(guī)模的擴(kuò)大，并行計(jì)算成為解決大型土體-結(jié)構(gòu)耦合系統(tǒng)動(dòng)力響應(yīng)的有效手段，各種并行算法成為研究熱點(diǎn)[4-7]，我國(guó)超級(jí)計(jì)算機(jī)的蓬勃發(fā)展與領(lǐng)先，也為并行計(jì)算提供了必需的硬件基礎(chǔ)。

土體-結(jié)構(gòu)耦合系統(tǒng)，就是將土體和所研究結(jié)構(gòu)視為彼此協(xié)調(diào)工作的整體，考慮結(jié)構(gòu)與土體動(dòng)力響應(yīng)的相互作用，進(jìn)而求解系統(tǒng)方程獲得變形和內(nèi)力等響應(yīng)。從研究方法上看，可分為理論研究[8]、實(shí)驗(yàn)研究及數(shù)值模擬[9]。三種方法各有優(yōu)勢(shì)，數(shù)值模擬技術(shù)因其研究周期短、可重復(fù)性強(qiáng)、可處理非常復(fù)雜的問題，應(yīng)用越來越廣泛。從研究?jī)?nèi)容上看，有單獨(dú)對(duì)各類土體-結(jié)構(gòu)耦合系統(tǒng)特性的研究[10]，有考慮土體-結(jié)構(gòu)耦合的大型工程分析[11-13]，前者的研究最終也是為輔助其在實(shí)際工程中得到合理的應(yīng)用?？紤]土體-結(jié)構(gòu)耦合的大型工程數(shù)值模擬帶來精細(xì)結(jié)果的同時(shí)，也伴隨著模型自由度及求解方程規(guī)模的急劇增加，為加快求解速度，提高計(jì)算效率，各種并行計(jì)算方法應(yīng)運(yùn)而生，YERLI等[4]在分布式計(jì)算機(jī)上利用子結(jié)構(gòu)法和有限-無限單元法提出土體-結(jié)構(gòu)耦合問題的并行算法，并加以程序?qū)崿F(xiàn)，通過算例證明了所提并行算法的有效性。GENES等[5-6]利用比例邊界有限元法模擬遠(yuǎn)場(chǎng)無限土體，采用有限元法模擬近場(chǎng)土體及結(jié)構(gòu)，在二者基礎(chǔ)上提出土體-結(jié)構(gòu)耦合模型并將其并行化，通過求解現(xiàn)有文獻(xiàn)中的算例，并和文獻(xiàn)中已有的結(jié)果進(jìn)行對(duì)比，驗(yàn)證其模型及并行程序的準(zhǔn)確性。MESCHKE等[7]利用有限元法提出隧道機(jī)械化掘進(jìn)中樁-土耦合作用的并行數(shù)值模擬方法，并對(duì)基于OpenMP的共享內(nèi)存式及基于MPI的分布內(nèi)存式并行技術(shù)進(jìn)行對(duì)比分析，結(jié)果表明基于區(qū)域分解的分布式并行計(jì)算能夠獲得更好的并行效率。目前針對(duì)土體-結(jié)構(gòu)耦合并行計(jì)算的研究多是基于不同的數(shù)值方法與一些特定應(yīng)用對(duì)象進(jìn)行，大多未考慮分區(qū)優(yōu)劣對(duì)并行效率的影響。

本文提出針對(duì)大型土體-結(jié)構(gòu)耦合系統(tǒng)動(dòng)力學(xué)分析的土體-結(jié)構(gòu)耦合負(fù)載均衡及基于對(duì)偶圖理論區(qū)域分解的并行計(jì)算方法。土體-結(jié)構(gòu)耦合采用基于對(duì)稱罰函數(shù)的方法處理，系統(tǒng)方程采用顯式中心差分有限元方法求解。通過工程實(shí)例與基于傳統(tǒng)的貪婪法和坐標(biāo)遞歸對(duì)分的并行計(jì)算方法進(jìn)行對(duì)比，評(píng)價(jià)各種方法的并行性能，為類似研究提供借鑒。

1 SSI動(dòng)力響應(yīng)問題基本數(shù)值方法

1.1 顯式有限元方法

土體-結(jié)構(gòu)耦合系統(tǒng)的動(dòng)力響應(yīng)，如地震動(dòng)響應(yīng)，列車動(dòng)載激勵(lì)下的響應(yīng)等，可利用顯式有限元方法求解，耦合系統(tǒng)動(dòng)力學(xué)方程可表示為：

(1)

土體-結(jié)構(gòu)耦合動(dòng)力學(xué)問題，因存在接觸而具有高度非線性特點(diǎn)，易采用顯式算法求解式(1)。顯式中心差分法無需組集總體剛度矩陣，當(dāng)質(zhì)量和阻尼陣為對(duì)角陣時(shí)無需矩陣分解和求逆，可節(jié)省內(nèi)存和計(jì)算機(jī)時(shí)，也更適合并行化。利用中心差分在時(shí)域內(nèi)離散：

(2)

(3)

式中:Δtn+1/2=(Δtn+Δtn+1)/2。當(dāng)0～t時(shí)刻節(jié)點(diǎn)位移、速度及加速度已知時(shí)，利用上述遞推公式可得到整個(gè)時(shí)域內(nèi)系統(tǒng)的響應(yīng)。顯式中心差分法是條件穩(wěn)定的，穩(wěn)定時(shí)間步長(zhǎng)由最小單元特征長(zhǎng)度及當(dāng)前波速?zèng)Q定。

1.2 基于對(duì)稱罰函數(shù)的SSI算法

土體與結(jié)構(gòu)耦合，存在相互擠壓和滑動(dòng)等行為[13]，屬邊界非線性問題，接觸算法比節(jié)點(diǎn)重合更符合物理實(shí)際。研究表明基于對(duì)稱罰函數(shù)的接觸算法適于土體-結(jié)構(gòu)耦合動(dòng)力響應(yīng)問題的處理。

利用基于對(duì)稱罰函數(shù)的接觸算法處理土體-結(jié)構(gòu)耦合問題，涉及接觸搜索算法和接觸力計(jì)算方法兩個(gè)方面。接觸搜索可采用基于段的Bucket分類搜索算法，其基本原理是將接觸空間劃分為多個(gè)很小的Bucket,然后，根據(jù)每個(gè)從節(jié)點(diǎn)所屬Bucket搜索并確定最近的主段，從而進(jìn)行接觸判斷。接觸力計(jì)算則采用對(duì)稱罰函數(shù)法，該方法動(dòng)量守恒準(zhǔn)確，易于程序?qū)崿F(xiàn)。其基本原理是：每一時(shí)步檢查各從節(jié)點(diǎn)是否穿透由接觸搜索算法確定的主段，若無穿透，則不作處理；若發(fā)生穿透，則在從節(jié)點(diǎn)與被穿透主段之間引入界面法向接觸力，如式(4)。從節(jié)點(diǎn)處理完成后，對(duì)各主節(jié)點(diǎn)進(jìn)行一遍相同的操作。

fi=lλkiηi

(4)

式中:l為判斷是否穿透的參數(shù)，λ為罰函數(shù)因子，ki為罰函數(shù)剛度，ηi為節(jié)點(diǎn)i的穿透量。接觸面的摩擦力利用經(jīng)典Coulomb摩擦理論計(jì)算。

2 SSI動(dòng)力響應(yīng)問題并行計(jì)算方法

2.1 區(qū)域分解技術(shù)

區(qū)域分解是實(shí)現(xiàn)并行計(jì)算的必要前處理，其思想最早由Schwarz提出，隨著高性能計(jì)算的發(fā)展，區(qū)域分解技術(shù)得到廣泛研究和應(yīng)用。其數(shù)學(xué)描述為，假設(shè)計(jì)算區(qū)域?yàn)棣福瑒t子區(qū)域Pi?Ω,i=1,…,NP的集合P={Pi}i=1,…,NP為Ω的一個(gè)區(qū)域分解，當(dāng)：

(5)

區(qū)域分解的邊界為集合：

(6)

區(qū)域分解技術(shù)可分為重疊和非重疊區(qū)域分解兩類，其中非重疊區(qū)域分解更適合并行計(jì)算體系結(jié)構(gòu),易于處理大型工程的并行數(shù)值分析[14]。有限元網(wǎng)格的非重疊區(qū)域分解有兩種即基于節(jié)點(diǎn)的區(qū)域分解和基于單元的區(qū)域分解，分別對(duì)應(yīng)于圖分區(qū)中對(duì)偶圖分解和節(jié)點(diǎn)圖分解，如圖1。基于節(jié)點(diǎn)的區(qū)域分解(圖1(a))，每個(gè)單元有唯一分區(qū)歸屬，相鄰分區(qū)的共享節(jié)點(diǎn)則復(fù)制并分配到所有共享分區(qū)，各鄰接分區(qū)間通過共享節(jié)點(diǎn)通信；基于單元的區(qū)域分解(圖1(b))，每個(gè)節(jié)點(diǎn)有唯一分區(qū)歸屬，相鄰分區(qū)的共享單元被復(fù)制并分配到所有共享分區(qū)，各鄰接分區(qū)通過共享單元通信。在不考慮接觸等額外計(jì)算負(fù)載時(shí)，顯式有限元?jiǎng)恿Ψ治龅臅r(shí)間基本消耗在單元計(jì)算上，如內(nèi)力更新等，同時(shí)相鄰分區(qū)間共享節(jié)點(diǎn)的復(fù)制更加容易，因而在土體-結(jié)構(gòu)耦合動(dòng)力分析中采用基于節(jié)點(diǎn)的非重疊區(qū)域分解技術(shù)。

圖1 有限元網(wǎng)格非重疊區(qū)域分解Fig.1 Non-overlapping domain decomposition of FE mesh

分區(qū)結(jié)果的優(yōu)劣直接關(guān)系后續(xù)并行求解的效率，通常可通過兩方面評(píng)判即各個(gè)子區(qū)域的計(jì)算負(fù)載是否均衡以及各分區(qū)間的通信消耗及整體通信消耗是否最小[15]。良好的負(fù)載均衡減少或者消除各處理器核的空閑等待時(shí)間，盡量少的通信可以降低通信在整個(gè)并行計(jì)算中所占時(shí)間比例，二者都可以提高并行計(jì)算效率。

2.2 區(qū)域分解方法

區(qū)域分解的實(shí)現(xiàn)可采用多種方法，如基于貪婪法的區(qū)域分解(Greedy Domain Decomposition,GDD)、遞歸坐標(biāo)對(duì)分方法(Recursive Coordinate Bisection,RCB)、最小和最大極慣性軸方法，遞歸譜對(duì)分方法等，此處重點(diǎn)介紹GDD方法、RCB方法及本文提出的土體-結(jié)構(gòu)耦合負(fù)載均衡區(qū)域分解方法(Soil-Structure Interaction Balanced Domain Decomposition,SSIB)和基于對(duì)偶圖的區(qū)域分解方法(Dual-Graph Based Domain Decomposition,DGB)。

貪婪算法是遞歸利用當(dāng)前局部最優(yōu)解,逐步得到或者逼近全局最優(yōu)解的一種求解最優(yōu)化問題的策略。基于貪婪法的區(qū)域分解GDD，即是利用貪婪算法的理念完成有限元網(wǎng)格等的區(qū)域分解。定義節(jié)點(diǎn)的權(quán)重為該節(jié)點(diǎn)所連接的單元數(shù)；模型總單元數(shù)為Nt；目標(biāo)分區(qū)數(shù)為Np。GDD的貪婪策略為每次選擇當(dāng)前最小權(quán)重的節(jié)點(diǎn)所連接的所有單元為初始分區(qū)，然后將已有單元的鄰接單元擴(kuò)展進(jìn)分區(qū)，遞歸執(zhí)行，直到分區(qū)內(nèi)單元數(shù)達(dá)到Nt/Np，然后進(jìn)行下一個(gè)分區(qū)。GDD的基本執(zhí)行流程為：首先，選擇上一分區(qū)Pi-1邊界中權(quán)重最小的節(jié)點(diǎn)ni，初始整體網(wǎng)格為分區(qū)P0；接著，將和該節(jié)點(diǎn)所有相鄰的未標(biāo)記的單元?jiǎng)潥w該分區(qū)Pi；然后，對(duì)分區(qū)Pi中所有已有的單元ek，遞歸執(zhí)行：①標(biāo)記ek，②將組成ek的所有節(jié)點(diǎn)的權(quán)重減1，③將所有和ek鄰接的未標(biāo)記的單元?jiǎng)澖o分區(qū)Pi；最后，當(dāng)分區(qū)Pi單元數(shù)達(dá)到Nt/Np，完成該分區(qū)，并進(jìn)入下一分區(qū)，直到分區(qū)完成。GDD方法原理簡(jiǎn)單，執(zhí)行高效，但局限于貪婪法的局部最優(yōu)特點(diǎn)，其分區(qū)難以達(dá)到全局最優(yōu)，對(duì)于土體-結(jié)構(gòu)耦合問題，并不能保證分區(qū)間土體-結(jié)構(gòu)耦合負(fù)載的均衡。

RCB方法由BERGER和BOKHARI在1987年提出，其分區(qū)基本原理為：首先，判斷模型在三個(gè)坐標(biāo)方向的長(zhǎng)度，沿垂直于最長(zhǎng)坐標(biāo)方向?qū)⒛Ｐ鸵环譃槎?，被分割成的兩部分分別包含網(wǎng)格總數(shù)的一半；然后，判斷新模型在三個(gè)坐標(biāo)方向的長(zhǎng)度，仍然沿垂直于最長(zhǎng)坐標(biāo)方向?qū)⑿履Ｐ鸵环譃槎蛔詈?，按照上面的策略，依次將模型分割直到滿足分割收斂條件。RCB方法簡(jiǎn)單高效，無需對(duì)整體模型進(jìn)行重新分區(qū)即能方便進(jìn)行局部負(fù)載重平衡。但對(duì)幾何形狀復(fù)雜的區(qū)域，RCB方法會(huì)產(chǎn)生不連續(xù)分區(qū)和極小分區(qū)，影響后續(xù)并行計(jì)算效率。同時(shí)，RCB方法僅根據(jù)幾何模型信息進(jìn)行分區(qū)，未考慮模型的載荷分布和接觸位置，對(duì)于土體-結(jié)構(gòu)耦合的大型結(jié)構(gòu)動(dòng)力響應(yīng)問題，由于存在大量接觸搜索計(jì)算，該方法并不能保證各子區(qū)域計(jì)算負(fù)載的均衡。

針對(duì)土體-結(jié)構(gòu)耦合問題的特點(diǎn)，本文提出基于耦合負(fù)載均衡的區(qū)域分解方法SSIB，其基本思想是保證各區(qū)域土體-結(jié)構(gòu)耦合計(jì)算負(fù)載均衡，從而使得各子區(qū)域總計(jì)算負(fù)載達(dá)到均衡。SSIB分區(qū)基本原理為：首先搜索所有參與土體-結(jié)構(gòu)耦合的接觸運(yùn)算單元及鄰接單元；接著根據(jù)參與并行計(jì)算的核數(shù)將所搜索的單元均分，并劃定區(qū)域邊界；然后，根據(jù)分區(qū)邊界，將邊界共享節(jié)點(diǎn)復(fù)制，并分配到各共享分區(qū)；進(jìn)而，將其余未參與耦合的單元按照空間坐標(biāo)位置分布到對(duì)應(yīng)子區(qū)域中，直到滿足分割收斂條件；最后將各分區(qū)模型導(dǎo)入求解器進(jìn)行并行計(jì)算。SSIB方法能在保證各子區(qū)域單元數(shù)目基本相同的前提下，將需要大量計(jì)算資源的土體-結(jié)構(gòu)耦合運(yùn)算平均分配，保證各子區(qū)域總體負(fù)載均衡。

基于圖理論的區(qū)域分解，能夠處理任意復(fù)雜的模型，同時(shí)易于并行化，在大規(guī)模并行有限元前處理中應(yīng)用愈加廣泛。本文提出基于對(duì)偶圖的區(qū)域分解方法DGB，其基本原理為：首先將有限元模型轉(zhuǎn)換為相應(yīng)的無向?qū)ε紙D；然后將生成的對(duì)偶圖利用圖分區(qū)理論進(jìn)行區(qū)域分解，本文采用K-路圖劃分算法[16]；最后將對(duì)偶圖分區(qū)結(jié)果映射還原為有限元網(wǎng)格的分區(qū)。圖2為有限元網(wǎng)格與其相應(yīng)對(duì)偶圖的映射關(guān)系。對(duì)偶圖在程序中利用在稀疏圖領(lǐng)域應(yīng)用廣泛的CSR (Compressed Storage)格式存儲(chǔ)，CSR利用兩個(gè)數(shù)組Index和Adjac存儲(chǔ)圖的鄰接信息，Index數(shù)組存儲(chǔ)每個(gè)頂點(diǎn)的開始和終了地址索引，Adjac數(shù)組存儲(chǔ)每個(gè)頂點(diǎn)的鄰接頂點(diǎn)，對(duì)于任意頂點(diǎn)，就可以根據(jù)Index中的地址索引在Adjac中找到其所有鄰接頂點(diǎn)。

圖2 有限元網(wǎng)格及其對(duì)偶圖Fig.2 FE mesh and itscorresponding dual-graph

2.3 基于區(qū)域分解的并行計(jì)算流程

利用Fortran語(yǔ)言結(jié)合MPI函數(shù)庫(kù)，實(shí)現(xiàn)所述區(qū)域分解方法的程序化，所開發(fā)程序可以讀入有限元輸入文件完成區(qū)域分解，輸出結(jié)果可以被后續(xù)求解程序讀入，完成并行計(jì)算求解。利用顯式求解技術(shù)以及各種區(qū)域分解方法，可以得到土體-結(jié)構(gòu)耦合問題動(dòng)力響應(yīng)并行求解流程，如圖3所示。

圖3 基于區(qū)域分解的顯式有限元并行計(jì)算流程Fig.3 Flow chart of the domain decomposition based parallel explicit FE solving procedure

本文研究工作并行計(jì)算環(huán)境為曙光5000A超級(jí)計(jì)算機(jī)(top500排名360，2014年11月)，其基本參數(shù)如表1所示。計(jì)算軟件為優(yōu)化的MPP版本LS-DYNA，采用所開發(fā)的區(qū)域分解程序?qū)ι婕巴馏w-結(jié)構(gòu)耦合的動(dòng)力學(xué)有限元模型進(jìn)行分區(qū)，然后進(jìn)行并行計(jì)算，評(píng)價(jià)各種區(qū)域分解方法對(duì)并行性能的影響。

表1 曙光5000A系統(tǒng)信息Tab.1 The system information of Dawning 5000A

本文采用強(qiáng)擴(kuò)展性評(píng)價(jià)各種方法的并行性能，評(píng)價(jià)指標(biāo)為加速比和并行效率，如式(7)所示，其中Sap為實(shí)際加速比，ts為單個(gè)核計(jì)算時(shí)耗費(fèi)時(shí)間，tm為采用多核計(jì)算時(shí)耗費(fèi)時(shí)間；Em為并行效率，Sip為理想加速比，對(duì)于起始為單核的情況下，Sip即為多核時(shí)的核數(shù)。

(7)

3 工程應(yīng)用實(shí)例

3.1 隧道地震動(dòng)響應(yīng)分析

隧道地震響應(yīng)中，隧道襯砌與周圍土體的作用屬于典型的土體-結(jié)構(gòu)耦合問題。本文以上海某輸水隧道為例，利用并行計(jì)算方法分析其在地震動(dòng)作用下的動(dòng)力響應(yīng)。圖4為有限元數(shù)值模型，隧道襯砌采用厚殼單元離散，土體采用8節(jié)點(diǎn)六面體單元模擬，模型節(jié)點(diǎn)總數(shù)約為190萬。

有限元模型中土體-結(jié)構(gòu)耦合采用前述基于對(duì)稱罰函數(shù)的接觸算法，土體-結(jié)構(gòu)耦合情況如圖4(c)所示，具體數(shù)值分析中土體和結(jié)構(gòu)接觸面分別設(shè)置為主和從接觸面，通過調(diào)整接觸控制參數(shù)，主要包括接觸剛度計(jì)算算法、接觸摩擦參數(shù)、罰因子、黏性接觸阻尼及最大穿透量等，提高土體結(jié)構(gòu)接觸模擬的準(zhǔn)確性。輸水隧道內(nèi)水采用附加質(zhì)量法模擬。土體模型根據(jù)實(shí)際勘測(cè)資料分為10層，基巖深度為300 m，土體本構(gòu)采用Biot滯回材料模型，隧道周圍土層的力學(xué)參數(shù)如表2所示，其中ρ為密度，E為彈性模量，ν為泊松比，γ為阻尼比，H為卓越頻率，w為剪切波速。采用PML人工邊界層，模擬無限區(qū)域，消除邊界截?cái)嗾`差。地震激勵(lì)采用50年超越概率10%的基巖加速度時(shí)程，時(shí)程曲線根據(jù)隧道場(chǎng)地基巖反應(yīng)譜等基巖地震動(dòng)參數(shù)擬合而成，并按照峰值加速度(PGA)0.1 g進(jìn)行調(diào)幅處理。地震動(dòng)采用一致激勵(lì)方式加載。

圖4 隧道地震動(dòng)響應(yīng)有限元數(shù)值模型Fig.4 FE numerical modelfor seismic response of tunnel

土層ρ/(kg·m-3)E/MPaνγH/Hzw/(m·s-1)11898380.260.0553.589.1121908690.280.0633.5118.8031796670.330.0783.5118.78417141160.350.0783.5158.40518372300.320.0703.5217.79618883710.280.0543.5277.18

圖5為地震激勵(lì)下，隧道某斷面直徑變化時(shí)程曲線，可以看到隧道截面呈現(xiàn)水平和豎直變化大小基本相同，符號(hào)相反的變化規(guī)律，這表明隧道呈現(xiàn)類橢圓變形，由圖5可知最大變形大約為1.2E-003，符合上海市《地基基礎(chǔ)設(shè)計(jì)規(guī)范》中隧道襯砌直徑變化的要求，表明隧道結(jié)構(gòu)地震作用下，直徑變化處于安全范圍內(nèi)。

圖5 典型斷面直徑變化率Fig.5 Diameter change rate of a typical section

圖6為采用不同的區(qū)域分解方法得到的16分區(qū)拓?fù)浣Y(jié)構(gòu)圖?？梢钥吹剑鞣N區(qū)域分解方法得到的分區(qū)拓?fù)洳槐M相同，GDD方法拓?fù)漭^不均勻，DGB方法產(chǎn)生較均勻的不規(guī)則分區(qū)，RCB方法沿兩個(gè)方向進(jìn)行遞歸對(duì)分，而SSIB方法基本沿隧道長(zhǎng)度方向分割，其根據(jù)土體-結(jié)構(gòu)耦合沿空間分布情況進(jìn)行區(qū)域分解，保證參與土體-結(jié)構(gòu)耦合單元在各子區(qū)域的平均分布。

圖6 土體-隧道耦合模型16分區(qū)拓?fù)銯ig.6 Topology of 16 partitions for the soil-tunnel coupling model

圖7給出了四種分區(qū)方法的加速比和并行效率的對(duì)比。由圖7(a)可見，在相同核數(shù)時(shí)，SSIB方法加速比最優(yōu)，最接近理想加速比，DGB及RCB方法次之，GDD方法加速比最小。隨著核數(shù)的增加，四種分區(qū)方法的差異也愈顯著，在32核時(shí)，SSIB和DGB表現(xiàn)出較明顯的高加速比。相應(yīng)地由圖7(b)可以看到，SSIB方法并行效率最高，DGB和RCB方法次之，GDD方法最差，可見考慮土體-結(jié)構(gòu)耦合對(duì)并行性能的影響。四種方法的并行效率，均隨核數(shù)的增加而降低，主要由于模型規(guī)模和拓?fù)涞脑?，隨著核數(shù)的進(jìn)一步提升，各分區(qū)間通信時(shí)間將增加，核數(shù)增加帶來的并行效率的提升被增加的通信時(shí)間所抵消，因而進(jìn)行并行計(jì)算時(shí)需要根據(jù)模型具體情況合理選擇并行核數(shù)。

圖7 不同分區(qū)方法加速比和并行效率對(duì)比Fig.7 Comparison of speed-up and parallel efficiency using different domain decomposition method

3.2 車-隧耦合動(dòng)態(tài)響應(yīng)分析

車-隧耦合振動(dòng)分析中也涉及到土體-結(jié)構(gòu)耦合問題。以上海某地鐵隧道為實(shí)例，研究各種區(qū)域分解方法的并行性能。根據(jù)設(shè)計(jì)資料和地質(zhì)勘測(cè)資料建立車-隧-土體耦合的三維有限元模型，模型節(jié)點(diǎn)數(shù)大約90萬。該隧道為雙管隧道，中間設(shè)置橫向聯(lián)絡(luò)通道，圖8為有限元數(shù)值模型。

圖8 車-隧-土體動(dòng)態(tài)耦合有限元數(shù)值模型Fig.8 Train-tunnel-soil dynamic coupling FE numerical model

土體-隧道耦合采用基于對(duì)稱罰函數(shù)的土體-結(jié)構(gòu)耦合算法，如圖8(c)為土體-結(jié)構(gòu)耦合模型，土體接觸面設(shè)為主接觸面，隧道接觸面設(shè)為從面，同時(shí)根據(jù)土體和結(jié)構(gòu)的相對(duì)尺寸及材料屬性，選擇最優(yōu)接觸控制參數(shù)，從而準(zhǔn)確模擬土體和隧道的耦合作用。模型采用分層土體(5層)，并利用D-P彈塑性本構(gòu)模型模擬土體非線性特性，關(guān)鍵土層物理力學(xué)參數(shù)如表3所示，其中ρ為密度，E為彈性模量，ν為泊松比，C為內(nèi)黏聚力，φ為內(nèi)摩擦角。采用黏彈性人工邊界模擬無限區(qū)域。列車采用10節(jié)編組，根據(jù)實(shí)際幾何信息，采用剛?cè)狁詈戏绞?車廂采用剛體模擬，懸架系統(tǒng)等采用柔性體模擬)建立有限元模型圖8(b)。軌道不平順采用美國(guó)六號(hào)軌道不平順譜模型。列車與軌道以接觸方式實(shí)現(xiàn)動(dòng)態(tài)耦合。公路車道載荷根據(jù)《公路隧道設(shè)計(jì)規(guī)范》施加。

表3 關(guān)鍵土層力學(xué)參數(shù)Tab.3 Mechanic parameters for key soil strata

模擬中首先采用動(dòng)力松弛方法獲得靜應(yīng)力場(chǎng)，并將其作為動(dòng)力分析的初始條件，然后軌道列車以90 km/h的速度勻速由外部駛?cè)氩⑼ㄟ^隧道。圖9為列車經(jīng)過時(shí)隧道某典型斷面拱頂和拱底的沉降曲線，可見車輛經(jīng)過時(shí)測(cè)點(diǎn)位移發(fā)生瞬時(shí)局部變化，引起隧道局部襯砌約1 mm的沉降，車輛經(jīng)過之后襯砌位移基本恢復(fù)到之前的狀態(tài)。圖10為列車經(jīng)過時(shí)某斷面的環(huán)向應(yīng)力變化情況，可見列車經(jīng)過時(shí)引起該斷面極值約0.13 MPa的附加應(yīng)力。由于10節(jié)車廂共20個(gè)懸架結(jié)構(gòu)，其中中間18個(gè)懸架兩兩在一起，因此軌道車輛經(jīng)過時(shí)，應(yīng)引起襯砌結(jié)構(gòu)最大主應(yīng)力約11次較明顯的震蕩，由圖10可明顯看出最大應(yīng)力出現(xiàn)了11次峰值，這與定性分析相符。

圖9 斷面沉降時(shí)程曲線Fig.9 Time history of the subsidence of a cross-section

圖10 斷面環(huán)向應(yīng)力時(shí)程曲線Fig.10 Time history of the hoop stress of a cross-section

圖11給出了32分區(qū)時(shí)不同分區(qū)方法產(chǎn)生的分區(qū)拓?fù)浣Y(jié)構(gòu)。和實(shí)例3.1類似，GDD的拓?fù)涑尸F(xiàn)不均勻的不規(guī)則分布，DGB的拓?fù)涑尸F(xiàn)較均勻的不規(guī)則分布，同時(shí)DGB和RCB方法呈現(xiàn)相似分區(qū)拓?fù)?；RCB每次對(duì)分均沿模型最長(zhǎng)軸方向，不考慮接觸負(fù)載分布，列車模型全部被分到一個(gè)子區(qū)域中，軌道與列車的接觸計(jì)算集中于一個(gè)區(qū)域內(nèi)；SSIB方法仍然根據(jù)土體-結(jié)構(gòu)耦合的分布情況將其均勻分部到各子區(qū)域，同時(shí)列車與隧道的耦合也均分到了相鄰幾個(gè)分區(qū)，使得總體計(jì)算量分配更加均衡。

圖11 列車-隧道-土體耦合模型32分區(qū)拓?fù)銯ig.11 Topology of 32 partitions for the train-tunnel-soil coupling model

圖12為四種不同方法的加速比和并行效率。由圖12(a)可見，四種方法的加速比呈現(xiàn)非線性變化，并逐漸遠(yuǎn)離相應(yīng)的理想加速比，ISSB方法在不同核數(shù)時(shí)均呈現(xiàn)較高的加速比，DGB和RCB方法次之，GDD方法加速比最小，并且這種趨勢(shì)隨著核數(shù)的增加而愈明顯，此外DGB和RCB方法的加速比比較接近。相應(yīng)地，圖12(b)中，并行效率由高到低依次為SSIB、DGB、RCB和GDD，同時(shí)可以看到四種方法的并行效率曲線的斜率逐漸變小，這也表明了隨著核數(shù)的增加，由于通信時(shí)間在總計(jì)算時(shí)間中所占比例的增加，而導(dǎo)致并行效率的逐漸降低。該實(shí)例的加速比和并行效率變化規(guī)律與實(shí)例3.1類似，這表明了本文所提出的基于土體-結(jié)構(gòu)耦合的負(fù)載均衡方法SSIB和基于對(duì)偶圖區(qū)域分解方法DGB的有效性。相比于圖7隧道地震動(dòng)響應(yīng)的加速比和并行效率。列車-隧道耦合動(dòng)力響應(yīng)的相應(yīng)加速比和并行效率更高，這是因?yàn)樵撃Ｐ椭羞€存在車輛和隧道的耦合接觸，其分布與土體-結(jié)構(gòu)耦合分布方向一致，采用SSIB方法也同時(shí)保證了車-隧耦合負(fù)載的均衡。

圖12 不同分區(qū)方法加速比和并行效率曲線Fig.12 Curves of speed-up and parallel efficiency using different domain decomposition method

4 結(jié) 論

對(duì)于含有大量土體-結(jié)構(gòu)耦合系統(tǒng)的動(dòng)力響應(yīng)數(shù)值模擬，并行計(jì)算應(yīng)用越來越多。本文針對(duì)此問題提出并程序?qū)崿F(xiàn)了基于土體-結(jié)構(gòu)耦合負(fù)載均衡SSIB及基于對(duì)偶圖理論區(qū)域分解DGB的并行計(jì)算方法。通過實(shí)例結(jié)合遞歸坐標(biāo)對(duì)分法RCB和貪婪法GDD，對(duì)這四種并行方法的并行性能進(jìn)行了研究，結(jié)果表明：

(1)SSIB方法充分考慮土體-結(jié)構(gòu)耦合單元負(fù)載的均衡并行效率最優(yōu)；

(2)DGB方法和RCB方法，雖未直接指定SSI負(fù)載的均衡，但其分區(qū)拓?fù)滹@示大部分子區(qū)域均有一定的SSI運(yùn)算單元，其并行效率次之；同時(shí)利用DGB方法分區(qū)時(shí)保證切邊最少，從而最大限度的降低通信量，其并行效率較優(yōu)于RCB方法；

(3)GDD方法負(fù)載分配最不均勻，其并行效率最差；

(4)隨核數(shù)的增多，分區(qū)間通信量增加，并行效率下降，需要根據(jù)問題規(guī)模合理選擇并行計(jì)算核數(shù)，獲得最優(yōu)的并行計(jì)算效益。

(5)通過基于罰函數(shù)的顯式有限元方法實(shí)現(xiàn)了兩類土體-結(jié)構(gòu)耦合問題動(dòng)力響應(yīng)的數(shù)值模擬，表明該方法能夠較好的處理這兩類問題。

未來的研究工作將考慮引入土體-結(jié)構(gòu)耦合負(fù)載均衡系數(shù)，對(duì)基于圖理論區(qū)域分解的并行計(jì)算方法進(jìn)行改進(jìn)，充分利用其界面通信量小的優(yōu)勢(shì)，通過調(diào)整土體-結(jié)構(gòu)耦合負(fù)載均衡系數(shù)使得其在各子區(qū)域負(fù)載均衡化，從而獲得更好的并行效率。

[1] LOU M, WANG H, CHEN X, et al. Structure-soil-structure interaction:literature review[J]. Soil Dynamics and Earthquake Engineering, 2011, 31(12): 1724-1731.

[2] 劉毅, 薛素鐸, 李雄彥. 土-結(jié)構(gòu)相互作用下網(wǎng)架結(jié)構(gòu)動(dòng)力性能研究[J]. 振動(dòng)與沖擊, 2014, 33(10): 22-28. LIU Yi, XUE Suduo, LI Xiongyan. Grid structure dynamic performance analysis considering soil-structure interaction [J]. Journal of Vibration and Shock, 2014, 33(10): 22-28.

[3] VASILEV G, PARVANOVA S, DINEVA P, et al. Soil-structure interaction using BEM-FEM coupling through ANSYS software package[J]. Soil Dynamics and Earthquake Engineering, 2015, 70: 104-117.

[4] YERLI H R, KACIN S, KOCAK S. A parallel finite-infinite element model for two-dimensional soil-structure interaction problems[J]. Soil Dynamics and Earthquake Engineering, 2003, 23(4): 249-253.

[5] GENES M C, KOCAK S. A combined finite element based soil-structure interaction model for large-scale systems and applications on parallel platforms[J]. Engineering Structures, 2002, 24(9): 1119-1131.

[6] GENES M C. Dynamic analysis of large-scale SSI systems for layered unbounded media via a parallelized coupled finite-element/ boundary-element/scaled boundary finite-element model[J]. Engineering Analysis with Boundary Elements, 2012, 36(5): 845-857.

[8] 鄒立華, 方雷慶. 考慮樁-土-結(jié)構(gòu)相互作用的結(jié)構(gòu)振動(dòng)控制研究[J]. 振動(dòng)與沖擊, 2010, 29(11): 100-104.

ZOU Lihua, FANG Leiqing. Structural vibration control considering pile-soil-structure dynamic interaction[J]. Journal of Vibration and Shock, 2010, 29(11): 100-104.

[9] ABATE G, MASSIMINO M R, MAUGERI M. Numerical modelling of centrifuge tests on tunnel-soil systems[J]. Bulletin of Earthquake Engineering, 2014,13(7): 1-25.

[10] LEE J H, KIM J K, KIM J H. Nonlinear analysis of soil-structure interaction using perfectly matched discrete layers[J]. Computers & Structures, 2014, 142: 28-44.

[11] 張巍, 陳清軍. 土-樁-結(jié)構(gòu)非線性相互作用體系行波效應(yīng)的并行計(jì)算分析[J]. 湖南大學(xué)學(xué)報(bào)(自然科學(xué)版), 2012, 39(6): 19-25. ZHANG Wei, CHEN Qingjun. Parallel computing analysis of the traveling wave effect of soil-pile-structure nonlinear interaction system finite element analysis of three-dimensional soil-structure interactionsystem[J]. Journal of Hunan University(Natural Sciences),2012, 39(6): 19-25.

[12] LI P, SONG E. Three-dimensional numerical analysis for the longitudinal seismic response of tunnels under an asynchronous wave input[J]. Computers and Geotechnics, 2015, 63: 229-243.

[13] DING J H, JIN X L, GUO Y Z, et al. Numerical simulation for large-scale seismic response analysis of immersed tunnel[J]. Engineering Structures, 2006, 28(10): 1367-1377.

[14] KRYSL P, BITTNAR Z. Parallel explicit finite element solid dynamics with domain decomposition and message passing: dual partitioning scalability[J]. Computers & Structures, 2001, 79(3): 345-360.

[15] FARHAT C. A simple and efficient automatic FEM domain decomposer[J]. Computers & Structures, 1988, 28(5): 579-602.

[16] KARYPIS G, KUMAR V. A fast and high quality multilevel scheme for partitioning irregular graphs[J]. SIAM Journal on Scientific Computing, 1998, 20(1): 359-392.

A study on parallel computing methods for dynamic response analysis of soil-structurenonlinear interaction systems

WANG Xiaoqing1,2, JIN Xianlong1,2, WANG Jianwei2

(1. State Key Laboratory of Mechanical System and Vibration, Shanghai Jiao Tong University，Shanghai 200240, China; 2. School of Mechanical Engineering, Shanghai Jiao Tong University，Shanghai 200240, China)

In order to solve the problem of massive amount of computation brought by the numerical dynamic response simulation of soil-structure nonlinear interaction (SSI) systems, a parallel computing method using SSI load balanced and dual-graph theory based domain decomposition (DD) algorithm was proposed. Combined with the traditional greedy and recursive coordinate bisection algorithm, the parallel performance of four algorithms was investigated. The SSI was dealt with the symmetric penalty method. The system equation was solved using the finite element method (FEM) with an explicit central difference scheme. The dynamic responses of typical engineering problems with SSI were simulated in parallel, and the scalability of the four algorithms was analyzed. The results indicate that the SSI load balanced algorithm which substantially balances the coupling loads of soil and structure shows the best parallel efficiency, followed by the dual-graph theory based algorithm and the recursive coordinate bisection algorithm. The greedy method has the lowest parallel efficiency. The parallel efficiency decreases with the number of cores, and the number of cores should be chosen properly according to the scale of the actual model to achieve the optimal parallel performance. The explicit FEM with the penalty method is a proper approach for SSI dynamic analysis.

soil-structure nonlinear interaction; domain decomposition; parallel computing; dynamic response; explicit FEM; symmetric penalty method

國(guó)家863高技術(shù)研究發(fā)展計(jì)劃(2012AA01AA307)；國(guó)家自然科學(xué)基金資助項(xiàng)目(11272214;51475287)

2015-08-11 修改稿收到日期：2015-12-08

王小慶 男，博士生，1988年生

金先龍 男，教授，1961年生

TB122

A

10.13465/j.cnki.jvs.2016.24.004