周 彪，謝雄耀，李永盛

（1.同濟(jì)大學(xué) 巖土及地下工程教育部重點(diǎn)實(shí)驗(yàn)室，上海200092；2.同濟(jì)大學(xué) 地下工程系，上海200092）

波數(shù)有限元法或2.5維有限元法是在假設(shè)實(shí)際的鐵路、地鐵、隧道結(jié)構(gòu)和土體地質(zhì)條件在長(zhǎng)度方向幾何形狀和材料性質(zhì)不變的前提下.將其在無限長(zhǎng)運(yùn)行方向用傅立葉變換等周期化方法將其分解成一系列沿該方向傳播的簡(jiǎn)諧波.在計(jì)算每個(gè)簡(jiǎn)諧波時(shí)，只需要對(duì)結(jié)構(gòu)－土體的橫截面進(jìn)行離散.從而在利用2.5維方法模擬列車荷載作用下軌道、隧道以及周邊環(huán)境振動(dòng)的研究中，既克服了二維平面模型不能模擬移動(dòng)荷載作用下周圍地基中波動(dòng)傳播的馬赫輻射效應(yīng)的瓶頸以及無法考慮列車運(yùn)行方向的荷載影響的缺點(diǎn)，其計(jì)算量又明顯低于3維有限元分析模型.因而被國(guó)內(nèi)外學(xué)者大量采用.

該方法是Hwang等[1]在研究地震行波作用下地下結(jié)構(gòu)動(dòng)響應(yīng)時(shí)提出的，之后被國(guó)內(nèi)外學(xué)者引入到交通荷載作用下地基振動(dòng)以及地下軌道交通環(huán)境振動(dòng)的研究中.同時(shí)該方法在基于有限元基礎(chǔ)之上，其邊界模擬上除了采用最初的粘滯透射邊界外，邊界元和無限單元法等技術(shù)也被引入并與之進(jìn)行耦合.其中 Yang和 Huang[2－4]用2.5 維有限元法分析了移動(dòng)荷載作用下地基的動(dòng)力響應(yīng)問題，并分析了高速列車在彈性地基上運(yùn)行的動(dòng)力響應(yīng)以及沿著鐵路線的空溝的隔振.邊學(xué)成等[5－6]采用2.5 維有限元結(jié)合邊界元法分析了列車和地鐵荷載作用下結(jié)構(gòu)－地基的動(dòng)力響應(yīng)問題.謝偉平等[7]指出列車運(yùn)行時(shí)的振動(dòng)問題主要是由移動(dòng)荷載的自振與移動(dòng)位置的變化引起的，并進(jìn)行了相關(guān)分析，應(yīng)用了2.5維有限元方法，對(duì)移動(dòng)方向進(jìn)行波數(shù)積分求解，并考慮了地基土不同性質(zhì)的影響.Sheng等[8]對(duì)波數(shù)有限元法和波數(shù)邊界元法應(yīng)用于列車運(yùn)行產(chǎn)生的大地振動(dòng)問題作了詳細(xì)的推導(dǎo).另外Gupta等[9]以及國(guó)內(nèi)劉維寧等[10]將Floquet技術(shù)應(yīng)用于有限元－邊界元（FEBE）耦合分析.

上述方法中除邊界元外，其他邊界處理方法均對(duì)有限元模擬區(qū)域及網(wǎng)格大小提出要求，Yang等[11]在利用有限元和無限元模擬二維振動(dòng)問題時(shí)分析了有限元模擬區(qū)域大小及網(wǎng)格質(zhì)量對(duì)模擬結(jié)果的影響，但未對(duì)使用2.5維技術(shù)過程中，上述參數(shù)的影響做出分析，因而本文基于2.5維方法，通過 Matlab編制程序，分別研究了粘滯阻尼邊界和無限元邊界條件下上述參數(shù)的影響，并對(duì)上述參數(shù)的合理選擇提出建議.

1 2.5維有限元方法

1.1 基本假設(shè)

如圖1所示，2.5維方法假設(shè)沿z方向以速度c移動(dòng)的荷載形式為

式中，Ψ（x，y）Φ（z－ct）反映荷載的空間特征，q（t）為其時(shí)間域特征，其可由諸多周期波組合而成.若假設(shè)q（t）＝exp（iω0t），經(jīng)由傅立葉變化其在頻域波數(shù)域內(nèi)的表達(dá)式為

ω0）／c，其中ω為角頻率，k為波數(shù).同理，可以通過式（3）所示逆變換將其重新轉(zhuǎn)化為時(shí)空域.

