非飽和土-結構動力響應的多耦合周期性有限元法1)

狄宏規(guī) *, 郭慧吉 , **, 周順華 *, 王炳龍 *, 何 超 *,

* (同濟大學道路與交通工程教育部重點實驗室，上海 201804)

? (同濟大學,上海市軌道交通結構耐久與系統(tǒng)安全重點實驗室,上海 201804)

** (上海申通地鐵集團有限公司技術中心,上海 201103)

引言

近年來,地鐵隧道車致環(huán)境振動問題日益突出,嚴重影響著線路周邊居民的正常生活、精密儀器的正常使用及古建筑與文物的保護等[1-3].為采取有效的減(隔)振措施,迫切需要可靠的車致動力響應預測模型.

常見的地鐵隧道動力響應預測模型主要包括解析(半解析)模型以及數(shù)值法模型[4].解析(半解析)模型主要有歐拉梁-地基模型[5-7]、Pip in Pip(PiP)模型[8-10]及波函數(shù)法模型[11-12]等.解析(半解析)模型計算效率較高,但難以處理不規(guī)則邊界.與解析(半解析)模型相比,數(shù)值模型的優(yōu)勢在于可以進行更為精細化的建模,便于處理復雜的邊界(界面),可考慮復雜的地基結構.常見的數(shù)值模型有二維有限元模型[13-14]、三維有限元模型[15-17].二維有限元忽略了彈性波沿隧道縱向的傳播、三維有限元則占用大量的計算機資源,為兼顧二維有限元計算效率,同時考慮隧道縱向的動力傳播規(guī)律,相關性學者假定模型沿縱向幾何和物理性質(zhì)的不變性,基于縱向Fourier或Floquet 變換,分別提出了2.5 維有限元[18-19]和周期性有限元模型[20-21].此外,為避免人工截斷模型邊界產(chǎn)生的誤差,相關學者提出了邊界元[22-23]、無限元[24-25]、比例邊界有限元[26]、多阻尼層[27]等邊界處理方法.然而上述既有的數(shù)值模型僅將土體視為單相彈性或兩相飽和多孔介質(zhì),鮮見考慮地基土體的三相特性,同時無法高效考慮模型的縱向變化特性.

為此本文基于非飽和土波動方程以及滲流運動方程,結合應力、滲流邊界條件,采用 Galerkin 法,推導了非飽和ub-pl-pg格式有限元方程.同時引入拉伸函數(shù)構建完美匹配層(PML)處理邊界,提出了非飽和地基-結構動力響應分析的有限元-完美匹配層方法.在此基礎上進一步引入自由波傳播理論,將結構沿一個方向的變化處理為多耦合周期性結構,基于自由波傳播常數(shù)及傳播向量進行各周期性結構間的耦合,最終提出了多耦合周期性有限元法,并進行了地基-隧道-排樁系統(tǒng)動力響應的算例分析.

1 多耦合周期性有限元法

1.1 非飽和ub-pl-pg格式有限元推導

為實現(xiàn)頻域內(nèi)的模型求解,首先定義關于時間變量t的傅里葉變化形式

式中,ω代表時間t對應的頻率;頂標“～ ”表示頻域中的變量.

基于式(1)中傅里葉變化形式,根據(jù)三相介質(zhì)理論[28],可得頻域內(nèi)三相介質(zhì)波動控制方程及滲流方程分別為

式中,B11～B25表達式詳見附錄A.其余各變量定義可參考文獻[28].

整理式(3)可得

引入應力、滲流邊界條件

式中,f,cl,cg分別表示邊界節(jié)點處外力及輸入的液體、氣體流速.

基于Galerkin 法,在式(2)和式(5)的基礎上,引入虛位移δub以及虛孔壓δpl,δpg,同時利用式(4)在計算過程中消去v,w分量,可得ub-pl-pg格式的虛原理表達式為

式中,各變量表達式詳見附錄B.

采用空間8 節(jié)點等參單元進行離散處理,可得ub-pl-pg格式的三相介質(zhì)有限元表達式為

式中,各變量表達式詳見附錄C.

1.2 完美匹配層單元

采用完美匹配層處理模型邊界,可在有限的單元內(nèi)實現(xiàn)彈性波無反射的迅速衰減,進而截斷無限域.因此引入拉伸函數(shù)[29],推導三相非飽和多孔介質(zhì)完美匹配層單元

式中,頂標“?”表示拉伸含義.

彈性波在傳播過程中可分為兩種波,即傳播波(propagating)與隱逝波(evanescent),因此定義拉伸函數(shù)為

式中,a0為無量綱頻率常數(shù),a0=ωLPML/c;LPML為PML 區(qū)間長度;ω為角頻率;c為彈性波波速.

