公鐵兩用隧道動(dòng)態(tài)響應(yīng)并行計(jì)算分析

張偉偉，金先龍，曹露芬，王建煒，王 新

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

隨著我國(guó)的城市化進(jìn)程的推進(jìn)，城市交通壓力激增，地鐵因其具有準(zhǔn)時(shí)、高效、運(yùn)量大以及占地少等優(yōu)勢(shì)，在大城市得到了迅速發(fā)展。同時(shí)，地鐵隧道建設(shè)也朝著多線、多層以及多用途的方向發(fā)展。對(duì)于這些具有復(fù)雜結(jié)構(gòu)和多種載荷作用的隧道結(jié)構(gòu)，如何更加準(zhǔn)確地評(píng)估其穩(wěn)定性，是隧道動(dòng)力響應(yīng)研究的熱點(diǎn)。

對(duì)于隧道結(jié)構(gòu)的動(dòng)力響應(yīng)的計(jì)算，數(shù)值模擬是一種高效的分析方法，國(guó)內(nèi)外學(xué)者已經(jīng)做出了許多有意義的研究。王祥秋等［1］提出了列車(chē)提速后的荷載模擬方法，并分析了京-廣線朱亭隧道以及珠三角地區(qū)典型地層隧道的動(dòng)力響應(yīng)特征［2-3］。白冰等［4］分析了地鐵振動(dòng)下平行隧道的彈塑性動(dòng)力響。龔倫等［5］用三維模型仿真了列車(chē)動(dòng)載引起的下穿隧道的振動(dòng)響應(yīng)。Gupta等［6］全面的模擬和分析了歐洲“大力士”列車(chē)穿越荷蘭Groene Hart隧道時(shí)，隧道及環(huán)境的動(dòng)力響應(yīng)情況。Degrande等［7］采用周期邊界元-有限元耦合模型預(yù)測(cè)了列車(chē)引起的隧道振動(dòng)。Andersena等［8］比較了二維與三維邊界元-有限元耦合模型的特點(diǎn)，與試驗(yàn)結(jié)果的對(duì)比表明：全三維模型可以更精確地得到振動(dòng)在隧道結(jié)構(gòu)中的傳遞情況，也能更好的反映不同埋深對(duì)隧道動(dòng)力響應(yīng)的影響。同時(shí)，有少數(shù)學(xué)者［9］也開(kāi)始著手研究準(zhǔn)靜態(tài)載荷作用對(duì)隧道結(jié)構(gòu)響應(yīng)的影響。目前，國(guó)內(nèi)研究者主要采用二維模型或三維模型的局部片段來(lái)模擬隧道，國(guó)外研究者主要采用周期性邊界元-有限元耦合模型模擬隧道，對(duì)于列車(chē)動(dòng)載荷的模擬多數(shù)采用簡(jiǎn)化的四分之一懸架模型。從研究現(xiàn)狀上看，隧道振動(dòng)數(shù)值模擬的模型往往是實(shí)際結(jié)構(gòu)的極度簡(jiǎn)化，不能夠全面的反映隧道以及地質(zhì)情況沿縱向的變化，載荷的模擬也主要采用解析解的形式。

本文以上海長(zhǎng)江隧道為背景，建立了軌道車(chē)輛-隧道-土體耦合非線性動(dòng)力學(xué)三維有限元模型，該模型考慮了聯(lián)絡(luò)通道、變形縫、車(chē)輛二級(jí)懸架等細(xì)節(jié)，節(jié)點(diǎn)數(shù)和單元數(shù)均超過(guò)百萬(wàn)。進(jìn)而，采用基于負(fù)載均衡的顯式并行有限元技術(shù)，在超級(jí)計(jì)算機(jī)上完成了求解。最后，分析了隧道結(jié)構(gòu)在公路車(chē)道準(zhǔn)靜態(tài)載荷以及列車(chē)動(dòng)載荷作用下的動(dòng)力響應(yīng)以及接觸負(fù)載均衡算法對(duì)計(jì)算效率的提升效果。

1 理論與方法