假設(shè)沿z方向地質(zhì)條件、幾何形狀和材料性質(zhì)均不變，且動(dòng)力響應(yīng)為線性，則在荷載Ψ（x，y）·exp（－ikz）exp（iωt）作用下，其頻域內(nèi)動(dòng)態(tài)響應(yīng)解可表示為D（ω）exp（－ikz）exp（iωt）.

由此，若可求取D（iω）值，則地基時(shí)域內(nèi)三維動(dòng)態(tài)響應(yīng)可通過式（4）所示變化求得.而其可以通過運(yùn)用數(shù)值方法對(duì)某一z值條件下的xoy截面進(jìn)行網(wǎng)格劃分求得，其具體方法將在節(jié)1.2中闡述.

圖1 2.5維模型示意圖Fig.1 The layout of the 2.5D method

1.2 波動(dòng)方程的2.5維數(shù)值表達(dá)式

根據(jù)1.1所述，在荷載 Ψ（x，y）exp（－ikz）·exp（iωt）作用下，將其位移響應(yīng)值定義為

根據(jù)Navier方程，在頻域內(nèi)的三維彈性土體介質(zhì)波動(dòng)方程可表示為

其中采用復(fù)系數(shù)的Lame常數(shù)λ′和μ′來考慮土體阻尼對(duì)波動(dòng)傳播的影響.

結(jié)合式（5）—（6），根據(jù)小變形假設(shè)，三維土體單元的幾何方程可表示為

同時(shí)根據(jù)應(yīng)力—應(yīng)變關(guān)系得到：

采用普通8節(jié)點(diǎn)等參單元對(duì)xoy平面進(jìn)行離散化處理，并運(yùn)用相關(guān)邊界條件，根據(jù)伽遼金（Galerkin）法即可得到在頻域中的離散方程：

式中，質(zhì)量矩陣M、剛度矩陣K和荷載矩陣F分別表示為

其中：

式中：ρ為單元材料密度；η及ζ為局部坐標(biāo)，上標(biāo)‘＊’表示共軛矩陣，上標(biāo)‘T’表示矩陣轉(zhuǎn)置.通過求取式（9），便可獲得U～，亦即D（iω）值.

2 邊界處理方法

在分析動(dòng)荷載作用下地基瞬態(tài)響應(yīng)過程中，由于有限元方法中網(wǎng)格劃分范圍以及單元大小是有限的，難以直接對(duì)無限地基進(jìn)行模擬，故必須建立一種合適的人工邊界來阻止外行波反射回有限元單元表示的近場(chǎng)區(qū)域中.在第一章所述邊界處理上，無限元和透射邊界都對(duì)有限元模擬區(qū)域或網(wǎng)格的大小有不同程度的要求.因而本文就各類邊界處理方法及其相關(guān)對(duì)近場(chǎng)有限元網(wǎng)格的要求做出闡述.

2.1 粘滯阻尼透射邊界

所謂粘滯阻尼透射邊界，是一種運(yùn)用與頻率相關(guān)的黏壺單元來模擬無限地基，如圖2所示，圖中，kx＝kz＝ρvs，ky＝ρvp，vs和vp分別為剪切波（S波）和壓縮波（P波）波速.盡管這種黏壺單元對(duì)外行波能量的吸收不但取決于地基土體的材料特性，同時(shí)也與外行波的頻率有關(guān)，近來的研究表明這種吸收單元在實(shí)際運(yùn)用中非常有效[6].

圖2 頻率相關(guān)黏壺透射邊界Fig.2 Frequency related sticky pot of transmitting boundary

該吸收單元的系數(shù)由半無限區(qū)域的材料特性來決定，在頻域內(nèi)，單元節(jié)點(diǎn)上的外加集中力與變形速率之間的關(guān)系可以表示為

2.2 無限元邊界

楊永斌[2]提出運(yùn)用2.5維FE－IFE耦合算法來模擬振動(dòng)波傳播，如圖3所示，其采用了較為精確的形函數(shù)，并且坐標(biāo)形函數(shù)與位移形函數(shù)不一致，坐標(biāo)形函數(shù)具體表達(dá)式如下：

位移形函數(shù)表示為

式中：α，k’分別為單元在局部坐標(biāo)系中的幅值衰減系數(shù)和名義波數(shù)，分別表示單元在局部坐標(biāo)系中的幅值衰減特性和相位滯后特性.楊永斌等[2]詳盡闡述了其數(shù)值的選擇，結(jié)合趙崇斌等[12]中關(guān)于中節(jié)點(diǎn)位置選擇原則可以定義上述參數(shù).其中波數(shù)k’的選擇有別于一般的二維無限元的取值，其值對(duì)應(yīng)于P波，S波，瑞利波（R波）的定義如下：