為驗證完美匹配層對流體中彈性波衰減的有效性,以孔隙流體影響的一維壓縮波的衰減為例,假定壓縮波傳播中位移函數(shù)為eik(k為壓縮波波數(shù)),計算從20 m 處開始坐標拉伸處理時,其位移幅值隨傳播距離的變化,并與無坐標拉伸時計算結果進行對比,如圖1 所示.可以發(fā)現(xiàn),經(jīng)坐標拉伸變化后,壓縮波迅速衰減,且對無拉伸處理區(qū)域內(nèi)彈性波傳播無影響.由此可見,合理選取坐標拉伸參數(shù)構建完美匹配層,同樣可實現(xiàn)氣相及液相中彈性波的無反射衰減.

圖1 一維壓縮波位移衰減對比圖Fig.1 The comparison diagram of displacement attenuation of onedimensional compression wave

將式(9)代入式(7)的推導過程中,可得完美匹配層有限元表達式為

1.3 多耦合周期性有限元模型

將結構根據(jù)縱向性質(zhì)的差異劃分為若干個周期性結構,如圖2 中的周期性結構1,2,j等.同時采用半無限周期性結構模擬縱向邊界.每個周期性結構由若干相同的單元組成.為便于模型耦合求解,將載荷作用結構j劃分為有載荷作用區(qū)域(jB),與無載荷作用區(qū)域(jA,jC).

圖2 多耦合周期性結構Fig.2 Multi coupling periodic structure

整理式(7),可得控制方程的矩陣表達式為

式中

對于施加外力的周期性區(qū)域(jB),例如圖2 中單元jB1,jB2,其控制方程矩陣表達式可寫為

沒有施加外力的周期性單元,例如圖2 中單元jAi,其單元控制方程矩陣表達式可寫為

根據(jù)自由波傳播理論[30],當彈性波經(jīng)過一個周期性單元后(圖3),其前后的位移、應力、孔壓等服從指數(shù)函數(shù)的傳遞規(guī)律,即

式中,?為自由波傳播常數(shù).

由圖3 可知,相鄰單元間存在如下連續(xù)性條件

圖3 自由波傳播示意圖Fig.3 Schematic diagram of free wave propagation

將式(14)和式(15)代入式(13),可得

進一步整理可得

觀察式(17),其為一元高次方程組,對其進行求解可得2m(m為單元左面(或右面)的自由度)個?.由自由波傳播理論可知,當?實數(shù)部為正數(shù)時,其為由單元左面向單元右面?zhèn)鞑サ恼蜃杂刹?當?實數(shù)部為負數(shù)時,其為由單元右面向單元左面?zhèn)鞑サ呢撓蜃杂刹?

求解式(17)時,每個?對應一個特征向量即φ,將自由波傳播常數(shù)?及其特征向量φ代入式(13),可得

式中,φ,ψ分別為Z,L對應的特征向量;p,n分別代表正向傳播波與負向傳播波分量.

而單元的位移及孔壓可表示為各自由波特征向量的線性組合,即

式中,α為特征向量的線性組合系數(shù).

根據(jù)式(19)推導各部分周期性結構剛度矩陣,主要分為4 類進行討論,即左邊界、右邊界、無載荷作用的周期性結構以及載荷作用的周期性結構.

(1)模型左邊界

由于模型左邊界左側為半無限空間,因此不存在波的反射和折射面,其只存在負向自由波,故

由式(20)整理可得

(2)模型右邊界

同理,由于模型右邊界右側為半無限空間,故

由式(22)整理可得

(3)無載荷作用的周期性結構

假定周期性結構j上無外載荷作用,如圖4 所示,根據(jù)式(19)整理可得

圖4 周期性結構自由波傳播示意圖Fig.4 Schematic diagram of free wave propagation in periodic structure

故周期性結構j左右應力、位移及孔壓關系矩陣為

(4)載荷作用周期性結構

如圖2 所示,載荷作用周期性結構可分為A,B,C,3 個區(qū)域,其中A,C 兩區(qū)參考無載荷作用的周期性結構處理方式即可,對于載荷作用的B 區(qū)采用傳統(tǒng)3 維有限元進行分析,即結合式(12),式(21),式(23)以及式(25)可得模型的總體剛度矩陣為

由Gaussian 消元法,可得

將式(27)代入式(26),化簡得

同理,可解得任意周期性結構左右界面位移、孔壓值,將之代入式(24)可解的各周期性結構模態(tài)坐標值,最終由式(24)可解得頻域內(nèi)周期性結構內(nèi)部單元的位移、孔壓值.在此基礎上進行Fourier 逆變化,可得時域中的解.

當載荷為移動載荷時

式中,v0為載荷移動速度,ω0為載荷激振角頻率

通過離散的Fourier 變化即可實現(xiàn)移動載荷作用下的模型求解.

2 模型的驗證

為了驗證模型的可靠性,首先采用將本文模型中非飽和土參數(shù)退化的方式,來計算單相彈性地基土-隧道系統(tǒng)的動力響應,并與文獻[31]中單相彈性地基土-隧道2.5 D FE-PML 模型的計算結果進行對比.2.5D FE-PML 模型長為14 m,寬為17 m,PML邊界厚度h=1 m.多耦合周期性有限元模型橫斷面尺寸與2.5 D FE-PML 模型相同,縱向一個周期單元長度為0.2 m,模型簡化示意如圖5 所示.模型計算參數(shù)及參數(shù)退化方式參考文獻[12].