1.1 顯式中心差分方法

根據(jù)達(dá)朗貝爾原理，振動(dòng)系統(tǒng)的動(dòng)力學(xué)方程為：

式中，M為結(jié)構(gòu)的質(zhì)量矩陣，C為結(jié)構(gòu)的阻尼矩陣，K為結(jié)構(gòu)的剛度矩陣，F(xiàn)ext為外界作用力矢量，a為加速度列陣，v為速度列陣，d為位移列陣。

對(duì)于線性問(wèn)題，可用解析方法求解該方程。若是非線性問(wèn)題時(shí)，一般運(yùn)用數(shù)值方法求解。目前，較多采用的數(shù)值方法［10］有兩種：Newmark方法（隱式）和中心差分方法（顯式）。而對(duì)于存在較多內(nèi)部接觸的高度非線性問(wèn)題，隱式方法往往難以保證收斂。在LS-DYNA軟件中，采用中心差分格式的直接積分方法進(jìn)行時(shí)間積分，在0～tn時(shí)間步結(jié)果已知的條件下，方程（1）可改寫(xiě)成：

加速度可通過(guò)質(zhì)量矩陣的轉(zhuǎn)置矩陣與剩余載荷矢量相乘得到。若M為對(duì)角陣，它的轉(zhuǎn)置矩陣仍為對(duì)角陣，且矩陣方程可看作是每一個(gè)自由度上的獨(dú)立方程組：

此時(shí)，在時(shí)域內(nèi)應(yīng)用中心差分法得：

滿足顯式有限元方法穩(wěn)定條件的最大時(shí)間步長(zhǎng)Δt的值為：

式中，ωmax為線性系統(tǒng)的最大自然角頻率；le是單元e的特征長(zhǎng)度；ce是單元e的當(dāng)前波速；取α為考慮非線性導(dǎo)致的不穩(wěn)定性影響的折減系數(shù)；當(dāng)Δt=αΔtcrit時(shí)，求解是穩(wěn)定的。因此，需要采用合理時(shí)間步長(zhǎng)來(lái)控制顯式有限元求解的穩(wěn)定性，對(duì)有限元建模要求也較高。

1.2 接觸算法

相對(duì)于隧道土體之間節(jié)點(diǎn)重合，隧道與土體之間采用接觸進(jìn)行模擬，更符合物理實(shí)際，且有利于相互作用力和摩擦力的計(jì)算。接觸計(jì)算的核心問(wèn)題主要是：① 接觸搜索算法，② 接觸力的計(jì)算方法，即接觸面之間力的傳遞。

在接觸搜索算法上，采用基于段的Bucket分類(lèi)主從接觸面搜索算法。它的基本原理是：對(duì)于每個(gè)從節(jié)點(diǎn)，搜索最近的主段，首先把接觸空間劃分成很多小的Bucket，對(duì)每個(gè)從節(jié)點(diǎn)，按照它的空間位置指定于某個(gè)Bucket中，然后在該Bucket內(nèi)搜索最近的主段，運(yùn)行一段時(shí)間后，根據(jù)該節(jié)點(diǎn)空間位置的變化再重新指定新的Bucket。該方法可以避免由于主面網(wǎng)格復(fù)雜或不連續(xù)，以及接觸網(wǎng)格差別較大等引起的錯(cuò)誤接觸搜尋結(jié)果。

在接觸力的計(jì)算上，采用對(duì)稱(chēng)罰函數(shù)算法。其基本原理是：每一時(shí)步先檢查各從節(jié)點(diǎn)是否穿透主表面，沒(méi)有穿透則對(duì)該從節(jié)點(diǎn)不作如何處理。如果穿透，則在該從節(jié)點(diǎn)與被穿透主表面之間引入一個(gè)界面接觸力，其大小與穿透深度、主片剛度成正比，系數(shù)為罰函數(shù)值。即：

式中，K為主接觸面材料的體積模量，V0為接觸片的初始體積，Ai為接觸片面積，ψ為罰函數(shù)值，λi為節(jié)點(diǎn)i的穿透量。

