郭慧吉 狄宏規(guī) 周順華 何 超 張小會

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

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

引言

地下列車運行時,輪軌沖擊振動會經由軌道、隧道傳至地基土體,弓起地基土和周邊建(構)筑物的振動[1].隨著地下鐵道的快速發(fā)展,地下列車運行弓起的環(huán)境振動問題和軟土隧道沉降問題日益凸顯[2-4].

Metrikine 等[5]將隧道視為埋置于彈性土體中的歐拉梁,建立了二維解析模型,探討了隧道及地表在恒定、移動簡諧載荷和隨機載荷作用下的動力響應.Forrest 和Hunt[6]采用薄壁圓柱殼模擬隧道,彈性圓土柱模擬地基,提出了隧道-土體系統(tǒng)動力響應三維半解析計算方法,即PiP 半解析法.Yuan 等[7]將隧道簡化為埋置在彈性半空間的中空圓柱體,基于平面波和柱面波轉換,獲得了彈性半空間埋置單洞隧道的動力響應閉合解.He 等[8]通過弓入波的平移公式,提出了單相彈性半空間中埋置雙洞隧道系統(tǒng)動力響應的基本解.數(shù)值法方面,Thiede 等[9]采用二維有限元法研究了地鐵列車振動響應問題.Bian 等[10]采用2.5 維有限元法研究了地鐵列車運行弓起的地基振動.Sheng 等[11]采用波數(shù)離散法,建立了2.5 維波數(shù)有限元-邊界元法,研究了地面和地下列車載荷弓起的振動響應問題.Degrande 等[12]基于Floquet 變換,提出了周期性有限元-邊界元法,并用于振動預測;Hung 和Yang[13]采用無限域模擬地基,提出了2.5 維有限元-無限元法,模擬地下列車運行弓起的地表振動.此外,針對富水地層的隧道,相關學者采用Biot 飽和多孔介質模擬飽和土體,對上述模型或方法進行了發(fā)展和改進:如飽和全空間的PiP 半解析法[14-15];飽和半空間的波面轉換閉合解[16-17]以及飽和半空間的2.5 維有限元-邊界元法[18-20]等.

然而,采用單相或兩相介質難以模擬地基土的實際情況,需考慮飽和度對動力響應的影響[21-22].徐明江[23]基于V-G 模型,建立了非飽和土的動力控制方程,求解了非飽和土地基與基礎的動力響應.高廣運等[24]、李偉華等[25]基于非飽和介質波動理論,采用有限元法研究了移動載荷作用下路基和地基的動力響應.郭鵬飛等[26]、章敏等[27]基于非飽和土的動力控制方程,分析了非飽和土中樁基的動力響應問題.但受計算方法限制,目前對于非飽和土-隧道系統(tǒng)動力響應的研究較少.劉洪波[28]等提出了簡諧載荷作用下非飽和土中深埋圓形隧洞系統(tǒng)動力響應的2 維解析解.考慮到2 維計算方法無法考慮系統(tǒng)的3維動力特性,狄宏規(guī)[29]等在文獻[15]的基礎上,進一步提出了非飽和土-隧道系統(tǒng)動力響應計算的3維PiP 半解析法.但上述2 維或3 維解[28-29]均假定隧道外的土體外徑為無窮大,即均為全空間算法,故無法獲得真實半空間地基表面的動力響應.

為研究非飽和半空間中埋置圓形隧道系統(tǒng)的動力響應問題,本文將圓形隧道視為Flügge 薄壁圓柱殼,隧道外的非飽和地基土簡化為由流、固、氣組成的三相介質,分別采用分離變量法求解圓柱殼的振動控制方程,Helmholtz 矢量分解定理求解非飽和土的波動方程,并基于隧-土交界面處及地表等邊界條件,利用平面波和柱面波的轉換進行耦合求解,提出了移動簡諧載荷下非飽和半空間地基中埋置圓形隧道動力響應計算的波函數(shù)法,基于該方法,探討了非飽和地基土-隧道系統(tǒng)的動力響應特征.

1 控制方程和求解

將隧道襯砌簡化為無限長的Flügge 薄壁圓柱殼,殼體由均質、各向同性、線彈性材料組成.將非飽和地基土簡化為由流、固、氣組成的三相介質.計算過程中包括圖1 所示的圓柱坐標系(r,θ,z)和直角坐標系(x,y,z)兩種坐標系.圖中,O1,O2,O3,O4為4 個觀測點:O1點位于載荷作用點正上方地表,O3點位于載荷作用點下方隧道底部,O2,O4點分別位于O1,O3左側20 m 處.

圖1 非飽和土-隧道系統(tǒng)及坐標系Fig.1 Unsaturated soil-tunnel system and coordinate

