高廣運，張繼嚴，謝 偉，徐晨曉

(1.同濟大學 地下建筑與工程系，上海 200092； 2.巖土及地下工程教育部重點實驗室(同濟大學)，上海 200092)

數(shù)值法可以計算分析復雜的動力模型，在高鐵引起的環(huán)境振動問題研究中得到廣泛應用。研究表明，2維模型忽略沿軌道方向的影響，計算結果不精確[1]；3維模型計算效率低[2]；2.5維模型假設沿軌道方向結構-土體的幾何形狀和材料性質(zhì)保持不變，對軌道方向和時間進行雙重Fourier變換，將空間-時間域的三維問題轉(zhuǎn)化為波數(shù)-頻域內(nèi)每個節(jié)點包含3個自由度的平面應變問題，只需選取垂直于軌道的結構-土體截面進行離散求解，再通過Fourier逆變換求得時域-空間域內(nèi)解答。2.5維有限元法兼顧精度和效率[3]。

2.5維有限元方法最早由Hwang等[4]用于地震行波作用下地下結構動力響應問題研究。楊永斌等[5]將該方法應用到列車運行引起的地面振動響應研究，提出2.5維FEM-IFEM，研究彈性半空間和分層地基的地面振動。之后，開發(fā)了包括多種地基類型和人工邊界條件的2.5維有限元動力計算模型[5]。Costa等[8]運用此法結合等效線性分析方法研究動荷載作用下土體的非線性特性并開發(fā)出2.5維FEM-BEM，建立輪軌相互作用耦合模型，研究分層路基上的動力響應。邊學成等[10]開發(fā)與頻率相關的黏壺單元吸收有限單元邊界上的外行波，發(fā)展了高鐵移動荷載作用下軌道和地基動力相互作用的2.5維有限元數(shù)值方法，模擬軌道和地面的振動及波的傳播，并探討了軌道不平順因素的影響。Francois等[12]提出完美匹配層方法(PML)，將其應用到2.5維有限元方法中，通過計算彈性半空間的格林函數(shù)位移驗證了2.5D PML能夠吸收所有不同入射角的波?；诖?，Lopes等[13]建立2.5D FEM-PML軌道-隧道-地基系統(tǒng)模型和振源附近的3D FEM建筑物模型，研究僅考慮地面土體-結構之間的相互作用，表明某些墊層剛度會放大建筑樓板的垂直振動。高廣運等[15]開發(fā)出2.5維黏彈性邊界條件，分析了高鐵運行時分層地基、飽和分層地基、橫觀各向同性地基、飽和橫觀各向同性地基和非飽和地基的振動響應及波的傳播特性。

綜上，基于2.5維有限元法的研究均針對地基的振動響應，將地基簡化為彈性體，忽略了地基因高鐵運行產(chǎn)生的塑性變形。由于真實土體組構的復雜性，尚沒有一個能真實描述土體塑性的本構模型，本文采用改進的摩爾-庫倫模型模擬地基土體。將問題視為材料非線性問題，引入切線剛度迭代法、后歐拉積分算法和一致切線模量算法，推導出適用于高鐵荷載下彈塑性地基振動與變形分析的2.5維有限元方法。數(shù)值計算中將彈性理論得到的土體位移作為試探位移，經(jīng)改進的摩爾-庫倫屈服準則判定屈服后，通過施加體荷載和更新土體剛度矩陣模擬彈塑性土體的變形。計算模型將軌道板視為鋪設在地基上的歐拉梁，采用黏彈性邊界模擬波在半無限空間中的傳播，列車荷載采用連續(xù)軸重荷載模型和經(jīng)英國軌道不平順譜修正的隨機激振荷載模型模擬。

1 彈塑性地基非線性理論

采用切線剛度迭代法表征材料的非線性。切線剛度迭代法采用定期更新總剛度矩陣和體荷載迭代相結合的方法獲得剛度方程的收斂，該方法能夠得到加載過程中的變形發(fā)展狀況，適用于動態(tài)問題求解，近年來得到廣泛應用。

1.1 體荷載生成

土體達到屈服條件后，通過迭代不斷修改加載在系統(tǒng)上的荷載，并重復利用彈性求解方法獲得剛度方程的收斂。在每一個荷載增量步中，位移增量值δi由剛度方程(1)求得，即

Kδi=pi

(1)

在進行應力重分布時，修改式(1)中的荷載增量矢量pi。通常，該荷載增量矢量表達式為

(2)

彈塑性土體的應力增量-應變增量顯式關系式如下:

Δσ=(De-Dp)Δε

(3)

式中：Δσ為應力增量，Δε為應變增量，De為彈性矩陣，Dp為塑性矩陣。