1.3 基于接觸負(fù)載均衡的并行算法

由于隧道、土體模型巨大，同時(shí)，在隧道與土體以及輪、軌之間存在大量接觸，外加土體存在的材料非線性，從而導(dǎo)致了仿真模型的計(jì)算量非常巨大。若采用普通的串行計(jì)算機(jī)進(jìn)行求解，將消耗大量的計(jì)算時(shí)間。采用并行計(jì)算可以較好的解決這一問(wèn)題［11，12］。而區(qū)域分解法的則是主流的并行計(jì)算技術(shù)。

區(qū)域分解是指：將復(fù)雜模型按一定的原則劃分為若干個(gè)子區(qū)域，然后將各個(gè)子區(qū)域交給不同的處理器進(jìn)行分布計(jì)算，并在處理器之間建立合理的數(shù)據(jù)交換機(jī)制。該方法非常有利于并行計(jì)算的實(shí)現(xiàn)。需要考慮的問(wèn)題主要是各計(jì)算節(jié)點(diǎn)負(fù)載均衡以及降低各節(jié)點(diǎn)之間的通信開(kāi)銷(xiāo)。

根據(jù)仿真模型中存在大量接觸的特點(diǎn)，本文提出采用基于接觸負(fù)載均衡（Contact Balance Bisection Algorithm，簡(jiǎn)稱(chēng)CBB）的并行計(jì)算方法進(jìn)行求解。首先，對(duì)模型進(jìn)行坐標(biāo)變換，提高坐標(biāo)對(duì)分的精度，以降低通信開(kāi)銷(xiāo)；進(jìn)而，從模型中搜索接觸定義中的所有單元及鄰接單元，根據(jù)參與計(jì)算的CPU數(shù)量將接觸定義中的單元均分，并對(duì)變換后的模型采用遞歸坐標(biāo)對(duì)分法劃定區(qū)域邊界；然后，根據(jù)上一步劃定的邊界進(jìn)行分區(qū)，邊界節(jié)點(diǎn)被復(fù)制，分布到對(duì)應(yīng)的區(qū)域中；從而，將模型中非耦合接觸單元按空間坐標(biāo)分布到對(duì)應(yīng)子區(qū)域中，直到滿足接觸負(fù)載均衡條件；最后將分區(qū)模型導(dǎo)入求解器，由各分區(qū)指定處理器進(jìn)行計(jì)算。

1.4 靜應(yīng)力場(chǎng)的加載方法

由于隧道與土體之間采用接觸進(jìn)行模擬，加入靜應(yīng)力場(chǎng)有利于隧道與土體的接觸計(jì)算（避免“脫開(kāi)”而導(dǎo)致接觸搜索不到而造成計(jì)算誤差）。在LS-DYNA中實(shí)現(xiàn)靜應(yīng)力場(chǎng)的加載有多種方法：重啟動(dòng)方法，隱式顯式序列求解的方法，“Dynain”文件方法和動(dòng)力松弛方法等。動(dòng)力松弛方法可以有效克服前3種方法的缺點(diǎn)，但需要合理的選擇阻尼c使系統(tǒng)接近于臨界阻尼ccr，從而快速逼近靜態(tài)解。

由矩陣特征值相關(guān)理論，系統(tǒng)阻尼矩陣取值為：

系統(tǒng)最高頻率由模型特征尺寸決定，同上文：

而系統(tǒng)阻尼可動(dòng)態(tài)取為：

2 計(jì)算模型

上海長(zhǎng)江隧道長(zhǎng)約8.95 km，單管外徑15 m，內(nèi)徑13.7 m。隧道整體斷面設(shè)計(jì)為上下層雙線隧道，兩單管間凈距約為16 m，隧道上層設(shè)三條公路車(chē)道，下層設(shè)一條鐵路軌道。

2.1 隧道－土體模型

本文采用實(shí)體和厚殼單元建立了隧道-土體的動(dòng)力耦合模型，采用1.2節(jié)所述的動(dòng)態(tài)面-面接觸方法，在隧道和土體之間建立主從接觸面來(lái)模擬它們的相互作用。其有限元模型及局部細(xì)節(jié)如圖1所示。

