横观各向同性饱和两相介质波动问题的时域显式有限元方法.pdf

第 34 卷 第 3 期 岩 土 工 程 学 报 Vol.34 No.3 2012 年 .3 月 Chinese Journal of Geotechnical Engineering Mar. 2012 横观各向同性饱和两相介质波动问题的时域显式 有限元方法 李 亮，翟 威，杜修力 北京工业大学城市与工程安全减灾教育部重点实验室，北京 100124 摘 要针对横观各向同性饱和两相介质的弹性波动问题，分别应用空间解耦技术和时域显式逐步积分计算格式完成 对该问题的有限元空间离散和时间离散，建立了横观各向同性饱和两相介质波动问题计算分析的时域显式有限元方法。 应用该方法进行了输入地震波作用下横观各向同性饱和两相介质动力反应的计算分析，并将计算结果与完全各向同性 饱和两相介质的计算结果进行对比研究，并对各向异性系数取值对于横观各向同性饱和两相介质动力反应计算结果的 影响进行了研究和讨论。计算结果表明，横观各向同性饱和两相介质与完全各向同性饱和两相介质的动力反应特性具 有较为显著的差异，各向异性系数取值对于横观各向同性饱和两相介质动力反应计算结果具有较为重要的影响。本文 的数值计算工作同时表明，本文建立的时域显式有限元方法是进行横观各向同性饱和两相介质动力反应计算分析的一 种有效方法。 关键词横观各向同性；饱和两相介质；弹性波动；时域显式有限元法 中图分类号TU435 文献标识码A 文章编号1000–4548201203–0464–07 作者简介李 亮1975– ，男，山西省太谷县人，博士，副教授，硕士生导师，主要从事土动力学理论与数值计算 方法的研究。E-mail liliang。 Time-domain explicit finite element for wave propagation of transversely isotropic fluid-saturated porous media LI Liang, ZHAI Wei, DU Xiu-li The Key Laboratory of Urban Security and Disaster Engineering Beijing University of Technology, Ministry of Education, Beijing 100124, China Abstract A time-domain explicit finite element for the elastic wave propagation of transversely isotropic fluid-saturated porous media is put forward. The space decoupling technology is adopted for the space discretization and the time-domain explicit step-by-step calculating at for the time discretization. Using the , the dynamic response of transversely isotropic fluid-saturated porous media is calculated and analyzed, and the calculated results are compared with those of isotropic fluid-saturated porous media. The effect of the value of anisotropic coefficient on the calculated results of the dynamic response of transversely isotropic fluid-saturated porous media is also studied. The calculated results show that the dynamic response of transversely isotropic fluid-saturated porous media has remarkable difference from that of isotropic fluid-saturated porous media, and the value of anisotropic coefficient has a significant effect on the calculated results of the dynamic response of transversely isotropic fluid-saturated porous media. Meanwhile, the present numerical calculation indicates that the time-domain explicit finite element is effective for the calculation and analysis of the dynamic response of transversely isotropic fluid-saturated porous media. Key words transverse isotropy; fluid-saturated porous medium; ealstic wave propagation; time-domain explicit finite element 0 引 言 饱和两相介质是一种由固相骨架与骨架孔隙中充 填的流体共同构成的二元混合体，其在自然界中广泛 ─────── 基金项目国家自然科学基金重大研究计划项目（90715035） ；国家自 然科学基金面上项目（51178011） ；北京市科技新星计划（A 类）项目 （2008A016） ；2011 年度北京市属高等学校人才强教深化计划中青年 骨干人才项目（PHR20110808） 收稿日期2011–07–29 第 3 期 李 亮，等. 横观各向同性饱和两相介质波动问题的时域显式有限元方法 465 存在，河流、海洋与湖泊底部的饱和沉积土层以及水 库大坝前方的饱和淤积泥砂层均应视为饱和两相介质 进行研究。饱和两相介质的波动问题是土动力学与岩 土地震工程学中的重要研究课题之一，具有较为重要 的理论和学术意义。同时，该问题也具有较为重要的 工程背景，诸如饱和土体与结构动力相互作用、海洋 平台与离岸结构的动力反应等工程问题中都要涉及到 饱和两相介质的波动问题。对该问题的研究开始于 Biot，他建立了饱和两相介质波动问题的数学分析模 型，推导得到各向同性线弹性饱和两相介质的波动方 程[1]，为饱和两相介质波动问题的研究奠定了基础。 描述饱和两相介质波动问题的方程为固相与液相 耦联，且时间与空间相耦合的二阶偏微分方程组，只 有针对少数具有特殊边界条件与荷载形式的简单问题 才能得到解析解，对于大多数问题必须采取数值方法 进行求解。求解两相介质波动问题的数值方法主要有 有限差分法[2]、有限单元法[3-7]、边界元法[8]，以及有 限元–无限元耦合的混合法[9]等，近年来又发展了无 网格法[10]和小波分析法[11-12]等方法。其中有限单元法 由于具有物理概念明确， 求解过程高度规范化的特点， 成为饱和两相介质波动问题求解中应用最为广泛的数 值方法。 有限单元法求解过程包括对求解系统在空间和时 间域的离散化。依据时间域离散所采用的计算格式不 同，求解两相介质波动问题的有限元法可以分为隐式 方法[3]，隐显式方法[5]及显式方法[6-7]。与隐式方法 相比，显式有限元方法实现了空间与时间域的双重解 耦，无须集成系统的总体刚度矩阵，无须求解耦联方 程组，可明显减少对计算机内存资源的占用，显著地 提高计算求解的效率，对于自由度数目巨大的复杂动 力问题的优势更为明显。 对于天然饱和土层而言，其在沉积形成过程中， 水平向与垂直向的力学性质具有显著差异，表现出明 显的各向异性性质，将其模拟为横观各向同性饱和两 相介质更接近其真实的物理力学性质。本文将以横观 各向同性饱和两相介质为研究对象，建立描述其动力 反应特性的弹性波动方程组，应用空间解耦技术实现 该方程组的空间离散， 应用中心差分法和 Newmark 常 平均加速度法相结合的时域积分计算方法实现有限元 时间离散，从而建立横观各向同性饱和两相介质弹性 波动问题计算分析的时域显式有限元方法，并应用该 方法对横观各向同性饱和两相介质在入射地震波作用 下的动力反应进行计算和分析。 1 横观各向同性饱和两相介质弹性波 动方程 1.1 运动方程 横观各向同性饱和两相介质的运动方程可由连续 介质力学理论建立，其写成矩阵形式如下 [ ] [ ]{ } T s 1{ } { }1LnpbuUnuσρ−∇ −−− ，1a f { } { }{ }n pbuUnUρ∇ − 。 1b 式中 s ρ和 f ρ分别为横观各向同性饱和两相介质中 固相和液相 （孔隙流体） 的质量密度； { }u 、{ }U 、 { }u 、 { }U 分别为固相和液相的速度与加速度向量；b为渗 流阻尼系数， 2 ff /bnkμ，其中 f μ为孔隙流体的粘滞 系数， f k为达西渗透系数，n为孔隙率。梯度算子 T xy ⎧⎫∂∂ ∇ ⎨ ⎬ ∂∂ ⎩⎭ ，微分算子矩阵[ ] T L 0 0 xy yx ∂∂⎡⎤ ⎢⎥ ∂∂ ⎢⎥ ∂∂⎢⎥ ⎢⎥ ∂∂ ⎣⎦ 。 1.2 本构关系 横观各向同性弹性饱和两相介质可视为由具有横 观各向同性性质的弹性骨架（固相）和孔隙流体（液 相）所构成的二元混合体。横观各向同性弹性体的力 学特性表现为在平行某一平面的所有各个方向上都具 有相同的弹性性质，而与之垂直的层面方向上的弹性 性质不同。横观各向同性弹性体的本构关系可表示为 111213 2223 33 44 55 66 000 000 000 00 0 xx yy zz xyxy yzyz xzxz aaa aa a a a a εσ εσ εσ γτ γτ γτ ⎧⎫⎡⎤⎧⎫ ⎪⎪⎢⎥⎪⎪ ⎪⎪⎢⎥⎪⎪ ⎪⎪⎢⎥⎪⎪ ⎪⎪⎪⎪ ⎢⎥⎨⎬⎨⎬ ⎢⎥⎪⎪⎪⎪ ⎢⎥⎪⎪⎪⎪ ⎢⎥⎪⎪⎪⎪ ⎢⎥⎪⎪⎪⎪ ⎩⎭⎣⎦⎩⎭ 对称 。2 式中 1122 1 1 aa E ， 1 12 1 a E μ −， 2 1323 2 aa E μ −， 33 2 1 a E ， 44 1 1 a G ， 5566 2 1 aa G 。其中， 1 E， 1 ν和 1 G分别为水平向的弹性模量、泊松比和剪切模量； 2 E， 2 ν和 2 G分别为竖直向的弹性模量、泊松比和剪 切模量。 土力学中研究的问题一般为平面应变问题，假设 土体在水平面上均质各向同性，水平方向和竖直方向 的不等性用各向异性系数 2 α表示，即 2 12 /E Eα。根 据文献[13]可知横观各向同性饱和土体土骨架各向异 性 系 数 的 大 致 范 围 为0.2 ≤ 2 α≤ 6.5 。 令 112 EEGGνν，，， 12 /ν να，于是针对平面问 题的横观各向同性饱和两相介质中固相骨架的刚度矩 阵改写为如下的形式 466 岩 土 工 程 学 报 2012 年 [ ] 1121 2122 33 0 0 00 kk Kkk k ⎡⎤ ⎢⎥ ⎢⎥ ⎢⎥ ⎣⎦ ， 3 式 中 ， 2 11 1E k ν λ − ， 1221 1/E kk ννα λ ， 22 22 1/E k να λ − ， 33 kG， 23 132λνν−−。 横观各向同性饱和两相介质中液相的本构方程可 写为 wvfvs 1 n pE n εε − 。 4 式中 w E为孔隙流体的体变模量； vs ε和 vf ε分别固相 骨架和孔隙流体的体积应变。 1.3 波动方程 将固相和液相的本构关系表达式（3）和（4）代 入运动方程（1） ，并引入小变形条件下的几何方程将 应变用位移表示，从而得到横观各向同性饱和两相介 质的弹性波动方程组如下 [ ] [][ ]{ }[ ] [][ ]{ }{ }{ } TT 12 LDLuLDL UbuU−− { } s 1nuρ− ， 5a [ ] [][ ]{ }[ ][][ ]{ }{ }{ } TT 23 LDLuLDL UbuU− { } f nUρ 。 5b 上述波动方程组以横观各向同性饱和两相介质的 固相位移{ }u和液相位移{ }U为基本未知量，式中的 [] 1 D，[] 2 D和[] 3 D为波动方程组的系数矩阵，其表达 式为 [] 22 11w12w 22 121w22w 33 11 0 11 0 00 nn kEkE nn nn DkEkE nn k ⎡⎤ −− ⎢⎥ ⎢⎥ ⎢⎥ −− ⎢⎥ ⎢⎥ ⎢⎥ ⎢⎥ ⎢⎥ ⎣⎦ ， .[] ww 2ww 110 110 000 n En E Dn En E −−⎡⎤ ⎢⎥ −− ⎢⎥ ⎢⎥ ⎣⎦ ， .[] ww 3ww 0 0 000 nEnE DnEnE ⎡⎤ ⎢⎥ ⎢⎥ ⎢⎥ ⎣⎦ 。 2 横观各向同性饱和两相介质弹性波 动问题的时域显式有限元方法 2.1 空间离散 采用有限元方法求解动力学问题的主要步骤可以 概括为对待求解问题在空间和时间两个方面的离散化 处理。首先进行的是空间离散，即采用一定的单元形 式将连续的计算区域划分为由若干单元和节点组成的 离散求解系统。 应用伽辽金Galerkin方法[14]， 可以得到对横观各 向同性饱和两相介质的弹性波动方程组 （5） 进行有限 元空间离散后的 Galerkin 弱式为 [ ][ ]{ }[ ][ ] { }{}[ ] [][ ]{ } [ ] [][ ]{} TT se T ee1e T 2e 1d d d dd d d d ee e e ee e e NnNux yNb N uUx yBDBux y BDB Ux y ρ ΩΩ Ω Ω − − ∑∑∫∫∫∫ ∑∫∫ ∑∫∫ [ ] [ ]{} T ue d e e NNf Γ Γ ∑∫ ， 6a [ ][ ]{}[ ][ ] { }{}[ ] [][ ]{ } [ ][][ ]{} TT fe T 2e T 3 d d d dd d d d ee e e ee ee e e e NnNUx yNb N uUx yBDBux y BDB Ux y ρ ΩΩ Ω Ω − − ∑∑∫∫∫∫ ∑∫∫ ∑∫∫ [ ] [ ]{} T Ue d e e NNf Γ Γ ∑∫ 。 6b 式中 {} ue f和{} Ue f分别为作用在单元节点上的固相 和液相边界外力，[ ]N为单元的形函数矩阵； [ ] [ ][ ]BLN，[ ][ ] [ ] TTT BNL。 依据波速有限性原理可知，在一定时间内波动的 传播范围有限，因此可认为波动对指定的某一节点的 影响范围有限。据此可实现波动问题显式有限元方法 的关键部分之一， 空间解耦技术， 即对于给定的节点， 认为其动力反应只受其周围邻近单元和节点的影响。 该节点与其周围相邻的单元和节点组成一个相对独立 的计算系统，如采用四节点矩形单元进行计算区域的 离散化，则该系统包含 4 个单元，9 个节点。对于指 定节点动力反应的计算求解，可以在上述局部单元和 节点系统内独立进行，即只考虑该节点周围的邻近节 点对其动力反应（如位移、速度和加速度等）的影响， 该系统之外的其他节点和单元的影响不予考虑。 因此， 计算过程中只需要集成反映相邻节点动力相互作用的 单元矩阵，不需要集成整个计算系统的总体动力相互 作用矩阵，从而大大减少了有限元前处理的工作量， 减少了对计算机内存资源的占用。 应用空间解耦技术可以得到关于指定节点i 动力 反应的解耦方程组为 { }{ } {} 4444 1111 LLLLL uijjijjj LjLj MuCuU − ∑∑∑∑ { }{}{} 44444 11111 LLLLL uuijjuUijjui LjLjL KuKUf ∑∑∑∑∑ ， 7a 第 3 期 李 亮，等. 横观各向同性饱和两相介质波动问题的时域显式有限元方法 467 {}{ } {} 4444 1111 LLLLL Uijjijjj LjLj MUCuU −− ∑∑∑∑ { }{}{} 44444 11111 LLLLL UuijjUUijjUi LjLjL KuKUf ∑∑∑∑∑ 。 7b 上述方程组对应于采用四节点矩形单元进行空间 离散的情况，对 L求和表示节点i 周围的 4 个单元对 其动力反应的贡献的叠加，对j求和表示同一单元内 的 4 个节点对节点i 动力反应的贡献的叠加。 式（7）中各个矩阵元素的表达式如下 11 11 11 11 [ , ][ , ]d d [ , ][ , ]d d L L uijij ij MN x yNx yx y NNJ ρ ρξ ηξ ηξ η Ω −− ∫∫ ∫ ∫ ， 22 11 22 11 [ , ][ , ]d d [ , ][ , ]d d L L Uijij ij MN x yNx yx y NNJ ρ ρξ ηξ ηξ η Ω −− ∫∫ ∫ ∫ ， 11 11 [ , ][ , ]d d [ , ][ , ]d d L L ijij ij Cb N x yNx yx y b NNJξ ηξ ηξ η Ω −− ∫∫ ∫ ∫ ， 1 11 1 11 [ , ] [][ , ]d d [ , ] [][ , ]d d L LT uuijij T ij KB x yDB x yx y BDBJξ ηξ ηξ η Ω −− ∫∫ ∫ ∫ ， 2 11 2 11 [ , ] [][ , ]d d [ , ] [][ , ]d d L LLT uUijUuijij T ij KKB x yDB x yx y BDBJξ ηξ ηξ η Ω −− ∫∫ ∫ ∫ ， 3 11 3 11 [ , ] [][ , ]d d [ , ] [][ , ]d d L LT UUijij T ij KB x yDB x yx y BDBJξ ηξ ηξ η Ω −− ∫∫ ∫ ∫ ， 式中， 11s 1nρρ−， 22f nρρ。 2.2 时间离散 有限元时间离散即建立计算系统中各节点动力反 应的时域逐步积分计算格式。忽略单元内部惯性力的 变化，假定在同一单元内惯性力为常量，于是可以得 到集中质量模型，从而可以实现节点动力反应计算格 式的显式化，即时间域解耦。应用中心差分法和 Newmark 常平均加速度法相结合的数值积分算法，建 立了横观各向同性饱和两相介质动力反应的时域逐步 积分计算格式，其表达式为 {} { }{ } {}{} {} 1 2 444 1,,, 111 2 ppp iii L pLL pL p uiuiijjj LLj uut u t MfCuU − Δ ⎡Δ −−− ⎢ ⎣∑ ∑∑ {}{} 4444 ,, 1111 LL pLL p uuijjuUijj LjLj KuKU ⎤ − ⎥ ⎦ ∑∑∑∑ ， 8a {} {}{} {}{} {} 1 2 444 1,,, 111 2 ppp iii L pLL pL p UiUiijjj LLj UUt U t MfCuU − Δ ⎡Δ −− ⎢ ⎣∑ ∑∑ {}{} 4444 ,, 1111 LL pLL p UuijjUUijj LjLj KuKU ⎤ − ⎥ ⎦ ∑∑∑∑ ， 8b {} { }{} {} {} {}{} {} 44 11,1, 11 4 ,1,,1, 1 2 ppLL pL p iiuiijjj Lj L pL pL pL p jjuiui L uuMCuu t UUff − ⎧ ⎡ −− ⎨ ⎣ ⎩ Δ ⎡ ⎤ −− ⎢ ⎦ ⎣ ∑∑ ∑ {} {}{} {} 4444 ,1,,1, 1111 LL pL pLL pL p uuijjjuUijjj LjLj KuuKUU ⎫⎤⎪ − ⎬⎥⎪ ⎦⎭ ∑∑∑∑ ， 9a {} {}{} {} {} {}{} {} 44 11,1, 11 4 ,1,,1, 1 2 ppLL pL p iiUiijjj Lj L pL pL pL p jjUiUi L UUMCuu t UUff − ⎧ ⎡ −− ⎨ ⎣ ⎩ Δ ⎡ ⎤ −− ⎢ ⎦ ⎣ ∑∑ ∑ {} {}{} {} 4444 ,1,,1, 1111 LL pL pLL pL p UuijjjUUijjj LjLj KuuKUU ⎫⎤⎪ − ⎬⎥⎪ ⎦⎭ ∑∑∑∑ 。 9b 由上述表达式可以看出，横观各向同性饱和两相 介质在当前步（第1p 步）的动力反应可以按照上述 计算格式由前一步（第p步）的动力反应直接计算得 到。表达式（8） 、 （9）共同组成横观各向同性饱和两 相介质动力反应的时域显式逐步积分计算格式，其类 似于一种显式差分计算格式。按照上述计算格式进行 逐步的迭代计算，就可以完成横观各向同性饱和两相 介质动力反应时程的求解，无须在计算的每一个时间 分步上求解耦联方程组。 对于自由度数目巨大的问题， 应用该显式计算格式将能显著地缩短计算时间，提高 计算效率。这正是显式计算方法的优势所在。 3 算 例 应用本文建立的时域显式有限元方法，对横观各 向同性饱和两相介质在入射地震波作用下的动力反应 进行计算和分析。算例的计算模型如图1所示，横观 图 1 横观各向同性饱和两相介质动力反应计算模型 Fig. 1 Calculating model for dynamic response of transversely isotropic fluid-saturated porous media 468 岩 土 工 程 学 报 2012 年 各向同性饱和两相介质计算区域取为200 m200 m， 单元网格划分尺寸为10 m10 m，共划分得到400个 单元和441个节点。采用Loma Prieta earthquake的地 震波作为输入地震波进行计算，其持时为48.69 s。输 入地震波的位移时程和速度时程分别如图2和图3所 示。迭代计算的时间步长取为0.003 s，对应于完整的 地震波输入，迭代计算共分16230步进行。计算所采 图 2 输入地震波位移时程 Fig. 2 Displacement time history of earthquake waves 图 3 输入地震波速度时程 Fig. 3 Velocity time history of earthquake waves 用的主要输入参数列于表1中。 对于横观各向同性饱和两相介质而言，可以采用 不同的各向异性系数 2 α的取值描述其各向异性程度 的不同。本文分别取各向异性系数 2 α的值为0.5，1.0 和5.0进行横观各向同性饱和两相介质动力反应的计 算， 其中 2 α取1.0对应于完全各向同性的情况。 2 α取 值不同时得到的输入地震波作用下横观各向同性饱和 两相介质计算区域自由面的固相和液相位移反应时程 分别如图4和图5所示，固相和液相速度反应时程分 别如图6和图7所示。由图中给出的计算结果可知， 横观各向同性饱和两相介质与完全各向同性饱和两相 介质的动力反应有较为显著的差异，这体现在反应峰 值的差异与动力反应时程波形的差异两个方面，与后 者相比，前者在峰值前后的动力反应幅值均有较为明 显的变化。同时，各向异性系数 2 α的取值对于横观各 向同性饱和两相介质动力反应的峰值影响显著，随着 2 α取值的增大，反应峰值逐渐增大。 2 α取值不同时 得到的横观各向同性饱和两相介质动力反应的峰值列 于表2中，由图中数据可知，与各向异性系数 2 α的取 值为0.5时的计算结果相比， 2 α的取值为5.0时位移 反应的正向与负向峰值均增大了约58，速度反应的 正向峰值增大了5.2，负向峰值增大了22.4。 对于本文所采用的横观各向同性本构关系而言， 可以通过调整各向异性系数 2 α的取值， 来描述具有不 同各向异性程度的饱和两相介质，实现对具有不同程 度各向异性的饱和两相介质动力反应的计算和分析。 4 结 论 本文建立了描述横观各向同性饱和两相介质弹性 波动问题的方程组，分别应用空间解耦技术和时域显 表 1 横观各向同性饱和两相介质动力反应计算参数 Table1 parameters for dynamic response calculation of transversely isotropic saturated media s ρ /kgm -3 f ρ /kgm -3 E /107 Pa G /107 Pa w E /109 Pa n ν 2 α 2600 1000 3.3 1.2 2.0 0.3 0.3 0.5 1.0 5.0 表 2 各向异性系数取值不同时动力反应峰值对比 Table 2 Peak values of dynamic response corresponding to different anisotropic coefficients 各向异性系数 2 α 固相位移/m 液相位移/m 固相速度/ms -1 液相速度/ms -1 正向峰值 0.141 0.141 0.406 0.406 0.5 负向峰值 0.162 0.162 0.540 0.540 正向峰值 0.170 0.170 0.387 0.387 1.0 负向峰值 0.195 0.195 0593 0593 正向峰值 0.223 0.223 0.427 0.427 5.0 负向峰值 0.257 0.257 0.661 0.661 第 3 期 李 亮，等. 横观各向同性饱和两相介质波动问题的时域显式有限元方法 469 图 4 计算区域自由面固相位移反应时程 Fig. 4 Solid-phase displacement time history of free surface 图 5 计算区域自由面液相位移反应时程 Fig. 5 Liquid-phase displacement time history of free surface 图 6 计算区域自由面固相速度反应时程 Fig. 6 Solid-phase velocity time history of free surface 图 7 计算区域自由面液相速度反应时程 Fig. 7 Liquid-phase velocity time history of free surface 式逐步积分计算格式完成对该问题的有限元空间离散 和时间离散，建立了横观各向同性饱和两相介质弹性 波动问题计算分析的时域显式有限元方法。应用该方 法进行了输入地震波作用下横观各向同性饱和两相介 质动力反应的计算分析，并将计算结果与完全各向同 性饱和两相介质的计算结果进行对比研究，并对各向 异性系数取值对于横观各向同性饱和两相介质动力反 应计算结果的影响进行了研究和讨论。 研究结果表明 （1） 横观各向同性饱和两相介质与完全各向同性 饱和两相介质的动力反应特性具有较为显著的差异， 表现为二者动力反应峰值与动力反应时程波形的差异 两个方面。 （2） 各向异性系数取值对于横观各向同性饱和两 相介质动力反应计算结果具有较为重要的影响。随着 各向异性系数取值的增大，横观各向同性饱和两相介 质动力反应的峰值逐渐增大。 本文的计算工作同时表明，本文所建立的时域显 式有限元方法能够对横观各向同性饱和两相介质在入 射地震波作用下的动力反应进行较为准确的计算，是 进行横观各向同性饱和两相介质动力反应计算分析的 一种有效方法。 参考文献 [1] BIOT M A. Theory of propagation of elastic waves in a fluid-saturated porous solid[J]. The Journal of the Acoustical Society of America, 1956, 282 168–191. [2] GARG S K. Wave propagation effects in a fluid-saturated porous solid[J]. Journal of Geophysical Research, 1971, 7632 7947–7962. [3] ZIENKIEWICZ O C, SHIOMI T. Dynamic behaviour of saturated porous media; the generalized Biot ulation and its numerical solution[J]. International Journal for Numerical and Analytical s in Geomechanics, 1984, 81 71– 96. [4] SIMON B R, WU J S S, ZIENKIEWICZ O C. uation of higher order, mixed and Hermitean finite element procedures for dynamic analysis of saturated porous media using one-dimensional models[J]. International Journal for Numerical and Analytical s in Geomechanics, 1986, 105 483–499. [5] PREVOST J H. Wave propagation in fluid-saturated porous media an efficient finite element procedure[J]. Soil Dynamics and Earthquake Engineering, 1985, 44 183– 202. [6] 赵成刚, 王进廷, 史培新, 等. 流体饱和两相多孔介质动力 反应分析的显式有限元法[J]. 岩土工程学报, 2001, 232 178–182. ZHAO Cheng-gang, WANG Jin-ting, SHI Pei-xin, et al. Dynamic analysis of fluid-saturated porous media by using explicit finite element [J]. Chinese Journal of 470 岩 土 工 程 学 报 2012 年 Geotechnical Engineering, 2001, 232 178 – 182. in Chinese [7] 王进廷, 杜修力, 赵成刚. 液固两相饱和介质动力分析的 一种显式有限元法[J]. 岩石力学与工程学报, 2002, 218 1199 – 1204. WANG Jin-ting, DU Xiu-li, ZHAO Cheng-gang. Explicit finite element for dynamic analysis of fluid-saturated porous solid[J]. Chinese Journal of Rock Mechanics and Engineering, 2002, 218 1199–1204. in Chinese [8] DOMINGUEZ J. Boundary element approach for dynamic poroelastic problems[J]. International Journal for Numerical s in Engineering, 1992, 352 307–324. [9] KHALILI N, YAZDCHI M, VALLIAPPAN S. Wave propagation analysis of two-phase saturated porous media using coupled finite-infinite element [J]. Soil Dynamics and Earthquake Engineering, 1999, 188 533– 553. [10] 张洪武, 王组鹏, 陈 震. 基于物质点方法饱和多孔介质 动力学模拟（I）耦合物质点方法[J]. 岩土工程学报, 2009, 3110 1505–1511. ZHANG Hong-wu, WANG Kun-peng, CHEN Zhen. Material point for dynamic analysis of saturated porous mediaI coupling material point [J]. Chinese Journal of Geotechnical Engineering, 2009, 3110 1505–1511. in Chinese [11] 张新明, 刘克安, 刘家琦. 流体饱和多孔隙介质二维弹性 波方程正演模拟的小波有限元法[J]. 地球物理学报, 2005, 485 1156–1166. ZHANG