假定應力狀態(tài)已處于彈塑性材料的破壞面上，隨后的應力改變能夠使應力狀態(tài)漂移到破壞面的另一個位置，且不超過原有破壞面，于是有

(4)

F為屈服函數(shù)。

假設應力改變僅由彈性應變分量引起，結合非關聯(lián)流動法則，應力增量為

(5)

Q為塑性勢函數(shù)。將式(5)代入式(4)得

(6)

(7)

式中f為塑性矩陣Dp的修正因子，由線性插值方法求得。

1.2 改進摩爾-庫倫屈服準則

摩爾-庫倫屈服準則物理意義明確、使用簡便，廣泛應用于巖土工程分析?，F(xiàn)有有限元法多采用簡化模型不能對土體介質(zhì)的塑性變形進行合理描述。因此，對能夠合理描述土體微小塑性變形的改進摩爾-庫倫模型進行2.5維有限元實現(xiàn)，以達到高鐵荷載作用下地基振動與變形的準確高效求解。本文規(guī)定拉應力為正，壓應力為負。由應力不變量表示的摩爾-庫倫屈服準則：

(8)

摩爾庫倫準則在π平面上為六邊形，存在6個棱角奇異點，造成屈服函數(shù)和勢函數(shù)不連續(xù)及梯度無法計算，導致奇異點附近數(shù)值計算繁瑣和收斂緩慢。采用分段函數(shù)描述K(θ)，旨在消除摩爾-庫倫屈服準則在π平面上的奇異點，使新屈服面在拐點處平滑連續(xù)。K(θ)的表達式如下[19]：

(9)

式中：θT為過渡角，取25°；A、B的表達式詳見文獻[19]。

在子午平面上，采用雙曲線近似模擬屈服面，消除頂點奇異性和子午面拐點不可導的問題。

修正后的摩爾-庫倫準則屈服函數(shù)表達式為[20]

(10)

式中參數(shù)m的取值范圍為[0,1]，通過調(diào)整m的大小，可以反映雙曲線對屈服函數(shù)的擬合程度和巖土介質(zhì)抗拉強度的大小。

取塑性勢函數(shù)和屈服函數(shù)的表示式一致[20]，即

(11)

式中ψ為膨脹角，K(ψ)的表達式與K(θ)相同。

令矢量a和b分別表示屈服函數(shù)F和塑性勢函數(shù)G關于應力的一階偏導數(shù)：

(12)

(13)

式中各分量的具體表達式詳見文獻[20]。

后向歐拉算法需通過塑性勢函數(shù)的二階偏導數(shù)對節(jié)點應力進行更新，以判斷土體是否達到屈服狀態(tài)。塑性勢函數(shù)的二階偏導數(shù)表達式為[20]

(14)

式中各分量的表達式詳見文獻[20]，當|θ|>25°時，K(θ)=A-Bsin 3θ，代入式(14)可消除分母中的cos 3θ項，方便數(shù)值計算。

1.3 彈塑性本構積分算法

彈塑性理論以屈服函數(shù)判斷土體所處狀態(tài)。當屈服函數(shù)F>0時，土體處于塑性狀態(tài)；當F<0時，土體處于彈性狀態(tài)。采用返回算法，進行初始彈性預測，設定增量步，判定土體應力狀態(tài)所處位置(屈服面內(nèi)或屈服面外)，對處于屈服面外的應力進行塑性修正，使之返回到屈服面上。

對屈服面外的應力進行塑性修正反映其真實的加載過程(偏應力隨增量步的增加而增加)。根據(jù)正交流動法則，等效塑性應變率可表示為

(15)

采用各向同性硬化法則，將黏聚力c看作等效塑性應變的函數(shù)，即

(16)

利用屈服函數(shù)的一致性條件結合式(15)、(16)整理為

(17)

對如圖1所示的彈性試探應力σB進行塑性修正。定義應力改變量計算公式為

圖1 后向歐拉積分算法迭代示意圖

Δσ=DeΔε-ΔλDeb=Δσe-ΔλDeb

(18)

式中Δσe為彈性應力增量。

將屈服函數(shù)F在B點進行一階泰勒展開得

(19)

采用后向歐拉積分算法，通過迭代將σB拉回屈服面。后向歐拉應力(C點應力)表示為

σC=σB-ΔλDeb

(20)

其初始估計值表示為

σC=σB-ΔλDebB

(21)

定義矢量r為當前應力σ與后向歐拉應力的差值，即

r=σ-(σB-ΔλDebC)=σ-σB+ΔλDebC

(22)

r=0表示應力返回到屈服面上。

將σB固定，對r進行泰勒展開得

(23)

令rn=0，則

(24)

式中I為單位矩陣。

屈服函數(shù)F在C點的一階泰勒展開式為

(25)

將式(23)代入式(25)，可得

(26)