圖1 隧道-土體耦合有限元模型Fig.1 Coupled tunnel-soil FE model

土體模型深度從地表直至海平面以下100 m，在土體寬度尺寸設(shè)計(jì)時(shí)，以隧道輪廓距土體邊界不小于2.5倍隧道截面寬度為標(biāo)準(zhǔn)，土體模型的寬度取為200 m。并根據(jù)實(shí)地測(cè)得的地質(zhì)，對(duì)土體模型了分層（5層）操作。采用德魯克-普拉格（Drucker-Prager）塑性本構(gòu)關(guān)系來(lái)模擬土體的非線性特性。土體相關(guān)參數(shù)見(jiàn)表1。

表1 土層物理力學(xué)參數(shù)Tab.1 Physical and mechanical parameters of soil strata

為了克服邊界對(duì)計(jì)算結(jié)果的影響，需要在土體四周施加人工邊界來(lái)模擬無(wú)限區(qū)域。LS-DYNA自帶了粘性邊界選項(xiàng)，但單獨(dú)使用粘性邊界具有精度不夠高和低頻飄移失穩(wěn)等缺點(diǎn)。本文采用劉晶波等［13］提出的粘彈性人工邊界對(duì)原有粘性邊界進(jìn)行改進(jìn)，其主要思想是在人工截?cái)噙吔缟显O(shè)置連續(xù)分布的并聯(lián)彈簧-阻尼器系統(tǒng)，該方法具有較好的效果和可行性。具體實(shí)現(xiàn)方法為在已經(jīng)建好的三維有限元模型周?chē)由煲欢ê穸鹊娜S單元作為等效一致粘彈性邊界單元，并將等效粘彈性人工邊界單元的最外層節(jié)點(diǎn)固定，并根據(jù)整體模型尺寸和等效單元厚度等參數(shù)設(shè)置其彈性模量。

在結(jié)構(gòu)動(dòng)力分析中，通常需要設(shè)置Rayleigh阻尼系數(shù)α和β。只要得到兩種振型下的振型頻率，并根據(jù)相關(guān)資料提供的阻尼比，便可得出α、β值。根據(jù)文獻(xiàn)［14］闡述的方法可求得土層的固有周期，進(jìn)而求得各土層的Rayleigh阻尼系數(shù)見(jiàn)表2。

表2 土層阻尼系數(shù)Tab.2 Damping parameters of soil strata

隧道模型含襯砌、車(chē)道板、道床、聯(lián)絡(luò)通道及變形縫等部分。隧道襯砌混凝土強(qiáng)度等級(jí)為C60，密度為2 500 kg/m3，彈性模量為 30 GPa，泊松比為 0.25。

2.2 車(chē)輛－軌道－道床動(dòng)力耦合模型

（1）車(chē)輛模型。根據(jù)的實(shí)際的幾何信息，本文建立了軌道車(chē)輛剛-柔耦合多體動(dòng)力學(xué)模型（僅考慮了車(chē)輛的垂向振動(dòng)），模型包括車(chē)體、轉(zhuǎn)向架、二級(jí)懸架以及輪對(duì)等部分。該軌道車(chē)輛包含10節(jié)編組，車(chē)廂采用剛性體模擬，二級(jí)懸掛系統(tǒng)采用等效線性彈簧和阻尼單元模擬。其主要?jiǎng)恿W(xué)參數(shù)見(jiàn)表3。

表3 軌道車(chē)輛模型參數(shù)Tab.3 Parameters of the train model

（2）軌道模型。軌道采用空間平行的兩排梁?jiǎn)卧M(jìn)行模擬，因輪對(duì)與軌道接觸剛度通常較大，軌道空間曲線的劇烈變化會(huì)導(dǎo)致較大的接觸力。所以，軌道軌跡曲線采用單元幾何位置與不平度疊加后擬合得到的光滑曲線來(lái)模擬［15］。軌道的不平順設(shè)置為美國(guó)6級(jí)軌道不平順譜，其高低不平順功率譜表達(dá)式為：