圖中:q代表激振載荷;v0代表載荷移動速度;d代表隧道中心線距離地表豎直距離;r1代表隧道襯砌半徑.

1.1 襯砌振動控制方程

隧道襯砌視為Flügge 薄壁圓柱殼,利用分離變量法對殼體控制方程進行求解,可得簡諧載荷作用下殼體平衡方程的矩陣表達式[6]

式中,上標“?”代表z方向所對應的波數(shù)域;上標“～”表示時間t所對應的頻域;m表示環(huán)向模態(tài);分別表示在頻域-波數(shù)域內位移和應力在單個環(huán)向模態(tài)數(shù)m下的物理量;式(1)各變量表達式詳見附錄A.

1.2 非飽和地基土波動方程

將非飽和地基土視為流、固、氣組成的三相介質,其實用波動方程如下[23]

式中,a=1-Kb/Ks,Kb=λ+2μ/3,bl=n0Srηl/(krlκ),;下標s,b,l,g 分別代表土顆粒、土骨架、孔隙水以及氣體物理量的分量;上標“′”、“′′”分別表示一階、二階導;u表示位移矢量;p表示流體壓力;ρ 表示密度;λ,μ為土骨架Lame 常數(shù);γ 為有效應力系數(shù).K為壓縮模量;n0,Sr分別為土體孔隙率與飽和度;η 為流體黏滯系數(shù);krl,krg分別為液體和氣體滲透系數(shù);κ 為土的滲透率.

非飽和土中流體和氣體的滲流連續(xù)性方程為

式中,a11,a12,a13,a21,a22,a23表達式見附錄B.

通過弓入非飽和多孔介質的平均密度

孔隙水與固體骨架的相對位移

氣體與固體骨架間的相對位移

將上述波動方程方程(2a)寫成ub-v-w的表達形式有

滲流連續(xù)性方程(2b)寫成ub-v-w的表達形式有

土體總應力表達式

式中,σij為土單元總應力;e為體積應變;εij為應變分量;δij為克羅內克符號;p為等效孔隙流體壓力.其余各參數(shù)具體表達式見附錄C.

根據Helmholtz 矢量分解定理,式(4)中土骨架位移ub和流體、氣體相對于骨架的位移v,w,用標量勢和矢量勢分別表示為

式中,下標SH,SV,P 分別代表SH,SV,P 波的分量φ,φ,χ 分別為土骨架部分、液體部分以及氣體部分的勢函數(shù).ez表示z方向單位向量.

考慮穩(wěn)態(tài)響應,將式(5)代入式(4a),對時間t進行Fourier 變換至頻域后,得到如下兩組方程

式中ω 為角頻率.

為滿足微分方程組(6)有非零解,則需滿足系數(shù)矩陣行列式為0 的條件,得

式中,Ba,Bb,Bc,Bd,Be表達式見附錄D.

根據式(7),可得如下兩組Helmholtz 方程

式中,kp1,kp2,kp3分別為非飽和地基中的縱波波數(shù);ks為非飽和地基中的橫波波數(shù).

利用式(6)、式(8),經推導整理,各勢函數(shù)可表達為

式中,μ1l,μ2l,μ3l,μ1g,μ2g,μ3g,μsl,μsg的表達式見附錄E.

由于隧道埋置的非飽和半空間中存在兩個散射面,即圓柱形隧道外表面以及地表水平面,因此總的波場中包括外行的柱面波以及下行的平面波兩種波.參考文獻[16]并結合式(9),可得到直角坐標系下下行面波的位移勢函數(shù)為

式中,“—”代表y向波數(shù)域;分別為SH,SV,P1,P2,P3 波直角坐標下頻域波數(shù)域中下行波的位移勢函數(shù);ky,kz分別為y向與z向波數(shù);為S 波x向波數(shù),且虛數(shù)部非負;為P 波x向波數(shù),且虛數(shù)部非負.直角坐標下上行波位移勢函數(shù)只需將式(10a)中ksx,kp1,2,3x替換為-ksx,-kp1,2,3x即可.圓柱坐標系下外行波的位移勢函數(shù)為

將式(10)代入式(4b)得直角坐標系與圓柱坐標系下孔壓的表達式,如式(11)所示

將式(10)、式(11)代入式(4c)得直角坐標系下SH,SV,P1,P2,P3 頻域波數(shù)域中下行波的總應力勢函數(shù)表達式,以及圓柱坐標系下SH,SV,P1,P2,P3 頻域波數(shù)域中外行波的總應力勢函數(shù)表達式,分別為

式中,E1,E2,E3,E4,E5,E6 表達式見附錄G.

根據圓柱坐標系與直角坐標系下各類波的位移,孔壓以及應力函數(shù)表達式,得位移總場,孔壓總場以及應力總場的頻域— 波數(shù)域表達式,如式(13)所示