一致切線模量法能有效避免土體從彈性突變?yōu)樗苄詴r可能引起的偽加載或偽卸載，保證計算程序的收斂性。

由式(21)可得標準隱式后向歐拉算法表達式

σ=σB-ΔλDeb

(27)

對式(27)求微分，可得

(28)

(29)

將式(28)代入式(29)并簡化得

(30)

將式(30)代入式(28)并簡化得

(31)

式中Dct為一致性剛度矩陣。當土體達屈服狀態(tài)后，通過不斷更新Dct即可實現(xiàn)地基塑性變形的求解。

2 彈塑性地基2.5維有限元模型構建

2.1 彈塑性地基2.5維有限元

假定列車運行方向為x正方向，時間為t。定義變量關于x，t的雙重Fourier變換為

(32)

式中：ω表示圓頻率，ξx表示在x方向上的波數(shù)，上標“～”和“-”分別表示波數(shù)域和頻域內(nèi)的變量。

由第1節(jié)知，彈塑性地基非線性振動與變形求解采用切線剛度法解決材料非線性問題，土體經(jīng)改進的摩爾-庫倫屈服準則判定屈服后，通過施加體荷載和更新土體剛度矩陣模擬彈塑性土體的塑性變形。因此，每一迭代步的計算均為彈性土體求解問題。

土體運動平衡微分方程為

(33)

式(33)結合幾何方程和物理方程，并對時間t作傅里葉變換，整理得頻域內(nèi)以位移表示的平衡微分方程

(34)

式中：λc和μc為考慮土體材料阻尼的拉梅常數(shù)，λc=λ(1+2βi)，μc=μ(1+2βi)，β為土體阻尼比，i為虛數(shù)單位。

頻域內(nèi)的應力邊界條件

(35)

(36)

對式(36)進行分部積分，運用高斯公式并對x方向進行傅里葉變換得頻域-波數(shù)域內(nèi)的控制方程

(37)

式中：δu*為虛位移，δε*為虛位移對應的虛應變，上標“*”表示共軛復數(shù)。

對幾何方程和物理方程作雙重傅里葉變換，得頻域-波數(shù)域內(nèi)的矩陣形式為

(38)

(39)

用4節(jié)點等參單元進行離散，形函數(shù)表達式如下：

(40)

對式(37)引入形函數(shù)式(40)，可得彈塑性土體2.5維有限元方程的矩陣形

(41)

2.2 無砟軌道模型

無砟軌道由鋼軌、軌枕、扣件、道床、道岔等部分組成，共同承擔列車輪軌的作用力。在列車荷載作用下，整個軌道系統(tǒng)可近似為整體變形，因此，可將軌道系統(tǒng)簡化為鋪設在地基上的歐拉梁。列車荷載等效為寬度3 m的均布荷載[21]。

軌道梁的動力方程在頻域-波數(shù)域內(nèi)可表達為

(42)

采用位移協(xié)調(diào)條件耦合式(41)和(42)得軌道-地基的控制方程。

2.3 黏彈性人工邊界

為在有限域內(nèi)實現(xiàn)半無限空間的準確求解，設置合適的人工邊界阻止外行波的反射。黏彈性人工邊界具有穩(wěn)定性好、易編程的優(yōu)點。在有限元中，黏彈性邊界是在邊界節(jié)點的每個方向上施加一端固定的彈簧-阻尼元件，以黏性阻尼的吸波作用和彈簧的剛性復原作用模擬半無限空間。該邊界對外行波的吸收取決于地基土體的材料特性和外行波的頻率[10]。吸收邊界節(jié)點的系數(shù)由半無限區(qū)域材料性質(zhì)決定，作用于邊界節(jié)點上的荷載表達式為

(43a)

作雙重傅里葉變換得

(43b)

2.4 列車荷載模型

列車運行引起的激振荷載包括準靜態(tài)荷載和隨機激振荷載兩部分。準靜態(tài)荷載由列車軸重引起，隨機激振荷載涉及列車懸掛體系、軌道不平順和地基的均勻性等方面。

2.4.1 準靜態(tài)列車荷載

1節(jié)列車由1個車廂、兩個轉(zhuǎn)向架和4個輪對等組成，和諧號CRH380AL動車組列車由2拖14動組成，具體參數(shù)如表1所示。假定列車以速度c沿x正方向運行，則頻域-波數(shù)域內(nèi)的連續(xù)軸重荷載表達式為

表1 CRH380AL型電力動車組參數(shù)

(44)

式中：δ(x)為Dirac函數(shù)，Pn1和Pn2分別表示第n節(jié)車體前輪對和后輪對軸重，Li表示車廂長度，an和bn分別為第n節(jié)車體輪軸間距。

2.4.2 修正列車荷載

