2. 北京军区 66469部队, 北京 100042
2. Militaryarea of Beijing, Unit 66469, Beijing 100042, China
地震波在从震源通过各种途径传播时,通常要穿过饱和土层和非饱和土层最后到达地面,地壳上部的连续固体骨架中包含着相互连通的孔隙和槽路,当将充满液体的土层进行动力分析时,按饱和多孔介质理论比按单相介质理论分析更为合理.与单相介质相比,饱和多孔介质中的波动问题要复杂得多,近年来饱和多孔介质人工边界及其数值模拟问题受到越来越多学者的重视.Modaressi等[1]首先将人工边界上的u-p表达式引入饱和多孔介质的有限元分析中;Degrande等[2]通过引入势函数、波场分解以及Fourier变换给出人工边界上频域内的应力和位移的表达式,实际上是DtN边界在法向上的表示;Akiyoshi等[3]通过波场分解和Fourier变换得到频域内固、液相位移、总应力和流体压力的表达式,利用旁轴近似和Fourier反变换得到人工边界上应力与u-p、u-U、u-W的黏性边界表达形式;Gajo[4]建立两种极端情形下(即黏性耦合很小和很大时的情形)饱和介质一维单侧波动方程,得到类似单相介质的黏性边界,然后根据Higdon的方法,得到高阶多向人工边界条件;Akiyoshi等[5]采用等效拉梅常数和平均渗透系数,通过线性变换,将各向异性波动方程等效为各向同性,又将此边界推广为考虑各向异性的人工边界;邵秀民等[6]基于透射边界,提出一种饱和介质的人工边界条件,并将其应用到二维有限元计算分析中;Zerfa[7]根据Bowen的公式[8]提出可考虑排水特性的黏性人工边界;刘光磊、宋二祥[9]根据Biot动力固结理论外行柱面波假设,提出了一种适用于饱和多孔介质时域动力分析的二维黏弹性边界条件;王子辉、赵成刚等[10]从二维平面外行柱面波出发,提出了一种饱和多孔介质黏弹性动力人工边界.杜修力、李立云[11]采用平面波和远场散射波经验叠加来反映外行波传播,并在推导过程中考虑了多角度透射的影响,提出了饱和多孔介质近场波动分析的一种黏弹性人工边界处理方法.
上述研究中,多数是针对二维问题展开的,有关饱和多孔介质的三维动力人工边界及其数值模拟的研究并不多见.本文基于饱和多孔介质波动理论,采用与文献[10]类似的思路提出饱和多孔介质三维黏弹性动力人工边界处理方法,同时根据本文作者已经建立的饱和多孔介质二维显式有限元数值计算方法[12],提出该理论的三维数值模拟方法,并开发了实现该三维方法的有限元程序.最后将建立的人工边界条件和饱和多孔介质三维动力反应分析的显式有限元法相结合,用于饱和多孔介质的动力研究.
2 饱和多孔介质三维显式有限元方程的建立Biot[13]提出了饱和多孔介质的波动方程
(1) |
式中λ,μ,Q,R为弹性常数;ρ11,ρ22,ρ12为动力质量系数,ρ11=(1-n)ρs+ρa,ρ22=nρf+ρa,ρ12=-ρa;ρs和ρf分别为固体骨架和孔隙流体的质量密度,ρa为描述流固两相耦合的附加质量密度;e和ε分别表示固相骨架和孔隙流体的体积应变,e=div u,ε=div U;b为与渗流有关的系数,b=ζn2/k,ζ为流体滞变系数,n为孔隙率,k为渗透系数.
2.1 波动方程的矩阵表达形式由式(1)可以得到该波动方程的矩阵表达式
(2) |
(3) |
其中
设N为单元形函数矩阵,则
(4) |
式中ue和Ue分别为单元e的节点固相和液相位移;ue和Ue分别为单元e的固相和液相位移.
根据文献[12],方程(2)、(3)的局部节点系GALERKIN弱式为
(5) |
其中应变矩阵B=L·N,L为三维问题的微分算子;fu e和fU e分别为作用在单元e边界上的固相和液相外力.
2.3 离散方程及解耦方法对于三维问题,采用8节点六面体等参单元,节点i和与其直接相邻的节点k所构成局部节点系,此节点系中包含27个节点.根据式(5)可以直接得到如下离散公式
(6) |
式中,L代表节点i周围的不同单元,因只有节点i周围8个单元对其有影响,所以只给出8个单元的求和,j代表同一单元内的8个节点.波动问题数值模拟的精度要求有限元的尺寸同波长相比应足够小,故在每一单元内惯性力的空间变化很小,可在一个单元内的加速度为常量,即等于其节点的加速度,
(7a) |
其中,
(7b) |
式(6),(7)中,
利用有限元方法建立离散的动力方程后,采用中心差分法与Newmark平均加速度近似格式相结合的方法求解节点动力方程式(6).
当t=pΔt时刻,将式(7)代入式(6)得到pΔt时刻i节点的运动方程:
(8) |
在(p+1)Δt时刻,将式(7)代入式(6),并与式(8)相加得到
(9) |
(1)节点位移动力反应的显式
根据时间的中心差分法
(10) |
有
(11) |
将式(11)代入式(8)可以得到用前一时刻的位移和速度表示的当前时刻位移的显式表达式:
(12a) |
(12b) |
(2)节点速度的动力反应显式
根据Newmark常平均加速度法,有
(13) |
将式(13)代入式(12),便可得到用前一时刻的位移、速度以及当前时刻的位移表示的当前时刻速度的显式表达式:
(14a) |
(14b) |
(3)节点加速度动力反应的显式
由式(13)可得用前一时刻的速度和此时刻的速度表示的加速度的显式表达式:
(15) |
式(13),(14),(15)为考虑耦合质量和流体黏性的饱和多孔介质局部节点系下某一节点的动力反应的显式表达式,该节点的动力反应只与相邻节点的动力反应有关.若不考虑耦合质量,即令ρa=0,则有
对线弹性饱和多孔介质,本构关系为[13]
(16) |
其中,σij为固体骨架应力张量,s为孔隙流体应力,e,ε分别为固体骨架和孔隙流体的体应变,eij=(ui,j+uj,i)/2,δij为Kronecker符号,λ,μ,Q,R为弹性常数.
3.2 人工边界应力条件的推导 3.2.1 人工边界切向应力因假设孔隙流体无黏性,所以,孔隙流体不承担剪力,在切向只有固体骨架承担剪应力.考虑扩散的球面剪切波,球坐标系球面剪切波动位移的近似解为
(17) |
由式(17)所确定的剪应变和剪应力分别为:
(18) |
(19) |
在坐标r处质点的运动速度为
(20) |
由式(17)、(18)、(19)可得在波阵面上
(21) |
式(21)即为三维切向人工边界方程,该方程给出了波阵面上切向应力与位移的关系式.
对于如式(21)所示的三维切向边界方程,容易证明该方程等价于并联的弹簧-阻尼器系统,对应物理元件的参数为
(22) |
其中,k为弹簧刚度系数,c为黏性阻尼系数.
3.2.2 人工边界法向应力忽略饱和多孔介质中慢P波的影响,只考虑球面快P波[9],球坐标系下其波动方程为
(23) |
其中ϕ为位移势函数,cp为饱和多孔介质的快P波波速,r为径向坐标.方程(23)的通解可表示为
(24) |
式中f(·)和g(·)为任意函数,分别表示外行扩散波和内行会聚波.
考虑外行扩散波
(25) |
首先分析固体骨架的法向应力,固体骨架法向位移为
(26) |
法向与切向应变为
(27) |
(28) |
由式(16)固体骨架法向应力为
(29) |
其中,η1为液相参与系数.
固体骨架速度、加速度为
(30) |
(31) |
对式(26)求导可得
(32) |
(33) |
令G=λ+Qη1+4μ,ρcx2=λ+2μ+Qη1,将式(26)、(27)和(28)代入式(29),得
(34a) |
(34b) |
将式(34b)乘以系数
(34c) |
将式(34c)与式(34a)相加,并结合式(30),(31)得
(35) |
式(35)与图 1所示的弹簧-阻尼-质量模型的动力系统运动方程式(36)等价
(36) |
在人工边界r=rb处,等效弹簧、黏性阻尼和质量系数为
(37) |
将图 1中的质量忽略,将阻尼的另一端直接固定,成为弹簧-阻尼模型,如图 2所示.这样简化并不明显影响人工边界的精度,但给实际应用带来极大的方便[15].简化后的三维黏弹性人工边界如图 3所示,图中坐标x、y为沿人工边界的切向方向,z为为沿人工边界的法向方向.
对于孔隙流体,因忽略慢P波,
(38) |
由式(29)得
(39) |
又由式(26)~(29)得
(40) |
将式(39)、(40)代入式(38),并注意到u=U/η1,
(41) |
其中,U,
假设在半空间问题分析中采用直角坐标系下的矩形人工边界.球面人工边界上的法向矢量与切向矢量构成局部坐标系,固体骨架应力矢量用矩阵表达为
(42) |
其中
zl为位移矢量,k φ,cφ,kθ,cθ分别为固体骨架切向弹簧刚度系数和黏性阻尼系数,kr,cr分别为固体骨架径向弹簧刚度系数和黏性阻尼系数.
局部坐标系下的应力与位移矢量与直角坐标系下的应力与位移矢量的变换关系为
(43) |
(44) |
其中,
T为坐标转换矩阵.将式(43),(44)代入式(42),得直角坐标系下的应力矢量为
(45) |
其中
孔隙流体应力在直角坐标系下表达为
(46) |
式(45)、(46)即为直角坐标系下矩形平面人工边界上的应力表达式.σ和si分别为人工边界上某一点的固体骨架和孔隙流体的应力矢量,在有限元计算中,可简单地将边界上节点应力乘以相邻四个单元面积的平均值,作为节点集中力列入到节点动力平衡方程中,并集装到总体刚度矩阵和阻尼矩阵中.
4 数值计算算例1、2分别给出两种不同荷载形式下地面一点的位移曲线,算例中把应用本文提出的人工边界条件的数值计算结果与应用黏性边界[3]及一、二阶透射边界[14]、扩展人工边界的解进行了对比分析.其中,扩展有限元解是将人工边界远置,使得在计算时间段内从人工边界反射回来的波不能到达所考虑的点或局部区域.
算例1:三维半空间情况下,在自由面中央A点作用一水平集中动荷载,见图 4,P(t)=P0sinωt,P0=1N,地震工程感兴趣的频段为:0.1~10 Hz.选择频率分别为0.1 Hz,1 Hz,3 Hz,5 Hz.计算区域:12 m×12 m×6 m,采用立方体单元离散计算区域,假设固相不可压缩,材料参数见表 1.计算在水平集中力作用下,作用点处水平位移时程曲线,计算时间为0.5 s,见图 5.图中分别给出了本文边界、黏性边界、一、二阶透射边界对应的数值解,精确解采用扩展有限元解.采用本文提出的有限元程序进行计算.图 5表明,在水平动荷载作用下,尽管应用本文提出的人工边界条件的数值计算结果与扩展有限元解之间略有误差,但它的数值解与扩展有限元解是最接近的,而且不论动载是在低频(f=0.1 Hz)或高频(f=10.0 Hz),其数值计算结果很稳定;一阶透射边界解与扩展边界解相差较远,这是由于本模型人工边界尺寸较小,一阶透射边界结果受不同模型影响明显,人工边界尺寸越大精度越高,这给透射边界应用带来不便.
算例2:如图 4所示,自由面有一宽度2 m竖向作用的条形简谐变化荷载,P(t)=P0sinωt,P0=1N.仍选择频率分为:0.1 Hz,1 Hz,3 Hz,5 Hz,计算域的尺寸和材料常数同算例1.计算得到地表正中央A点在不同荷载频率和人工边界条件下的竖向位移时程曲线,计算时间为0.5 s,见图 6.由图可见,在竖向动荷载作用下,尽管本文黏弹性人工边界的数值计算结果与扩展有限元解之间略有误差,但它的数值解与扩展有限元解是接近的,低频(f=0.1 Hz)时比二阶透射边界的精度高,较高频率(f=10.0 Hz)时与二阶透射边界的精度相当.而且不论动载是在低频或较高频率,其数值计算结果很稳定.
本文提出了饱和多孔介质的三维黏弹性人工边界和饱和多孔介质三维动力反应分析的显式有限元法,给出计算公式,开发了实现该有限元方法的程序.数值计算表明:本文提出饱和多孔介质三维黏弹性人工边界具有较好的精度和稳定性,应用简单,便于有限元数值计算;有限元程序具有较好的精度和稳定性,运算速度快,能达到较好的数值分析效果.
[1] | Modaressi H, Benzenati I. An absorbing boundary element for dynamic analysis of two-phase media. In:Proceedings of the 10th World Conference on Earthquake Engineering, 1992. 1157~1161 |
[2] | Degrande G, De Roeck G. Absorbing boundary condition for wave propagation in saturated poroelastic media-Part I:Formulation and efficiency evaluation. Soil Dynamics and Earthquake Engineering , 1993, 12: 411-421. DOI:10.1016/0267-7261(93)90004-B |
[3] | Akiyoshi T, Fuchida K, Fang H L. Absorbing boundary conditions for dynamic analysis of fluid-saturated porous media. Soil Dynamics and Earthquake Engineering , 1994, 13: 387-397. DOI:10.1016/0267-7261(94)90009-4 |
[4] | Gajo A, Saetta A, Vitaliani R. Silent boundary conditions for wave propagation in saturated porous media. International Journal for Numerical and Analytical Methods in Geomechanics , 1996, 20: 253-273. DOI:10.1002/(ISSN)1096-9853 |
[5] | Akiyoshi T, Sun X, Fuchida K. General absorbing boundary conditions for dynamic analysis of fluid-saturated porous media. Soil Dynamics and Earthquake Engineering , 1998, 17: 397-406. DOI:10.1016/S0267-7261(98)00026-8 |
[6] | 邵秀民, 蓝志凌. 流体饱和多孔介质波动方程的有限元解法. 地球物理学报 , 2000, 43(2): 264–278. Shao X M, Lan Z L. Finite element methods for the equations of waves in fluid-saturated porous media. Chinese J. Geophys. (in Chinese) , 2000, 43(2): 264-278. |
[7] | Zerfa Z, Loret B. A viscous boundary for transient analyses of saturated porous media. Earthquake Engineering and Structural Dynamics , 2004, 33: 89-110. DOI:10.1002/(ISSN)1096-9845 |
[8] | Bowen R M. Incompressible porous media by use of the theory of mixtures. Inter. J. Eng. Sci. , 1980, 18: 19-45. DOI:10.1016/0020-7225(80)90004-X |
[9] | 刘光磊, 宋二祥. 饱和无限地基数值模拟的黏弹性传输边界. 岩土工程学报 , 2006, 28(12): 2128–2133. Liu G L, Song E X. Visco-elastic transmitting boundary for numerical analysis of infinite saturated soil foundation. Chinese Journal of Geotechnical Engineering (in Chinese) , 2006, 28(12): 2128-2133. |
[10] | 王子辉, 赵成刚, 董亮. 流体饱和多孔介质黏弹性动力人工边界. 力学学报 , 2006, 38(5): 605–611. Wang Z H, Zhao C G, Dong L. A viscous-spring dynamic artificial boundary for saturated porous media. Chinese Journal of Theoretical and Applied Mechanics (in Chinese) , 2006, 38(5): 605-611. |
[11] | 杜修力, 李立云. 饱和多孔介质近场波动分析的一种黏弹性人工边界. 地球物理学报 , 2008, 51(2): 575–581. Du X L, Li L Y. A viscous spring artificial boundary for near field wave analysis in saturated porous media. Chinese J. Geophys. (in Chinese) , 2008, 51(2): 575-581. |
[12] | 赵成刚, 王进廷, 史培新, 等. 流体饱和两相多孔介质动力反应分析的显式有限元法. 岩土工程学报 , 2001, 23(2): 178–182. Zhao C G, Wang J T, Shi P X, et al. Dynamic analysis of fluid-saturated porous media by using explicit finite element method. Chinese Journal of Geotechnical Engineering (in Chinese) , 2001, 23(2): 178-182. |
[13] | Biot M A. Theory of propagation of elastic waves in a fluid-saturated porous solid.I. Low-frequency range. The Journal of the Acoustical Society of America , 1956, 28: 168. DOI:10.1121/1.1908239 |
[14] | 廖振鹏. 工程波动理论导论. 北京: 科学出版社, 2002 . Liao Z P. Introduction to Wave Motion Theories in Engineering (in Chinese). Beijing: Science Press, 2002 . |
[15] | 王振宇.大型结构-地基系统动力反应计算理论及其应用研究.北京:清华大学, 2002 Wang Z Y.Computational theory of dynamic response of large structure-soil systems and its application.Beijing:Tsinghua University, 2002 http://www.oalib.com/references/18563986 |