飽和多孔介質的二維Green函數及地震反應的BEM計算

魏 綱， 宋宥整， 丁伯陽

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

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

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

1 動力學方程的解

飽和多孔介質的Biot動力學控制方程可寫為[11-12]：

(1)

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

(2)

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

(3)

式中：λc，μ，α和M可以看作是4個獨立的材料常數。

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

(4)

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

(5)

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

(6)

2 利用De Hoop法求二維Green函數

2.1 二維固相Green函數gij

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

(7)

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

2.2 二維流相Green函數,g3i,gi3和g33

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

(8)

(9)

(10)

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

(11)

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

(12)

(13)

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

圖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 二維應力Green函數σikj,σik3

根據本構關系，應力Green函數可寫為式(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), 得到二維應力Green函數如式(17)與式(18)所示。

(17)

(18)

將表1的參數分別代入式(17)和式(18)中，σikj,σik2的部分結果如圖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函數BEM的數值實現

3.1 時域邊界積分方程的離散形式

Somigliana表象用于BEM中Green函數的數值積分實現，關于兩相飽和介質的Somigliana積分方程已經由Chen[22]在1994年推導出，可以寫成如下形式：

(19)

(20)

(21)

(22)

為了實現式(21)和(22)的邊界積分數值計算，邊界上的變量(如位移，應力和孔隙壓力)需要進行離散。總時間t應該在時域上分成N段等時間步長Δt。對于平面域，邊界Γ也應分成包括結點的q個單元。因此需要由單元中的相關結點處的位移，應力和孔隙壓力的值，通過插值函數獲得在時間τ時，邊界中任意點處的位移，應力和孔隙壓力值。這些插值公式分別是：

(23)

(24)

(25)

(26)

(27)

(28)

(29)

(30)

(31)

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

(32)

(33)

(34)

式中：Γj表示積分區域。

3.2 BEM積分的奇異性處理

時間和區域積分可以用Gauss積分求解，這在有關計算平臺(如，MathLab等)幾乎不會有困難。

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

3.3 工程實例的數值計算

(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平臺中求解得到一系列代數聯立方程組。用于積分計算的6×6標準Gauss積分可以滿足精度要求。相關位移結果分別如圖12和13所示。

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

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

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

(36)

式中：A(t)是相關的地震記錄，gzz(X/V,t)是圖12和13所表達的函數。分別選取EI Centro和汶川2008地震兩條形態相差較大的典型記錄，如圖14和15，代入式(36)積分，結果見圖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地震表面的位移響應Fig.16 Displacement response on surface for EI Centro quake

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

由計算結果可以發現隧道的地震放大作用明顯；對于EI Centro或汶川2008地震記錄，隧道底部的響應均要大于中間和頂部的響應。另外，由于在本計算中土參數取為線性，響應結果的頻率與輸入記錄的原頻率相差不大。

4 結 論

本文通過對飽和多孔介質三維U-P形式Green函數進行積分，得到了相應的二維Green函數。在應用Somigliana表象的BEM數值計算中，得到了可信的結果，從而得出以下結論：

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

(2)本文提出的U-P形式二維Green函數形式簡潔，可完成兩相飽和介質動態響應的計算，也有助于多孔飽和介質響應的計算。

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

附錄A

由坐標變換z=rtanθ, 有：

(A1)

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

(A2)

結合δ函數的性質, 可以推導出：

(A3)

那么

(A4)

同樣，可以得到

(A5)

(A6)

(A7)