圖5 驗證模型示意圖Fig.5 Schematic diagram of validation model

圖6 給出了隧道仰拱處作用單位簡諧載荷(f=30 Hz,v0=0 m/s)時,不同埋深處豎向動位移ux沿水平向變化曲線.由圖6 可得,兩者計算結果吻合較好,一定程度上驗證了本文提出的多耦合周期性有限元方法的可靠性.

圖6 多耦合周期性有限元方法與2.5 D FE-PML 模型計算結果豎向動位移ux對比圖Fig.6 Comparison of vertical dynamic displacement uxof multi coupling periodic finite element model and 2.5 D FE-PML model

為了進一步驗證本文方法中對于非飽和土考慮的正確性,將本文方法的計算結果與既有的解析算法[28]的計算結果進行對比,計算參數(shù)參考表1.

表1 驗證模型計算參數(shù)Table 1 Calculation parameters of validation model

圖7 給出了兩個非飽和地基-隧道模型在隧道仰拱處(x=0 m,y=-2.75 m,z=0 m)作用單位簡諧載荷(f=10 Hz,20 Hz,30 Hz,40 Hz)時,豎向動位移沿y=-3 m 的分布規(guī)律.從圖7 可以看出,兩者結果較為一致,進一步驗證了本文提出的多耦合周期性有限元方法的可靠性.

圖7 多耦合周期性有限元方法與解析算法計算結果對比圖Fig.7 Comparison of calculation results of multi coupling periodic finite element model and analytical algorithm

需要指出的是,本文方法計算時主要耗時為自由波特征值的求解,以驗證算例二為例,一個周期性結構單個面的自由度為5158,計算自由波傳播常數(shù)耗時7065 s,后續(xù)耗時為203 s.總的來說,計算效率比解析方法略低,但與有限元-邊界元方法相比,該方法在計算效率上具有較大的優(yōu)勢.

3 數(shù)值計算與分析

基于上述方法,以地基、隧道、隔離樁系統(tǒng)為例,討論隔離樁對地表振動響應的影響,在隧道仰拱處作用固定單位簡諧載荷(f=20 Hz),如圖8 所示.地基、隧道模型計算參數(shù)如表1 所示,隔離樁計算參數(shù)如表2 所示.

圖8 模型計算示意圖Fig.8 Schematic diagram of model calculation

表2 隔離樁計算參數(shù)Table 2 The calculation parameters of isolation pile

圖9 給出了無隔離樁、單排咬合隔離樁以及單排不咬合隔離樁(存在樁間距)情況下地表振動響應分布特征.可以發(fā)現(xiàn),隔離樁的存在使得地表振動響應規(guī)律發(fā)生了明顯變化,這主要是由于隔離樁的存在使得地基中產(chǎn)生了新的波反射(折射)面.靠近隔離樁附近,地表動力響應有所衰減,但距離隔離樁較遠的局部位置存在振動響應放大的現(xiàn)象.還可以看到,與單排不咬合隔離樁相比,咬合樁的隔離效果更好.由此可見,地基中結構的縱向多耦合周期性會影響系統(tǒng)動力響應預測結果,因此為準確預測系統(tǒng)動力響應,應考慮地基結構的縱向變化特性.

圖9 不同工況下地表振動響應分布Fig.9 The distribution characteristics of surface vibration response of double row pile wall

4 結論

(1)基于非飽和土運動微分方程、連續(xù)性方程以及滲流運動方程,結合應力、滲流邊界等條件,采用 Galerkin 法,推導了ub-pl-pg格式的固、液、氣三相非飽和介質(zhì)有限元表達式.相比于ub-v-w格式有限元單節(jié)點的9 個自由度,ub-pl-pg格式的有限元單節(jié)點減少至5 個自由度,大大節(jié)省了計算資源.

(2)基于自由波傳播理論,求解自由波傳播常數(shù)及自由波特征向量,基于自由波傳播常數(shù)及自由波特征向量求解各周期性結構剛度矩陣,在此基礎上結合傳統(tǒng)有限元、完美匹配層單元,提出了多耦合周期性有限元法,將該方法分別與2.5 維有限元法及解析算法計算結果相對比,驗證了該算法的可靠性.

(3)多耦合周期性有限元法具有可考慮結構沿縱向變化特性的優(yōu)點.其與解析方法相比,計算效率略低,但與有限元-邊界元方法相比,該方法在計算效率上具有一定優(yōu)勢.地基中結構的縱向變化特性會影響系統(tǒng)的動力響應,因此,進行動力響應預測時,應考慮地基中結構沿縱向的變化特性.

附錄

附錄A

附錄B

附錄C