基于DQEM的分层流体饱和热弹性多孔介质轴对称问题的动力响应分析

叶东生 朱媛媛 王笑梅 王玉善

摘 要： 研究了分层不可压流体饱和热弹性多孔介质轴对称问题的动力响应问题，基于多孔介质理论（PMT），给出了该问题的数学模型.在空间域内采用微分求积单元法（DQEM）设置离散的控制微分方程、边界条件和连接条件，在时间域内采用二阶向后差分格式处理时间导数.在离散化的初始条件下，运用Newton-Raphson法进行迭代求解，得到各离散点处未知物理量的数值结果.研究表明：该方法有效、可靠，且具备精度较高、计算量较小、数值稳定等优点.

关键词： 流体饱和热弹性多孔介质; 多孔介质理论（PMT）; 微分求积单元法（DQEM）; 动力学响应

中图分类号： TU 311文献标志码： A文章编号： 1000-5137（2019）02-0141-10

0 引 言

饱和多孔介质热-流-固耦合系统的研究不仅在土力学、水文学等经典应用领域发挥重要作用，也已成为许多新兴学科和应用技术发展的关键，其相关理论和数值方法的研究工作具有重要的理论意义和广泛的应用背景.

1955年，BIOT[1]建立了饱和多孔介质热弹性和热动力理论.随后，CARTER等[2]研究了饱和土球形热源附近的固结问题.CUI等[3]引入热孔隙度状态表面的概念，提出了饱和土热-水-力耦合分析的理论模型.刘干斌等[4]通过对Biot波动方程的修正，研究了简谐均布荷载作用下，多孔弹性地基土体的热-水-力耦合动力响应问题.白冰[5]对循环温度荷载下饱和多孔介质热-水-力耦合响应的一维情形进行了研究.戴清晨等[6]研究了热局部非平衡条件下，横观各向同性饱和多孔介质中柱形空洞的热应力.de BOER等[7-8]基于连续介质混和物公理体系和体积分数概念，建立了较完整的多孔介质理论.de BOER等[9]利用拉普拉斯变换，给出了多孔介质一维动力学响应的解析解.刘占芳等[10]等利用拉普拉斯变换和卷积定理，得到了边界自由排水时，任意应力和位移边界条件下，瞬态波动过程的解析表达.李向约等[11]推导了描述饱和多孔介质在热固结耦合作用下的数学模型.HE等[12]建立了热局部非平衡条件下饱和不可压多孔弹性介质热-流-固耦合模型.YANG[13]建立了相应的Gurtin型广义变分原理.朱媛媛等[14]利用微分求积法（DQM）分析了流体饱和热弹性多孔介质圆柱体的动力响应.严俊等[15]利用连续介质理论，提出了非饱和多孔介质的热本构关系以及孔隙流体的热运动规律.

本文作者分析了分层流体饱和热弹性多孔介质的动力学行为，基于多孔介质理论，建立了一维流体饱和热弹性多孔介质轴对称问题的数学模型，采用微分求积单元法（DQEM）和微分-代数方法求解耦合系统的动力学问题，获得耦合系统一些定性的结果.将DQEM用于求解分层流体饱和热弹性问题分析中，研究工作有一定的探索性，建模思路和研究方法可为工程实践提供参考.

1 问题的数学描述

1.1 控制微分方程

图1为由n层介质组成的分层轴对称不可压流体饱和多孔热弹性体.其中，（x，y）表示直角坐标系，（r，θ）为极坐标系，n表示第n层介质，q（t）表示外加垂直机械载荷，φ（t）表示外加温度载荷.