其中，Ω為不平順的空間角頻率，Ωc為截?cái)囝l率，取值為 0.824 5，k為安全系數(shù)，取值為 0.75，Av為粗糙度常數(shù)，取值為0.033 9。空間域的高低不平順由傅里葉逆變換得到。

（3）道床模型。由于隧道結(jié)構(gòu)凈空間限制，道床可用空間較小，故選用了現(xiàn)澆混凝土無(wú)軌枕式整體道床，其混凝土強(qiáng)度為C40。仿真計(jì)算時(shí)采用體單元對(duì)其進(jìn)行建模，并與隧道仰拱固連。

（4）耦合模型。輪軌之間采用單點(diǎn)罰函數(shù)接觸進(jìn)行模擬，以保證車(chē)輪沿軌道行進(jìn)，垂向和橫向力單獨(dú)計(jì)算，不考慮軌道蠕滑等效應(yīng)［15］。由于采用了無(wú)軌整式整體道床，其對(duì)軌道的支承基本可視為連續(xù)性的。建模時(shí)，梁?jiǎn)卧c道床單元采用節(jié)點(diǎn)重合，將墊片、扣件等的彈性和阻尼效應(yīng)等效到道床單元的參數(shù)中，類(lèi)似于文獻(xiàn)［16］描述的連續(xù)彈性支承無(wú)限長(zhǎng)梁模型。車(chē)輛-軌道-道床動(dòng)力學(xué)仿真模型如圖2所示。整體動(dòng)力學(xué)模型適合求解低頻條件下隧道的動(dòng)力響應(yīng)。

2.3 公路載荷

根據(jù)文獻(xiàn)［17］，公路車(chē)道載荷為均布載荷外加一個(gè)移動(dòng)載荷。均布載荷標(biāo)準(zhǔn)值為q=10.5 kN/m，公路為單向三車(chē)道，取車(chē)道載荷橫向分布系數(shù)為0.78。同時(shí)，根據(jù)規(guī)范和隧道限行規(guī)定，取公路移動(dòng)載荷為180 kN，其移動(dòng)速度為90 km/h，與軌道車(chē)輛中部同步前進(jìn)。

圖2 車(chē)輛-軌道-道床耦合模型Fig.2 Coupled train-rail-bed model

3 計(jì)算結(jié)果分析

首先根據(jù)施工步驟進(jìn)行靜應(yīng)力計(jì)算，采用釋放荷載法［18］，得到襯砌結(jié)構(gòu)和土層的最終的靜應(yīng)力場(chǎng)。在加載靜應(yīng)力場(chǎng)時(shí)，可能造成的殘余不平衡節(jié)點(diǎn)力會(huì)影響動(dòng)力計(jì)算結(jié)果。所以，計(jì)算隧道結(jié)構(gòu)動(dòng)力響應(yīng)之前，0～6s，加載靜應(yīng)力場(chǎng)與重力場(chǎng)，運(yùn)用動(dòng)力松弛法，使結(jié)構(gòu)達(dá)到平衡狀態(tài)。系統(tǒng)能量變化情況如圖3所示，動(dòng)能在較快的時(shí)間內(nèi)即衰減，單元節(jié)點(diǎn)力約在4秒后重新達(dá)到平衡，能量趨于穩(wěn)定。同時(shí)，系統(tǒng)一直處于臨界阻尼狀態(tài)，避免了過(guò)大的瞬態(tài)沖擊造成殘余塑性應(yīng)力。這表明，采用動(dòng)力松弛方法加載靜應(yīng)力場(chǎng)是非常有效的。

圖3 系統(tǒng)能量變化Fig.3 Energy curve of the system

6 s～7 s，以斜坡函數(shù)形式加載公路車(chē)道準(zhǔn)靜態(tài)載荷；8 s以后，列車(chē)由獨(dú)立于分析模型的外部軌道駛?cè)胨淼馈Ｍ瑫r(shí)，還分別計(jì)算了單獨(dú)考慮公路載荷和列車(chē)載荷的工況，分析它們的疊加效應(yīng)。

