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

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

（1.機械系統(tǒng)與振動國家重點實驗室，上海 200240；2.上海交通大學(xué) 機械與動力工程學(xué)院，上海 200240）

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

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

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

1 理論與方法

1.1 顯式中心差分方法

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

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

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

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

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

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

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

1.2 接觸算法

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

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

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

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

1.3 基于接觸負載均衡的并行算法

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

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

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

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

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

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

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

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

2 計算模型

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

2.1 隧道－土體模型

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

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

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

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

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

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

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

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

2.2 車輛－軌道－道床動力耦合模型

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

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

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

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

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

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

2.3 公路載荷

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

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

3 計算結(jié)果分析

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

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

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

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

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

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

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

圖6 襯砌底部法向應(yīng)力時程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)力時，在該斷面處建立局部柱坐標(biāo)系，以襯砌圓心為原點，襯砌半徑方向為徑向，圓周方向為環(huán)向，以隧道襯砌環(huán)的縱向為軸向。兩種工況時的環(huán)向應(yīng)力云圖如圖7所示。含靜應(yīng)力工況下，襯砌底部內(nèi)側(cè)受拉應(yīng)力，外側(cè)受壓應(yīng)力。而不含靜應(yīng)力工況下，底部內(nèi)外側(cè)同受拉應(yīng)力，只是外側(cè)較大。

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

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

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

4 并行計算效率分析

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

式中，ts為采用單個CPU的計算時間，tm為多CPU并行計算時所需時間。

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

圖9 32分區(qū)拓撲圖Fig.9 Topological diagram of 32 partitions

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

5 結(jié)論

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

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

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

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

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

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

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

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

［5］龔 倫，鄭余朝，仇文革.列車動載引起下穿隧道振動三維數(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］劉晶波，李 彬.三維黏彈性靜–動力統(tǒng)一人工邊界［J］.中國科學(xué)（E輯），2005，35（9）：966-980.

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

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

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

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

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