式（1）中第1，2个方程分别为第i层固相和流相介质的动量守恒方程，第3个方程为质量守恒方程，第4个方程为能量守恒方程;uir为第i层固相的位移分量，u·ir为其相应的速度分量，u··ir为其相应的加速度分量;wir为第i层流相相对速度，w·ir为其相应的相对加速度;θi=Θi-Θi0为第i层系统的变温，其中Θi为第i层系统的绝对温度，Θi0为第i层系统的初始绝对温度，θ·i表示第i层温度的变化速度; σSEir，σSEiθ分别表示第i层系统极径r上的固相有效应力和极角θ的固相有效应力，pi为第i层流体的有效孔隙压力，nαi（α=S，F）分别代表固相和流相的体积分数，且nSi+nFi=1，ραi （α=S，F）分别代表固相和流相的宏观质量密度，它与微观质量密度ραRi之间有关系ραi=nαiραRi，在两相不可压缩的情况下ραRi为常数;KSi=λSi+2μSi/3为体积模量，μSi，λSi为固相材料的Lame系数，αSi为固相材料的热膨胀系数;Siv=（niF）2γFRiκFi为固相与流相间的耦合系数，其中，γFRi为流相的有效比重，κFi为达西渗透系数;βi为与流速有关的附加热交换因子;ki为物相之间的热传导系数;ρic=ρSicSi+ρFicFi为第i层系统的总比热系数，其中，cαi （α=S，F）分别代表固相和流相的比热系数;ρir=ρSirSi+ρFirFi为第i层系统的总热源强度，其中，rαi（α=S，F）分別代表固相和流相的热源强度.

1.3 界面之间的连接性条件

对于分层不可压流体饱和多孔热弹性体而言，必须满足界面之间的连接性条件.本研究的连接性条件为：1） 界面之间的固相位移分量必须相等;2） 界面之间的总应力分量必须相等;3） 界面之间的流体有效孔隙压必须相等;4） 界面之间的流体流量必须相等;5） 界面之间的变温必须相等;6） 界面之间的热流强度必须相等.

1.4 初始条件

2 DQEM和控制方程的微分求积（DQ）离散化

DQM是BELLMAN等[16-17]在20世纪70年代初提出的一种求解偏微分方程的数值算法.该算法的基本思路是用解区域中所有离散点处沿某个方向函数值的线性加权和作为该未知函数和它的各阶导数在某一离散点的近似值，权系数只与解区域中所选择的离散点和试函数有关.因此，任何一个微分方程都可以转化成相应的代数方程.

DQM具有公式简单、使用方便、计算量少、精度高等优点.传统的DQM对于求解具有非规则区域和间断性条件的问题，存在一些局限性.因此，研究人员构建了微分求积单元法（DQEM），并取得了一系列的研究成果[18-19].DQEM基本步骤是：1） 将求解区域分割成若干个子区域或单元;2） 利用DQM，将各子区域的微分方程和边界条件转化为离散的代数方程组或者常微分方程组;3）将各单元的离散化方程连接起来，组成一个整体离散化的代数方程组或者常微分方程组;4） 采用适当方法求解，从而得到各节点的未知量.

考虑在区域Ω={x0≤x≤a}内的未知函数ψ（x），设沿x方向布置Nx个节点，根据DQM，函数ψ（x）在节点x=xξ（ζ=1，2，3，…，Nx）处对自变量x的n阶导数可近似表示为：

其中，ψk=ψ（xk），为相应节点的函数值，A（n）ζk为试函数对于x的n阶偏导数的权系数.本研究中，权系数由Lagrange插值多项式决定.

按照分层不可压流体饱和多孔热弹性体，耦合系统被划分为n个单元，每个单元内布置N个节点（图2），节点坐标由Chebyshev-Lobatto多项式的零点决定.

2.1 空间域内控制方程的DQ离散化

2.2 边界条件和对称性条件的DQ离散化

2.4 对称轴上奇异性条件的处理

3 分层不可压轴对称流体饱和多孔热弹性体的动力学特性

3.1 数值结果的验证

耦合系统被分别划分为2，3，4个单元，每个单元内布置Ni=11个节点（图2），节点坐标由Chebyshev-Lobatto多项式的零点来决定.