為了減小“邊界效應(yīng)”對(duì)結(jié)果分析的影響，選取隧道中部某普通斷面結(jié)果進(jìn)行分析。車(chē)輛通過(guò)時(shí)，斷面位移響應(yīng)如圖4所示：由公路車(chē)道均布載荷引起的襯砌垂向位移約1.15 mm，列車(chē)動(dòng)載引起的襯砌垂向位移峰值約為0.8 mm，公路移動(dòng)載荷引起的垂向位移很小，約0.2 mm。圖中還表明，公路車(chē)道均布載荷引起的隧道位移大于列車(chē)動(dòng)載的影響，而列車(chē)引起的動(dòng)態(tài)效應(yīng)更強(qiáng)；由公路載荷和列車(chē)動(dòng)載引起的振蕩峰值稍大于它們的簡(jiǎn)單疊加。

圖4 斷面位移響應(yīng)Fig.4 Displacement response of section

通過(guò)含靜應(yīng)力工況和無(wú)靜應(yīng)力工況的對(duì)比來(lái)分析靜應(yīng)力的影響。含靜應(yīng)力工況斷面襯砌外表面的法向應(yīng)力如圖5所示，底部單元接觸面的法向應(yīng)力時(shí)程曲線如圖6所示。列車(chē)經(jīng)過(guò)時(shí)，隧道襯砌的上下表面都承受了來(lái)自土體的較大壓力。然而，在不加靜應(yīng)力和土體重力的工況下，隧道襯砌上半部分沒(méi)有受到任何接觸壓力，相當(dāng)于改變了動(dòng)力系統(tǒng)的特性，列車(chē)經(jīng)過(guò)時(shí)土體對(duì)襯砌上半部分沒(méi)有約束。對(duì)比圖6與圖4可以發(fā)現(xiàn)，接觸面法向應(yīng)力時(shí)程曲線與襯砌的位移的時(shí)程曲線變化規(guī)律相似。因而，考慮靜應(yīng)力有利于使模型的接觸計(jì)算和隧道響應(yīng)計(jì)算更準(zhǔn)確。

圖5 接觸面法向應(yīng)力云圖Fig.5 Contour of contact pressure

圖6 襯砌底部法向應(yīng)力時(shí)程Fig.6 Time history of normal contact pressure on the bottom of tunnel lining

圖7 斷面環(huán)向應(yīng)力云圖Fig.7 Contours of hoop stress on choose section

在分析斷面應(yīng)力時(shí)，在該斷面處建立局部柱坐標(biāo)系，以襯砌圓心為原點(diǎn)，襯砌半徑方向?yàn)閺较颍瑘A周方向?yàn)榄h(huán)向，以隧道襯砌環(huán)的縱向?yàn)檩S向。兩種工況時(shí)的環(huán)向應(yīng)力云圖如圖7所示。含靜應(yīng)力工況下，襯砌底部?jī)?nèi)側(cè)受拉應(yīng)力，外側(cè)受壓應(yīng)力。而不含靜應(yīng)力工況下，底部?jī)?nèi)外側(cè)同受拉應(yīng)力，只是外側(cè)較大。

在列車(chē)經(jīng)過(guò)時(shí)，隧道襯砌的應(yīng)力分布規(guī)律較靜應(yīng)力場(chǎng)并沒(méi)有發(fā)生明顯的變化，只是在數(shù)值上產(chǎn)生了明顯的波動(dòng)，襯砌環(huán)底部環(huán)向應(yīng)力變化情況如圖8所示。列車(chē)通過(guò)時(shí)，隧道底部應(yīng)力振蕩幅度最大，底部外表面單元壓應(yīng)力振蕩幅值約0.1 MPa，無(wú)靜應(yīng)力工況下，底部外表面單元拉應(yīng)力振蕩幅值約0.12 MP；10節(jié)列車(chē)引起的振蕩次數(shù)為11次，相鄰兩懸架共同引起一次振蕩，首尾懸架各引起一次振蕩，且幅度較小；列車(chē)經(jīng)過(guò)時(shí)，環(huán)向應(yīng)力也包含一些小幅度的振蕩，這是由軌道不平度引起，顯然，由車(chē)輛所受重力引起的應(yīng)力振蕩要遠(yuǎn)遠(yuǎn)大于軌道不平度引起的應(yīng)力振蕩。