為滿足柱面波與平面波在同一邊界條件的耦合求解,弓入Bessel 函數(shù)與指數(shù)函數(shù)間的轉化關系[8],如式(14)所示

式中β=arcsin(ky/kr).

將式(14)代入式(10)～式(12)得柱面波函數(shù)與平面波函數(shù)間的轉化關系,如式(15)所示

式中,Imj(ky,ω),hj表達式見附錄H.

假定地表應力為0 且地表為透水邊界,結合式(13)和式(15),可得

由式(16)得

式中,Kjj′見附錄I.

將式(17)代入式(13)得圓柱坐標下總場表達式如下

1.3 耦合求解

根據式(1)殼體的平衡方程,利用隧道與地基土交界面處的位移與應力的連續(xù)性條件,可得

假設隧道襯砌不透水,不透氣,則隧,土交界面土體孔壓和氣壓的法向導數(shù)為0,即

2 算法的驗證

在本文提出的非飽和地基土-隧道動力響應計算方法中,若將非飽和地基土參數(shù)Sr,Se趨近于1,As趨近于0,則地基土可退化為兩相飽和地基土;若將非飽和地基土參數(shù)Sr趨近于Sw0,孔隙率n0趨近于0,則地基土退化為單相彈性地基土[30].為驗證本文算法的可靠性,將本文方法的計算結果分別與既有的2.5 維有限元-邊界元法以及全空間PiP 半解析法的計算結果進行對比分析.

2.1 單相彈性地基土

令Sr=0.06,n0=0.001,將本文計算方法中的地基土體退化為單相彈性地基土,與既有文獻[31-32]中的2.5 維有限元-邊界元法進行對比,計算參數(shù)參考文獻[31].固定單位簡諧載荷作用于隧道仰拱(x=-2.75 m,y=0 m,z=0 m)處,分別取地表兩個觀察點(O1:x=10 m,y=0 m,z=0 m;O2:x=10 m,y=20 m,z=0 m)進行對比驗證,土體豎向位移響應隨頻率的變化曲線如圖2 所示.可以看到,本文算法的計算結果與既有單相彈性地基土的2.5 維有限元-邊界元法的計算結果吻合較好,驗證了將非飽和地基土退化為單相彈性地基土時,本文方法計算結果的可靠性.

2.2 兩相飽和地基土

圖2 退化后算法的驗證(非飽和地基土退化為單相彈性地基土)Fig.2 Verification of the method after the unsaturated soil degenerate into single-phase elastic soil

令Sr=0.999,As=0,將本文的非飽和地基土退化為兩相飽和地基土,與既有的飽和地基土2.5 維有限元-邊界元法進行對比,計算參數(shù)參考文獻[20].計算時移動單位簡諧載荷(v0=16.7 m/s,f0=20 Hz)作用于隧道仰拱(x=-2.75 m,y=0 m,z0=z-v0t=0 m)處.圖3 給出了隧道底部土體位移與孔壓幅值沿隧道軸向的分布.從圖3 可以看到,本文方法的計算結果與既有的飽和土2.5 維有限元-邊界元法的計算結果吻合較好,驗證了將本文計算方法中非飽和地基土退化為兩相飽和地基土時,本文方法計算結果的可靠性.

圖3 退化后算法的驗證(非飽和地基土退化為飽和地基土)Fig.3 Verification of the method after the unsaturated soil degenerate into two-phase saturated soil

2.3 全空間非飽和地基土

進一步驗證本文提出的非飽和半空間地基土-隧道系統(tǒng)動力響應計算方法的可靠性.由于既有文獻中未見關于非飽和半空間地基土-隧道系統(tǒng)動力響應的研究,為此,選取既有的非飽和全空間地基土的PiP 半解析法[29]進行對比驗證.分別選取隧道中心線距離地表豎直距離d=25 m 與d=5 m 時進行計算和對比分析,隧道半徑r1=3 m,隧道襯砌厚度h=0.25 m.隧道襯砌為混凝土材料,其楊氏模量Et=50 GPa,泊松比υt=0.3,密度ρt=2500 kg/m3,材料阻尼為βt=0.03.非飽和土的參數(shù)參考文獻[23,33].計算移動單位簡諧載荷(v0=16.7 m/s,f0=20 Hz)作用于隧道仰拱(x=-2.75 m,y=0 m,z0=0 m)處時,隧道底部土體的動力響應,結果如圖4 所示.由圖4 可以發(fā)現(xiàn),本文方法計算的隧道底部土體的動力響應與全空間PiP 半解析法的計算結果基本吻合,驗證了本文計算方法的可靠性.

圖4 非飽和地基土對比驗證Fig.4 Verification by comparison with the results obtained by existing tunnel model

圖4 非飽和地基土對比驗證(續(xù))Fig.4 Verification by comparison with the results obtained by existing tunnel model(continued)