图3给出了热弹性体不同深度处的位移ur曲线，Δt=1 s.其中，ra表示饱和多孔介质半径.实线和圈划线分别对应n=2个单元和n=4个单元的情形下，利用DQEM得到的数值解;点线为利用文献[14]中DQM模型得到的结果.

图4给出了二阶向后差分格式的步长对热弹性体不同深度处位移ur的影响，单元数n=3.实线和圈划线分别对应Δt=1 s和Δt=2 s的情形下，利用DQEM得到的数值解.点线为利用文献[14]中DQM模型得到的结果.

通过计算发现：每个单元内布置Ni=7个节点，能得到令人满意的结果.为了节省篇幅，不在此给出示例.

从图3和4中可以看到：采用两种模型求得的解趋于一致，证明本方法具有较高的精度和收敛性.

3.2 分层不可压轴对称流体饱和多孔热弹性体的动力学特性

3.2.1 热交换系数βi对热弹性体动力学特性的影响

从图5中可以看到：ur随着时间的增加趋于相同稳定值;wr随时间增加逐渐趋于0;p由初始值逐渐消散至0，且表面附近孔隙压的消散速度快于内部的消散速度;θ随时间的增加逐渐上升并由表面向纵深处传导和扩散，最后达到等温状态.同时，当βi较小时，流、固两相之间相互作用力较小，初始阶段固相的热体积膨胀效应被抑制，热弹性体沉降大.当βi较大时，热传导过程快于机械载荷下的固结过程，初始阶段表现为固相的热体积膨胀效应，而后固结作用才逐渐显现.

从图6中可以看到，在分层不可压轴对称流体饱和多孔热弹性体中，由于各个层中βi不同，p和θ在界面处不连续.

3.2.2 体积分数对热弹性体动力学特性的影响

图7为不同的体积分数ci对分层热弹性体动力学特性的影响.实线为分层（n=3）多孔热弹性体的实验结果，此时c1=0.6，c2=0.8，c3=0.6;虚线为均匀多孔热弹性体的实验结果，即c1=c2=c3=0.6.从图7中可以看到，在分层不可压轴对称流体饱和多孔热弹性体中，由于各个层中体积分数不同，wr在界面处明显不连续.

4 结 论

在热局部平衡条件下，基于PMT，研究了分层轴对称流体饱和多孔热弹性体在表面温度载荷作用下的动力学特性，提出了问题的数学模型，采用DQEM、二阶向后差分法及Newton-Raphson迭代法模拟问题的数值结果.为了验证本方法的正确性，研究了不可压流体饱和多孔弹性体的动力固结问题，并与

现有结果进行比较，二者能良好地吻合，证明DQEM具备精度高、计算量小、数值稳定等优点.研究和比较了一维分层轴对称流体饱和多孔热弹性体在表面受到温度载荷时的动力学特性，考察了材料参数对热弹性体动力学特性的影响.

参考文献：

[1] BIOT M A.Theory of elasticity and consolidation for a porous anisotropic solid [J].Journal of Applied Physics，1955，26：182-185.

[2] CARTER J R，SAVVIDOU C.Consolidation around a spherical heat source [J].International Journal of Solids and Structures，1984，20：1079-1090.

[3] CUI Y J，SULTAN N，DELAGE P.A thermomechanical model for saturated clays [J].Canadian Geotechnical Journal，2000，37：607-620.

[4] 劉干斌，姚海林，杨洋，等.考虑热-水-力耦合效应多孔弹性地基的动力响应 [J].岩土力学，2007，28（9）：1784-1795.

LIU G B，YAO H L，YANG Y，et al.Coupling thermo-hydro-mechanical dynamic response of a porous elastic medium [J].Rock and Soil Mechanics，2007，28（9）：1784-1795.

[5] 白冰.循环温度荷载作用下饱和多孔介质热-水-力耦合响应 [J].工程力学，2007，24（5）：87-92.