圖8 底部單元環(huán)向應(yīng)力時(shí)程Fig.8 Time history of hoop stress of the bottom element

綜上，隧道在列車(chē)和公路載荷作用下的位移響應(yīng)屬于很小量級(jí)，應(yīng)力響應(yīng)值也遠(yuǎn)小于靜應(yīng)力狀態(tài)，隧道在公路、鐵路同時(shí)運(yùn)營(yíng)的工況下是穩(wěn)定的。

4 并行計(jì)算效率分析

本文在曙光5000A超級(jí)計(jì)算平臺(tái)上，采用默認(rèn)的遞歸坐標(biāo)對(duì)分方法（RCB）和基于接觸負(fù)載均衡的分區(qū)方法（CBB），分別對(duì)仿真模型進(jìn)行了分區(qū)并求解。對(duì)于不同的CPU個(gè)數(shù)設(shè)置，采用加速比Sp來(lái)衡量并行計(jì)算的效率，其定義如下：

式中，ts為采用單個(gè)CPU的計(jì)算時(shí)間，tm為多CPU并行計(jì)算時(shí)所需時(shí)間。

圖9為不同分區(qū)方法對(duì)應(yīng)的拓?fù)湟晥D，可以發(fā)現(xiàn)，CBB方法的接觸分布更加合理，且網(wǎng)格之間交匯面較少，減少了信息傳遞開(kāi)銷(xiāo)。圖10為不同CPU個(gè)數(shù)時(shí)兩種分區(qū)方法的加速比情況。縱向比較：當(dāng)CPU個(gè)數(shù)為16個(gè)時(shí)，加速比較CPU數(shù)量較少時(shí)有明顯提升，而與32個(gè)CPU時(shí)差距不大，計(jì)算效率比較理想。這是由于隧道模型在空間布局上為細(xì)長(zhǎng)型，當(dāng)CPU個(gè)數(shù)大于16時(shí)，導(dǎo)致了不同處理器之間信息傳遞的開(kāi)銷(xiāo)迅速增大。橫向比較：當(dāng)CPU數(shù)量較多時(shí)，接觸負(fù)載均衡方法比默認(rèn)方法有約15%的提升，因而較適合于求解這類(lèi)大規(guī)模、含較多接觸計(jì)算的有限元問(wèn)題。

圖9 32分區(qū)拓?fù)鋱DFig.9 Topological diagram of 32 partitions

圖10 兩種分區(qū)方法加速比Fig.10 Speedups of the two kinds of decompositions

5 結(jié)論

（1）靜應(yīng)力場(chǎng)存在的條件下，可以更準(zhǔn)確的計(jì)算隧道與土體的接觸。采用隨時(shí)步變阻尼的動(dòng)力松弛方法加載靜應(yīng)力場(chǎng)，可以取得較好的效果。該方法也可應(yīng)用于其他含有靜應(yīng)力加載的相關(guān)工程應(yīng)用中。

（2）由列車(chē)重力及車(chē)道載荷引起的隧道襯砌的位移響應(yīng)大于軌道不平度的影響。拱頂及拱底的垂向位移比拱腰位置略大，隧道襯砌位移以及底部環(huán)向應(yīng)力的振蕩次數(shù)與列車(chē)編組數(shù)相關(guān)。

（3）公路車(chē)道均布載荷引起的隧道位移大于列車(chē)動(dòng)載引起的隧道位移，而列車(chē)載荷的動(dòng)態(tài)效應(yīng)則更強(qiáng)。對(duì)于非線性系統(tǒng)，公路載荷和列車(chē)動(dòng)載的綜合效應(yīng)不能通過(guò)二者響應(yīng)的簡(jiǎn)單疊加得到。