其整體坐標(biāo)與局部坐標(biāo)之間的轉(zhuǎn)換如圖3所示.

圖3 無限單元整體和局部坐標(biāo)系Fig.3 Global and local co－ordinates of the infinite element

按照有限元的伽遼金（Galerkin）法，最終也得推得如式（7）—（10）所示頻域－波數(shù)域中與2.5維有限元一致的形式，但由于無限元在ζ方向上積分范圍為[0，∞]，因而其積分方法不能采用常用的高斯積分準(zhǔn)則，張楚漢等[13]給出了相應(yīng)的積分方法，本文采用該方法對(duì)式（10）進(jìn)行積分得到相應(yīng)的質(zhì)量、剛度及荷載矩陣.

3 不同邊界條件對(duì)有限元網(wǎng)格選擇的影響

為了檢驗(yàn)文中所提方法在Matlab中實(shí)現(xiàn)的準(zhǔn)確性以及網(wǎng)格劃分質(zhì)量和范圍的影響，本文利用Ansys建立如圖4所示模型，圖中L為單元最大尺寸，R為荷載作用點(diǎn)距離有限元邊界距離，中間白色區(qū)域?yàn)?節(jié)點(diǎn)有限元單元網(wǎng)格，邊界灰色區(qū)域?yàn)闊o限元單元.如圖4，在模型中央施加移動(dòng)荷載f（x，y，z＝0，t），根據(jù)1.1節(jié)所述原理，ω 取為32 Hz，ω0取為0 Hz，移動(dòng)速度取為120 m·s－1.通過從 Ansys中輸出荷載、單元和節(jié)點(diǎn)信息并輸入到Matlab中，基于上述2.5維FE－IFE耦合方法再將有相關(guān)單元信息進(jìn)行重組，得到相應(yīng)的質(zhì)量、剛度及荷載矩陣，根據(jù)式（9）即可求取其在頻域范圍內(nèi)的振動(dòng)響應(yīng)解.

圖4 有限元和邊界元模型及網(wǎng)格劃分Fig.4 The modeling of the finite element method（FEM）and infinite element method（IFE）

圖4所示模型中材料彈模取為2×107MPa，泊松比取為0.25，密度為2 000 kg·m－1，阻尼比為0.05，通過改變模型尺寸和網(wǎng)格大小，計(jì)算地基土的動(dòng)力響應(yīng).選取如圖4所示距離地表0.9 m所示測(cè)線上頻域范圍內(nèi)土體三向位移解，并乘以2πG／c以求得量綱—值，其中G為剪切模量.將其結(jié)果與Eason所提出的解析解[14]進(jìn)行比較，以檢驗(yàn)網(wǎng)格質(zhì)量對(duì)模擬結(jié)果的影響.

3.1 不同邊界下有限元模擬區(qū)域尺寸的選擇

針對(duì)粘滯邊界和無限元邊界兩種情況，根據(jù)節(jié)2.1—2.2所述方法求解在簡(jiǎn)諧荷載作用下地基響應(yīng)問題.如圖5，當(dāng)有限元模擬區(qū)域R取為1.5倍瑞利波波長(zhǎng)時(shí)（根據(jù)材料參數(shù)，瑞利波波速約為4.48 m，以下均表示為λ），測(cè)線處粘滯邊界和無限元邊界的計(jì)算結(jié)果與Eason解析解的對(duì)比情況如圖5所示.

圖5 不同邊界對(duì)豎向位移解的影響Fig.5 Effect of different boundary condition on vertical displacement

從圖5中可以看出無限元邊界較粘滯邊界有更好的模擬效果.在此工況下，無限元解已經(jīng)與解析解基本吻合，而粘滯邊界仍有較大的偏差，特別是距離邊界處這種趨勢(shì)更加明顯.同時(shí)可以看出虛部解偏差更為顯著，所以下文分析中對(duì)某些工況只比較虛部解的差異來確定影響程度的大小.

為更好分析粘滯邊界條件下模擬尺寸的影響，如圖6所示，分別選取R為0.8λ，1.5λ，2λ以及4λ，比較測(cè)線豎向位移的虛部解的差異.從圖中可以看出所有曲線都能反映波動(dòng)的大致趨勢(shì).但小于1λ時(shí)偏差較大，當(dāng)模擬距離大于2λ時(shí)，其結(jié)果才與解析解較吻合，但仍有偏差，所以如運(yùn)用粘滯邊界條件時(shí)，應(yīng)將模擬區(qū)域大于2λ，如需精確解要將范圍擴(kuò)至3λ～4λ.