3 算例分析

基于本文提出的算法,重點研究固定簡諧載荷作用下飽和度對地基土-隧道系統(tǒng)動力響應的影響.圓形隧道襯砌和非飽和地基土計算參數(shù)如表1 所示,其中,隧道襯砌計算參數(shù)參考文獻[33]中上海地鐵參數(shù),非飽和地基土參考文獻[23]給出的砂土參數(shù),如表1 所示.

圖5 給出了固定簡諧載荷作用于隧道仰拱(x=-2.75 m,y=0 m,z0=z-v0t=0 m)處時,不同位置處(O1,O2,O3,O4)土體豎向位移幅值隨激振頻率f0的變化曲線,計算時考慮了不同的飽和度(Sr=1,0.9,0.7,0.5).由圖5 可以看到,在相同的激振頻率但不同的飽和度下,土體豎向位移幅值存在差異,這是由于飽和度的改變一方面會弓起土體有效應力的改變,另一方面會弓起土體動剪切模量的變化,兩者共同作用使得不同飽和度狀態(tài)下土體位移幅值存在差異.

表1 隧道襯砌和非飽和土計算參數(shù)Table 1 Calculation parameters of tunnel lining and unsaturated foundation soil

圖5 不同地表觀測點位移頻響曲線Fig.5 Frequency response curves of soil displacement of different observation points on ground surface

圖5 不同地表觀測點位移頻響曲線(續(xù))Fig.5 Frequency response curves of soil displacement of different observation points on ground surface(continued)

從圖5 中還可以看到,土體豎向位移隨載荷激振頻率的變化曲線呈現(xiàn)出明顯的震蕩現(xiàn)象,這主要是由于地表瑞雷面波和橫,縱波的干涉效應所致[6].同時在不同飽和度下,豎向位移頻響曲線極值點出現(xiàn)位置的不同.根據Forrest 等[6]的研究,頻響曲線極值點出現(xiàn)的位置與土體中波的傳播速度存在相關性.而從非飽和土控制方程求解中可得,彈性波波速表達式c=ω/Re(k)中k值與飽和度存在關聯(lián).因此,土體飽和度的變化會弓起頻響函數(shù)極值點位置的改變.

圖6 為不同飽和度下隧道下方觀測點O3和O4處土體超孔隙水壓力的頻響曲線.由圖6 可以看到,飽和度越小,土體超孔隙水壓力越小.此外,還可以看到,當飽和度從1.0 變?yōu)?.9 時,超孔隙水壓力下降幅度大,其原因在于當土體飽和度接近1 時,氣體以閉氣泡形式存在于液體中,由于氣體的體積模量遠小于液體的體積模量,因此隨著氣體的含量的增大,等效孔隙流體的體積模量迅速下降,使得等效孔隙流體分擔的力減小,造成超孔隙水壓力迅速下降.綜合圖5 和圖6 可知,飽和度對系統(tǒng)動力響應存在較大影響,對于非飽和地基土,計算隧道動力響應時應考慮飽和度的影響.

圖6 不同觀測點超孔隙水壓力頻響曲線Fig.6 Frequency response curves of pore water pressure at different observation points

4 結論

(1)本文提出了非飽和土-隧道系統(tǒng)動力響應計算的波函數(shù)法.該方法的計算結果與既有的2.5 維有限元-邊界元法以及全空間PiP 半解析法的計算結果均吻合較好,驗證了本文計算方法的可靠性.

(2)算例分析結果表明,飽和度對土體位移與超孔隙水壓力的幅值有較大影響.因此,對于非飽和地層中的隧道系統(tǒng),計算其系統(tǒng)動力響應時,將土體視為三相介質而非單相介質或兩相介質是有必要的.

(3)本文算法中的三相介質(非飽和地基土)可退化為兩相介質(飽和地基土)或單相介質(彈性地基土).通過參數(shù)退化,結合位移應力協(xié)調以及滲流連續(xù)等條件,該方法可進一步發(fā)展為單相彈性地基土,飽和地基土,非飽和地基土并存條件下分層地基土-隧道系統(tǒng)動力計算方法.

附錄

附錄A

其中,ξ 為殼z方向波數(shù);h為殼的厚度;E為殼體楊氏模量;ν 為殼體泊松比;ρt為殼體密度;qr,qθ,qz分別為襯砌殼中心面沿r,θ,z向的凈應力;f0為激振頻率.

附錄B

其中,Sw0為液體的約束飽和度;α1,α2,α3分別為擬合參數(shù),υs為土體泊松比;φ 土飽和時的內摩擦角;μs分別為土體飽和狀態(tài)下的動剪切模量.

附錄C

附錄D

附錄E

附錄F

附錄G

附錄H

附錄I