魏 綱， 宋宥整， 丁伯陽

(1.浙江大學(xué) 建筑工程學(xué)院，杭州 310058； 2.浙江大學(xué)城市學(xué)院 土木工程系，杭州 310015)

巖土工程中的很多問題，都可以簡化為平面應(yīng)變問題討論[1-5]，因此飽和多孔介質(zhì)二維Green函數(shù)在土動(dòng)力學(xué)中有廣泛應(yīng)用。目前國內(nèi)外關(guān)于二維Green函數(shù)求解的方法主要有：首先，De Hoop等[6]提出二維Green函數(shù)可以由對(duì)應(yīng)的三維Green函數(shù)沿z軸在(-∞,+∞)上積分得到，Manolis等[7]證明了此法精確可行。其次，Chen[8]利用無量綱參數(shù)研究了U-P形式下飽和多孔介質(zhì)二維Green函數(shù)。由于丁伯陽等[9-10]在早期文獻(xiàn)中已經(jīng)得到了飽和多孔介質(zhì)的三維位移Green函數(shù)，所以本文依然采納De Hoop的方法得到二維Green函數(shù)，并首次推導(dǎo)出了它的流相Green函數(shù)。

事實(shí)上，飽和多孔介質(zhì)Biot動(dòng)力方程只需要4個(gè)未知量的U-P(位移和孔壓)形式解答，6個(gè)未知量的U-W形式(W流相相對(duì)位移)解答存在多余的增根。另外，考慮到動(dòng)力響應(yīng)評(píng)價(jià)的需要，受集中力δ作用的二維Green函數(shù)U-P解答也在本文中導(dǎo)出；利用它和Somigliana積分方程可實(shí)現(xiàn)BEM數(shù)值方法對(duì)隧道地震響應(yīng)的計(jì)算。從而在理論或工程領(lǐng)域便于飽和多孔介質(zhì)動(dòng)力響應(yīng)的研究。

1 動(dòng)力學(xué)方程的解

飽和多孔介質(zhì)的Biot動(dòng)力學(xué)控制方程可寫為[11-12]：

(1)

γ(ω)=α0(ω)ρf/β0+ηi/[ωkd(ω)]

(2)

如果ks，kf和kb分別為固相，流相和兩相介質(zhì)的體積模量，那么有[19-20]：

(3)

式中：λc，μ，α和M可以看作是4個(gè)獨(dú)立的材料常數(shù)。

動(dòng)力方程(1)，在三維狀態(tài)受到脈沖力δ(t)的Green函數(shù)在頻域中的解答可寫為[21]：

(4)

式中：R是場點(diǎn)x和源點(diǎn)V之間的距離，R=[(xi-Vi)(xi-Vi)]1/2。xi,Vi分別是x和V在i方向的坐標(biāo)。Kα1,Kα2和Kβ分別是快，慢縱波和橫波的波數(shù)。α1,α2和β分別是在飽和多孔介質(zhì)中快，慢縱波和橫波的速率。δij是Kronecker-δ。定義(丁伯陽，1999)：

(5)

式中:n=1, 2。λ1和λ2被定義為縱波解耦系數(shù)(丁伯陽，1999；Ding, 2013, 2016)。式(4)表示源點(diǎn)V在j方向受到單位集中力作用,在場點(diǎn)x的i向位移。用xi替換xi-ζi通過頻-時(shí)域的Stokes公式互換，可得到時(shí)域三維Green函數(shù)如式(6)，基于Ding等的工作，式(6)的表達(dá)已得到了改進(jìn)[23]。

(6)

2 利用De Hoop法求二維Green函數(shù)

2.1 二維固相Green函數(shù)gij

De Hoop提出，在(x1,x2)平面內(nèi)的二維位移場Green函數(shù)gij，可以從其對(duì)應(yīng)的三維解答Gij沿z軸在無窮域積分得到。因此，時(shí)域中受脈沖力δ(t)的二維Green函數(shù)可以通過對(duì)式(6)沿z軸在(-∞,+∞)積分得到：(詳見附錄A)。

(7)

式中：r=[(xi-Vi)(xi-Vi)]1/2(i,j=1, 2)是二維平面上場點(diǎn)x和源點(diǎn)V之間的距離，并且R=[r2+(x3-V3)2]1/2。當(dāng)ρf=0,λ1=0,λ2=1時(shí)，式(7)中飽和多孔介質(zhì)Green函數(shù)容易退化到單相介質(zhì)Green函數(shù)解答。

2.2 二維流相Green函數(shù),g3i,gi3和g33