圖6 粘滯邊界條件下有限元模擬尺寸對(duì)豎向解的影響（虛部）Fig.6 Effect of mesh size R on vertical displacement（image part） on condition of using visco boundary

以上可以看出，在相同條件下運(yùn)用無限元邊界將有更好的效果.因而本文以下著重分析在無限元邊界條件下網(wǎng)格大小及模擬區(qū)域?qū)τ?jì)算結(jié)果的影響.

3.2 無限元邊界有限元模擬區(qū)域尺寸（R）的選擇

在運(yùn)用無限元邊界條件下，分別選取模擬范圍為0.5λ、0.8λ、1.5λ以及2λ，通過比較測(cè)線x，y，z方向上位移分析模擬尺寸的影響.

圖7分別描述了不同長(zhǎng)度模擬區(qū)域下豎向（y向）實(shí)虛部解.從結(jié)果上看，模擬區(qū)域的大小對(duì)虛部解影響較大.從結(jié)果上看當(dāng)模擬范圍取為0.8λ時(shí)，結(jié)果雖與解析解偏差已經(jīng)很小.所以當(dāng)僅關(guān)注作用方向上位移時(shí)，模擬區(qū)域可超過1λ即可滿足精度要求.

圖8—9分別顯示了x方向和z方向上位移虛部解與解析解的差異，從結(jié)果上看，模擬區(qū)域?qū)ζ溆绊懸黠@大于豎向位移.特別是z方向即荷載移動(dòng)方向上，當(dāng)L＜λ時(shí)，其偏差較大，僅當(dāng)模擬范圍超過1.5λ時(shí)才發(fā)現(xiàn)數(shù)值解收斂并與解析解吻合.同時(shí)發(fā)現(xiàn)x方向上邊界有反射現(xiàn)象發(fā)生，這主要是因?yàn)樵跓o限元邊界參數(shù)中波數(shù)參量的選擇只能考慮一種波，而實(shí)際上地基土中存在有P，S以及R波，從而會(huì)對(duì)結(jié)果造成影響.因而在模擬x和z向土體響應(yīng)時(shí)，其模擬范圍應(yīng)該超過1.5λ.

圖7 無限元邊界條件下有限元模擬尺寸R對(duì)豎向解的影響Fig.7 Effect of mesh size R on vertical displacement on condition of using infinite element boundary

圖8 無限元邊界條件下有限元模擬尺寸R對(duì)x向解的影響（虛部）Fig.8 Effect of mesh size R on x direction displacement（image part）on condition of using infinite element boundary

圖9 無限元邊界條件下有限元模擬尺寸R對(duì)z向解的影響（虛部）Fig.9 Effect of mesh size R on z direction displacement（image part）on condition of using infinite element boundary

3.3 無窮元邊界有限元網(wǎng)格大小的選擇

從以上分析可以看出，當(dāng)模擬范圍為2λ時(shí)，數(shù)值解能滿足精度要求，因而本文在模擬范圍為2λ前提下，對(duì)不同網(wǎng)格大小的影響做出分析.

圖10 無限元邊界條件下有限元網(wǎng)格大小L對(duì)y向解的影響（虛部）Fig.10 Effect of element size L on y direction displacement（image part）on condition of using infinite element boundary

圖11 無限元邊界條件下有限元網(wǎng)格大小L對(duì)x向解的影響（虛部）Fig.11 Effect of element size L on x direction displacement（image part）on condition of using infinite element boundary

圖12 無限元邊界條件下有限元網(wǎng)格大小L對(duì)z向解的影響（虛部）Fig.12 Effect of element size L on z direction displacement（image part）on condition of using infinite element boundary

圖10—12分別為x，y，z向位移的虛部解，從結(jié)果上看，網(wǎng)格大小的影響要小于模擬范圍，其影響大多集中在荷載作用點(diǎn)附近.從圖中結(jié)果可以得出結(jié)論，當(dāng)網(wǎng)格大小超過λ／10時(shí)即能滿足精度要求.

4 結(jié)論

本文基于Matlab平臺(tái)，編寫了適用于計(jì)算和預(yù)測(cè)列車振動(dòng)波傳播和衰減的2.5維計(jì)算模塊，討論了粘滯邊界以及無限元邊界條件在不同有限元模擬范圍以及網(wǎng)格大小對(duì)結(jié)果的影響，并得出以下結(jié)論：

（1）在運(yùn)用2.5維計(jì)算方法模擬地基響應(yīng)過程中，無限元邊界較粘滯邊界有更好的適用性，其所需要的有限元近場(chǎng)模擬區(qū)域更小.

