饱和多孔介质的二维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)