三維流相Green函數(shù)可以分別寫成如式(8)～(10)的形式：

(8)

(9)

(10)

式中:G4i被定義為三維中由于源點(diǎn)V處i方向作用單位δ(t)力，在場點(diǎn)x處的孔壓[22]，如果g3i表示二維中的上述定義，那么式(11)能夠在z軸(-∞,+∞)區(qū)間對(duì)式(8)積分得到。積分過程可詳見附錄A。

(11)

同理，二維Green函數(shù)gi3被定義為在源點(diǎn)處V將單位δ(t)流體注入到孔隙中，在場點(diǎn)x處i方向上的位移。對(duì)式(9)關(guān)于z軸在(-∞,+∞)進(jìn)行積分,可以得到對(duì)應(yīng)的二維Green函數(shù)gi3如下：

(12)

(13)

采用Chen在表1中的無量綱參數(shù)，代入式(11)～(13)，可得流相二維Green函數(shù)脈沖曲線見圖1～5。另外，Chen用Laplace變換獲得過Heaviside力作用的二維流相Green函數(shù)。由于Heaviside力作用的Green函數(shù)可由δ(t)力作用結(jié)果的時(shí)間積分得到，也將表1無量綱參數(shù)分別代入式(7)及式(11)～(13)的時(shí)間積分結(jié)果中，并與Chen的結(jié)果進(jìn)行對(duì)比，發(fā)現(xiàn)兩者吻合且穩(wěn)定，表明本文的二維Green函數(shù)可靠。

圖1 式(11)中的脈沖曲線g31Fig.1 Impulse curve g31of in Eq.(11)

圖2 式(11)中的脈沖曲線g32Fig.2 Impulse curve g32of in Eq.(11)

圖3 式(12)中的脈沖曲線g13Fig.3 Impulse curve g13of in Eq.(12)

圖4 式(12)中的脈沖曲線g23Fig.4 Impulse curve of g23in Eq.(12)

圖5 式(13)中的脈沖曲線g33Fig.5 Impulse curve g33of in Eq.(13)

2.3 二維應(yīng)力Green函數(shù)σikj,σik3

根據(jù)本構(gòu)關(guān)系，應(yīng)力Green函數(shù)可寫為式(14)和(15)所示的公式：

σikj=λcδikgmj,m+μ(gij,k+gkj,i)-αδikg3j

(14)

σik3=λδikgm3,m+μ(gi3,k+gk3,i)-αδikg33

(15)

式(14)和(15)中：

(16)

將式(7), (11), (12)和(13)代入式(14)和(15), 得到二維應(yīng)力Green函數(shù)如式(17)與式(18)所示。

(17)

(18)

將表1的參數(shù)分別代入式(17)和式(18)中，σikj,σik2的部分結(jié)果如圖6～8所示。

圖6 式(17)中的脈沖曲線σ111Fig.6 Impulse curve σ111of in Eq.(17)

圖7 式(17)中的脈沖曲線σ121Fig.7 Impulse curve σ121of in Eq.(17)

圖8 式(18)中的脈沖曲線σ123Fig.8 Impulse curve σ123of in Eq.(18)

3 二維Green函數(shù)BEM的數(shù)值實(shí)現(xiàn)

3.1 時(shí)域邊界積分方程的離散形式

Somigliana表象用于BEM中Green函數(shù)的數(shù)值積分實(shí)現(xiàn)，關(guān)于兩相飽和介質(zhì)的Somigliana積分方程已經(jīng)由Chen[22]在1994年推導(dǎo)出，可以寫成如下形式：

(19)

(20)

(21)

(22)

為了實(shí)現(xiàn)式(21)和(22)的邊界積分?jǐn)?shù)值計(jì)算，邊界上的變量(如位移，應(yīng)力和孔隙壓力)需要進(jìn)行離散?？倳r(shí)間t應(yīng)該在時(shí)域上分成N段等時(shí)間步長Δt。對(duì)于平面域，邊界Γ也應(yīng)分成包括結(jié)點(diǎn)的q個(gè)單元。因此需要由單元中的相關(guān)結(jié)點(diǎn)處的位移，應(yīng)力和孔隙壓力的值，通過插值函數(shù)獲得在時(shí)間τ時(shí)，邊界中任意點(diǎn)處的位移，應(yīng)力和孔隙壓力值。這些插值公式分別是：

(23)

(24)

(25)

(26)

(27)

(28)

(29)

(30)

(31)

式中：ξ3=-ρf/γ。對(duì)于區(qū)域積分，有式(32)～(34)：

(32)

(33)

(34)

式中：Γj表示積分區(qū)域。

