非飽和土-隧道系統(tǒng)動(dòng)力響應(yīng)計(jì)算的波函數(shù)法1)

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

(同濟(jì)大學(xué)道路與交通工程教育部重點(diǎn)實(shí)驗(yàn)室,上海 201804)

(同濟(jì)大學(xué)上海市軌道交通結(jié)構(gòu)耐久與系統(tǒng)安全重點(diǎn)實(shí)驗(yàn)室,上海 201804)

引言

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

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

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

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

1 控制方程和求解

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

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

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

1.1 襯砌振動(dòng)控制方程

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

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

1.2 非飽和地基土波動(dòng)方程

將非飽和地基土視為流、固、氣組成的三相介質(zhì),其實(shí)用波動(dòng)方程如下[23]

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

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

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

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

孔隙水與固體骨架的相對(duì)位移

氣體與固體骨架間的相對(duì)位移

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

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

土體總應(yīng)力表達(dá)式

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

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

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

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

式中ω 為角頻率.

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

式中,Ba,Bb,Bc,Bd,Be表達(dá)式見附錄D.

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

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

利用式(6)、式(8),經(jīng)推導(dǎo)整理,各勢(shì)函數(shù)可表達(dá)為

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

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

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

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

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

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

根據(jù)圓柱坐標(biāo)系與直角坐標(biāo)系下各類波的位移,孔壓以及應(yīng)力函數(shù)表達(dá)式,得位移總場(chǎng),孔壓總場(chǎng)以及應(yīng)力總場(chǎng)的頻域— 波數(shù)域表達(dá)式,如式(13)所示

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

式中β=arcsin(ky/kr).

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

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

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

由式(16)得

式中,Kjj′見附錄I.

將式(17)代入式(13)得圓柱坐標(biāo)下總場(chǎng)表達(dá)式如下

1.3 耦合求解

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

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

2 算法的驗(yàn)證

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

2.1 單相彈性地基土

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

2.2 兩相飽和地基土

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

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

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

2.3 全空間非飽和地基土

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

圖4 非飽和地基土對(duì)比驗(yàn)證Fig.4 Verification by comparison with the results obtained by existing tunnel model

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

3 算例分析

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

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

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

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

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

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

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

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

4 結(jié)論

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

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

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

附錄

附錄A

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

附錄B

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

附錄C

附錄D

附錄E

附錄F

附錄G

附錄H

附錄I