準靜態(tài)連續(xù)軸重荷載僅考慮列車軸重，不能考慮列車的隨機激振部分。梁波等[21]根據(jù)英國鐵路技術中心的理論分析和試驗研究，提出一個與高、中、低頻相應的，反映軌道不平順、附加動荷載和軌面波形磨耗效應的激勵力模擬列車荷載，同時考慮了鋼軌、軌枕的分散傳遞因素，該激勵力模型包含了列車隨機激振的主要組分。由該激勵力表示的修正列車連續(xù)軸重荷載為

(45a)

式(45a)涵蓋了列車動荷載的一系列重要因素(車輛因素、軌道及軌下基礎因素等)，能較好地表征真實列車荷載。其中：

表2 英國軌道不平順管理值

對式(45a)作雙重傅里葉變換，得修正列車荷載在頻域-波數(shù)內(nèi)的表達式為

(45b)

式中：當ω>0時為“+”，ω<0時為“-”。

j=0,1,2,3,4

3 數(shù)值方法驗證

依托VS2013開發(fā)平臺，編制高鐵荷載下彈塑性地基2.5維有限元Fortran計算程序。計算模型如圖2所示。通過以下兩種方法驗證理論及計算模型的正確性：1)與移動點荷載作用下彈性半空間振動的解析解的退化驗證；2)與京滬線上高鐵運行引起的地面振動實測結果進行驗證。

3.1 彈性半空間解析解驗證

Eason[23]推導出移動點荷載fd(x,t)=P0δ(x-ct)作用下彈性半無限空間振動問題的解析解。建立基巖層上厚30 m、寬40 m的2.5維有限元計算模型，如圖2(b)，忽略路堤，假定荷載作用在有限元模型土體表面中心處，采用軸對稱模型，模型兩側(cè)分別設置軸對稱邊界和黏彈性人工邊，土體及荷載參數(shù)同文獻[23]。當彈塑性模型中黏聚力趨近于無窮大時，土體不會屈服進入塑性狀態(tài)，彈塑性土體退化為彈性土體。

圖2 彈塑性地基2.5維有限元模型

由圖3可知，本文算法得到的在荷載移動方向和豎向的振動幅值與解析解吻合較好，證明本文理論的正確性。

圖3 移動點荷載作用下(0,0,-1)處歸一化振動位移幅值

3.2 實測驗證

翟婉明等[24]實測了京滬高鐵蘇州段的地面振動。如圖4所示，選取近軌道處(1.7 m)和遠軌道處(25 m)測點進行對比驗證。

圖4 測點布置示意

建立2.5維有限元計算模型如圖2(b)所示，對于雙線路堤取一半結構建模分析，模型兩側(cè)分別設置軸對稱邊界和黏彈性人工邊界，模型底部為固定邊界。計算模型路堤高6.3 m、面寬4.3 m，地基寬37 m、深31 m。樁采用等代樁墻處理，列車荷載模型采用1拖7動，其表達式如式(45)，參數(shù)取值同文獻[21]。土體及樁體參數(shù)同文獻[24]。

測點1.7 m處的振動加速度實測結果及本文計算結果如圖5所示。可以看出，本文計算結果頻率成分稀疏，但能表征主要的頻率成分。近軌道處，正向加速度峰值計算結果約為實測結果2倍。這是因為本文采用軸對稱模型，默認另一半結構也有同樣的荷載施加，對應雙線軌道同時有高鐵同向運行的最不利情況，而實測時只有一條軌道上有列車運行。另外，近軌道處的負向振動受荷載施加方式的抑制作用而使其負向加速度峰值小于正向峰值。

圖5 測點1.7 m處加速度時程對比

距軌道中心25 m處的實測加速度值及本文計算結果如圖6所示。遠軌道處振動響應受荷載施加方式及振動疊加效應影響較小，計算結果接近實測結果。由圖6可知，正負加速度峰值和主頻特性與實測結果吻合較好?？傮w而言，本文建立的計算模型能準確描述地基的振動峰值及主頻特性，證明了計算模型應用的可靠性。

圖6 測點25 m處加速度時程對比

4 結 論

1)以彈性地基2.5維有限元法為基礎，地基視為彈塑性，以切線剛度法表征地基土體非線性，引入改進的摩爾-庫倫屈服準則，以彈性理論得到的位移為試探位移，通過定期更新總體剛度矩陣和體荷載的方法實現(xiàn)高鐵荷載下彈塑性地基變形求解。

2)通過與移動點荷載引起的彈性半空間振動響應解析解對比，驗證了本文理論的正確性。將修正的列車荷載模型應用到本文建立的彈塑性地基2.5維有限元模型中，構建高鐵荷載下地基振動與變形分析模型，通過與實測結果對比，驗證了本文模型應用的可靠性。