3.2 BEM積分的奇異性處理

時(shí)間和區(qū)域積分可以用Gauss積分求解，這在有關(guān)計(jì)算平臺(tái)(如，MathLab等)幾乎不會(huì)有困難。

至于奇點(diǎn)的處理，BEM中對(duì)于單相介質(zhì)已有成熟的方法。對(duì)于二維情況，lnr型的奇異性可直接通過高斯積分得到；1/r型的奇異性可以通過雅可比行列式的坐標(biāo)變換來去除；1/r2型的奇異性可以通過柯西主值積分來解決。由于這里快，慢縱波已經(jīng)解耦，所以兩相問題完全可以根據(jù)單相的方法進(jìn)行處理，此處不再贅述。

3.3 工程實(shí)例的數(shù)值計(jì)算

(35)

圖9 飽和土半空間中圓形斷面隧道Fig.9 Tunnel of a circle section inthe half space of the saturated soil

圖10 隧道頂部影響曲線Fig.10 Influence curve in the top of tunnel

圖11 隧道頂部影響曲線Fig.11 Influence curve in the top of tunnel

通過離散邊界積分方程(26)和(27)，在MatLab平臺(tái)中求解得到一系列代數(shù)聯(lián)立方程組。用于積分計(jì)算的6×6標(biāo)準(zhǔn)Gauss積分可以滿足精度要求。相關(guān)位移結(jié)果分別如圖12和13所示。

圖12 表面位移隨時(shí)間的變化曲線Fig.12 Curves of surface displacement vs time

圖13 隧道頂部，中間和底部位移隨時(shí)間的變化曲線Fig.13 Curves of displacement vs time in the bottom,the middle and top of the tunnel

由于圖12和13是δ力作用下的位移反應(yīng)，根據(jù)Duhamel積分公式，得到：

(36)

式中：A(t)是相關(guān)的地震記錄，gzz(X/V,t)是圖12和13所表達(dá)的函數(shù)。分別選取EI Centro和汶川2008地震兩條形態(tài)相差較大的典型記錄，如圖14和15，代入式(36)積分，結(jié)果見圖16～19。

圖14 EI Centro的地震加速度記錄Fig.14 Acceleration record of EI Centro quake

圖15 汶川2008年地震加速度記錄Fig.15 Acceleration record of Wenchuan 2008 quake

圖16 EI Centro地震表面的位移響應(yīng)Fig.16 Displacement response on surface for EI Centro quake

圖19 對(duì)于汶川2008地震的隧道頂部，中部和底部的位移響應(yīng)Fig.19 Displacement response in the top, middleand bottom of the tunnel for Wenchuan 2008 quake

由計(jì)算結(jié)果可以發(fā)現(xiàn)隧道的地震放大作用明顯；對(duì)于EI Centro或汶川2008地震記錄，隧道底部的響應(yīng)均要大于中間和頂部的響應(yīng)。另外，由于在本計(jì)算中土參數(shù)取為線性，響應(yīng)結(jié)果的頻率與輸入記錄的原頻率相差不大。

4 結(jié) 論

本文通過對(duì)飽和多孔介質(zhì)三維U-P形式Green函數(shù)進(jìn)行積分，得到了相應(yīng)的二維Green函數(shù)。在應(yīng)用Somigliana表象的BEM數(shù)值計(jì)算中，得到了可信的結(jié)果，從而得出以下結(jié)論：

(1)本文提出的U-P形式二維Green函數(shù)可作為Somigliana表象邊界元積分的核函數(shù)。

(2)本文提出的U-P形式二維Green函數(shù)形式簡潔，可完成兩相飽和介質(zhì)動(dòng)態(tài)響應(yīng)的計(jì)算，也有助于多孔飽和介質(zhì)響應(yīng)的計(jì)算。

(3)本文利用Green函數(shù)求得相關(guān)隧道單位脈沖力δ作用下動(dòng)力響應(yīng)的BEM計(jì)算，通過Duhamel公式積分求解了EI Centro和2008年汶川地震加速度記錄下的飽和土隧道的地震響應(yīng)計(jì)算，應(yīng)用此種方法可以較為方便地進(jìn)行隧道的地震動(dòng)位移反應(yīng)計(jì)算。

附錄A

由坐標(biāo)變換z=rtanθ, 有：

(A1)

假設(shè)t-rsecθ/αi=x, 式(A1)成為

(A2)

結(jié)合δ函數(shù)的性質(zhì), 可以推導(dǎo)出：

(A3)

那么

(A4)

同樣，可以得到

(A5)

(A6)

(A7)