BAI B.Thermo-hydro-mechanical response of saturated porous media under cycle thermal loading [J].Engineering Mechanics，2007，24（5）：87-92.

[6] 戴清晨，何录武.热局部非平衡条件下含柱形空洞横观各向同性饱和多孔介质的热应力分析 [J].力学季刊，2014，35（1）：1-9.

DAI Q C，HE L W.Thermal stresses around a cylindrical hole in a transversely isotropic poroelastic medium considering local thermal non-equilibrium [J].Chinese Quarterly of Mechanics，2014，35（1）：1-9.

[7] de BOER R.Theoretical poroelasticity：a new approach [J].Chaos Solitons & Fractals，2005，25（4）：861-878.

[8] de BOER R，KOWALSKI S J.Thermodynamics of fluid-saturated porous media with a phase change [J].Acta Mechanica，1995，109（1/2/3/4）：167-189.

[9] de BOER R，EHLERS W，LIU Z.One-dimensional transient wave propagation in fluid-saturated incompressible porous media [J].Archive of Applied Mechanics，1993，63（1）：59-72.

[10] 劉占芳，姜乃斌，李思平.饱和多孔介质一维瞬态波动问题的解析分析 [J].工程力学，2006，23（7）：19-24.

LIU Z F，JIANG N B，LI S P.An analysis on one-dimensional transient wave motion in saturated porous media [J].Engineering Mechanics，2006，23（7）：19-24.

[11] 李向约，李向维.饱和多孔介质的热固结理论 [J].固体力学学报，1990，11（4）：330-338.

LI X Y，LI X W.Theory of thermo-consolidation for saturated porous elastic media [J].Acta Mechanica Solida Sinica，1990，11（4）：330-338.

[12] HE L W，JIN Z H.A local thermal non-equilibrium poroelastic theory for fluid saturated [J].Journal of Thermal Stresses，2010，33：799-813.

[13] YANG X.Gurtin-type variational principles for dynamics of a non-local thermal equilibrium saturated porous medium [J].Acta Mechanica Solida Sinica，2005，18（1）：37-45.

[14] 朱媛媛，胡育佳，程昌钧，等.基于DQM的空间轴对称流体饱和多孔热弹性柱体动力学特性研究 [J].振动与冲击，2017，36（23）：83-91.

ZHU Y Y，HU Y J，CHENG C J，et al.The study on dynamic characteristics for a space-axisymmetrical fluid-saturated porous thermo-elastic cylinder based on DQM [J].Journal of Vibration and Shock，2017，36（23）：83-91.

[15] 严俊，魏迎奇，蔡红，等.非饱和多孔介质水-热-力耦合数学模型研究 [J].水利学报，2014（增刊2）：152-160.

YAN J，WEI Y Q，CAI H，et al.A mathematical thermal hydraulic-mechanical coupling model for unsaturated porous media [J].Journal of Hydraulic Engineering，2014（Suppl.2）：152-160.

[16] BELLMAN R E，CASTI J.Differential quadrature and long term integration [J].Journal of Mathematical Analysis & Applications，1970，34（2）：235-238.

[17] BELLMAM R E，KASHEF B G，CASTI J.Differential quadrature：a technique for the rapid solution of nonlinear partial differential equations [J].Journal of Computational Physics，1972，10（1）：40-52.

[18] ZHU Y Y，HU Y J，CHENG C J.DQEM for analyzing dynamic characteristics of layered fluid-saturated porous elastic media [J].Acta Mechanica，2013，224（9）：1977-1998.

[19] 聂国隽，仲政.用微分求积法求解梁的弹塑性问题 [J].工程力学，2005，22（1）：59-62.

NIE G J，ZHONG Z.Elasto-plastic analysis of beams by differential quadrature method [J].Engineering Mechanics，2005，22（1）：59-62.

（责任编辑：包震宇）