（4）采用基于負(fù)載均衡的分區(qū)并行計(jì)算方法，可以明顯提高含有大量接觸有限元模型的求解速度。同時(shí)，需要綜合考慮模型的規(guī)模和空間拓?fù)浣Y(jié)構(gòu)來(lái)決定CPU的個(gè)數(shù)，以取得更好的并行效率。

［1］王祥秋，楊林德，高文華.鐵路隧道提速列車(chē)振動(dòng)測(cè)試與荷載模擬［J］.振動(dòng)與沖擊，2005，24（3）：99-102.

［2］王祥秋，楊林德，周治國(guó).列車(chē)振動(dòng)荷載作用下隧道襯砌結(jié)構(gòu)動(dòng)力響應(yīng)分析［J］.巖石力學(xué)與工程學(xué)報(bào)，2006，25（7）：1337-1342.

［3］王祥秋，張玉紅，楊林德.珠三角地區(qū)典型地層地鐵區(qū)間隧道動(dòng)力特性分析［J］.振動(dòng)與沖擊，2008，27（1）：151-154.

［4］白 冰，李春峰.地鐵列車(chē)振動(dòng)作用下近距離平行隧道的彈塑性動(dòng)力響應(yīng)［J］.巖土力學(xué)，2009，30（1）：123-127.

［5］龔 倫，鄭余朝，仇文革.列車(chē)動(dòng)載引起下穿隧道振動(dòng)三維數(shù)值分析［J］.現(xiàn)代隧道技術(shù)，2008，45（4）：23-27.

［6］ Gupta S，Vanden B H，Lombaert G，et al.Numerical modelling of vibrations from a thalys high speed train in the Groene Hart tunnel［J］.Soil Dynamics and Earthquake Engineering，2010，30：82-97.

［7］Degrande G，Clouteaub D，Othmanb R，et al.A numerical model for ground-borne vibrations from underground railway traffic based on a periodic finite element boundary element formulation［J］.Journal of Sound and Vibration，2006，293：645-666.

［8］Andersena L，Jones C J.Coupled boundary and finite element analysis of vibration from railway tunnels-a comparison of twoand three-dimensional models［J］.Journal of Sound and Vibration，2006，293：611-625.

［9］Lombaert G，Degrande G.Ground-borne vibration due to static and dynamic axle loads of InterCity and high-speed trains［J］.Journal of Sound and Vibration，2009，319：1036-1066.

［10］ Belytschko T，Liu W K，Moran B.Nonlinear finite elements for continua and structures［M］.Chichester：John Wiley ＆Sons Ltd，2000.

［11］ Li Z，Jin X L，Chen X D.Parallel numerical simulation with domain decomposition for the fluid-filled drum crash［J］.Proceedings of the Institution of Mechanical Engineering，Part C：JournalofMechanicalEngineeringScience， 2008，222（3）：319-328.

［12］Guo Y Z，Jin X L，Ding J H.Parallel computing for seismic response analysis of immersed tunnel with domain decomposition ［J］. Engineering Computations， 2007，24（2）：182-199.

［13］劉晶波，李 彬.三維黏彈性靜–動(dòng)力統(tǒng)一人工邊界［J］.中國(guó)科學(xué)（E輯），2005，35（9）：966-980.

［14］鄧亞虹.層狀自由場(chǎng)地固有頻率的求解方法、特性及應(yīng)用研究［D］.杭州：浙江大學(xué)，2007.

［15］ Livermore Software Technology Corporation.LS-DYNA keyword user’s manual 971［Z］.California：Livermore，2007.

［16］劉學(xué)毅，王 平.車(chē)輛-軌道-路基系統(tǒng)動(dòng)力學(xué)［M］.成都：西南交通大學(xué)出版社，2010.

［17］JTG D70-2004.公路橋涵設(shè)計(jì)通用規(guī)范［S］.北京：人民交通出版社，2004.

［18］潘昌實(shí).隧道力學(xué)數(shù)值方法［M］.北京：中國(guó)鐵道出版社，1995.