（2）在運(yùn)用無限元模擬地基響應(yīng)時(shí)，為獲取高精度解，其有限元近場(chǎng)模擬區(qū)域邊緣距離荷載作用點(diǎn)必須超過某一數(shù)值，對(duì)于豎向響應(yīng)，其必須大于1倍瑞利波波長(zhǎng).研究橫向和側(cè)向響應(yīng)，應(yīng)大于1.5倍瑞利波波長(zhǎng).

（3）在模擬過程中，網(wǎng)格大小對(duì)計(jì)算結(jié)果精度影響較小，當(dāng)網(wǎng)格大小大于0.1倍瑞利波波長(zhǎng)時(shí)即可滿足計(jì)算精度要求.

[1] Hwang R N，Lysmer J.Response of buried structures to traveling waves [J].Journal of Geotechnical Engineering，ASCE，1981，107（GT2）：183.

[2] Yang Y B，Hung H H.A 2.5D finite／infinite element approach for modelling visco－elastic bodies subjected to moving loads[J].International Journal for Numerical Methods in Engineering，2001，51：1317.

[3] Yang Y B，Hung H H.Chang D W.Train－induced wave propagation in layered soils using finite／infinite element simulation[J].Soil Dynamic and Earthquake Engineering，2003，23（4）：263.

[4] Hung H H，Yang Y B，Chang D W.Wave barriers for reduction of train－induced vibrations in soils[J].Journal of Geotech and Geoenvironment Engineering，ASCE，2004，130（12）：1283.

[5] 邊學(xué)成，陳云敏.基于2.5維有限元方法分析列車荷載產(chǎn)生的地基波動(dòng) [J].巖石力學(xué)與工程學(xué)報(bào)，2006，25（11）：2335.BIAN Xuecheng，CHEN Yunmin.Ground vibration generated by train moving loadings using 2.5D finite element method[J].Chinese Journal of Rock Mechanics and Engineering，2006，25（11）：2335.

[6] 邊學(xué)成，陳云敏，胡婷.基于2.5維有限元方法模擬高速列車荷載產(chǎn)生的地基振動(dòng)[J].中國(guó)科學(xué)，2008，38（5）：600.BIAN Xuecheng， CHEN Yunmin， HU Ting. Numerical simulation of high－speed train induced ground vibrations using 2.5D finite element approach[J].Science in China Press（G），2008，38（5）：600.

[7] 謝偉平，孫洪剛.地鐵運(yùn)行時(shí)引起的土的波動(dòng)分析 [J].巖石力學(xué)與工程學(xué)報(bào)，2003，22（7）：1180.XIE Weiping，SUN Honggang.FEM analysis on wave propagation in soils induced by high speed train load [J].Chinese Journal of Rock Mechanics and Engineering，2003，22（7）：1180.

[8] Sheng X，Jones C J C，Thompson D J.Modelling ground vibration from railways using wave number finite－and boundary element methods[J].Proceedings of the Royal Society A，2005，461：2043.

[9] Gupta S，Degrande G，Lombaert G..Experimental validation of a numerical model for subway induced vibrations [J]，Journal of Sound and Vibration，2009，321（3－5）：786.

[10] 劉衛(wèi)豐，劉維寧 地鐵振動(dòng)預(yù)測(cè)的周期性有限元－邊界元耦合模型[J].振動(dòng)工程學(xué)報(bào)2009，22（5）：480.LIU Weifeng，LIU Wenning.A coupled periodic finite elementboundary element model for prediction of vibrations induced by metro traffic[J].Journal of Vibration Engineering，2009，22（5）：480.

[11] Yang Y B，Kuo S R，Hung H H.Frequency－independent infinite element for analyzing semi－infinite problems[J].International Journal for Numerical Methods in Engineering，1996，39：3553.

[12] 趙崇斌，張楚漢.映射動(dòng)力無窮元及其特性研究[J].地震工程與工程動(dòng)，1987，3（22）：1.ZHAO Congbin，ZHANG Chuhan.Studies on the characteristics of dynamic mapping infinite elements [J].Journal of Earthquake Engineering and Engineering Vibration，1987，3（22）：1.

[13] 張楚漢，趙崇斌.復(fù)雜地基中波動(dòng)問題的數(shù)值模擬[J].土木工程學(xué)報(bào)，1987，20（4）：83.ZHANG Chuhan，ZHANG Congbin.A numerical simulation of wave propagation for complex foundation [J].China Civil Engineering Journal，1987，20（4）：83.

[14] Eason G.The stresses produced in a semi－infinite solid by a moving surface force[J].International Journal of Engineering Science，1965，2：581.