流体饱和多孔热弹性对称平面的动力学分析

朱媛媛,杨 骁,吴海涛

(1.上海师范大学信息与机电工程学院,上海 200234;2.上海大学力学与工程科学学院,上海 200444)

饱和多孔介质热-流-固耦合系统由于具有特殊的结构和材料性质,不仅在传统的应用领域(如土力学、水文学等)得到广泛应用,而且也逐渐成为许多新兴学科和高新技术(如核废料污染物处置、地热资源、热能贮存、人体关节软骨组织分析等)发展的关键.因此,研究饱和多孔介质热力学性能的理论、数值方法及其应用具有重要的意义.

基于Biot[1]的研究工作建立的有关饱和多孔介质的热弹性理论,目前已有许多研究成果.在连续介质混和物公理体系和体积分数概念的基础上,de Boer[2]建立了较完整的多孔介质热弹性理论,该理论可以直接将一些微观性质用来描述宏观性质,而且容易反映非线性效应.利用混合物理论,de Boer等[3]为一维流体饱和不可压缩多孔介质的动力学固结问题提供了解析解.Heider等[4]研究了波在饱和多孔半空间中的传播.Hu等[5]利用微分求积法对黏弹性流体饱和多孔介质的动力学特性进行了研究.刘林超等[6]讨论了饱和土中群桩的竖向动力相互作用问题.在de Boer等[7]提出的多孔介质热力学本构关系的基础上,He等[8]给出了一种在热局部非平衡条件下的流体多孔介质的数学模型,Yang[9]建立了该模型的相应的Gurtin型变分原理,朱媛媛等[10]对空间轴对称流体饱和多孔热弹性柱体动力学特性进行了分析.但是,目前还较少看到有关饱和多孔介质热-流-固耦合系统非线性理论和应用以及数值模拟的报道.

本工作基于多孔介质理论(porous media theory,PMT)[2,7],研究了在有限变形和热局部平衡条件下,不可压流体饱和多孔热弹性半平面受到表面温度载荷作用下的动力学特性,并在建立问题非线性数学模型的基础上,通过采用微分求积(differential quadrature,DQ)法在空间域内离散问题的非线性数学模型,二阶后向差分格式来处理时间导数,Newton-Raphson法求解非线性代数方程组,从而可得到问题的数值结果,分析半平面的动力学特性.为了验证本方法的正确性,在无热效应和忽略几何非线性影响时,计算了不可压流体饱和多孔弹性介质一维动力学问题.与解析解[2]对比发现,本方法是有效可靠的,且具有计算量小、精度高等优点.最后,研究了在受到表面温度载荷作用下,流体饱和多孔热弹性半平面的动力学特性,考察了达西渗透系数、热传导系数、热交换系数等材料参数以及几何非线性的影响.

1 数学模型

图1为平面直角坐标系中流体饱和不可压多孔热弹性半平面动力学模型.设半平面处于热局部平衡状态(固相介质和液相介质的变温相同).不失一般性,假设作用于半平面表面的载荷或者温度具有某种对称性,取Oz轴为对称轴,其正方向为由表面指向半平面内部(见图1(a)),可认为沿Oy轴方向的固相位移uy≈0和流相的相对流速wy≈0,并且各物理量均与坐标y无关.假设平面介质为各向同性线性热弹性多孔固相骨架和理想流体,根据多孔介质理论,在不可压和小流速的假设下,满足如下动力学方程.

图1 流体饱和多孔热弹性半平面动力学问题模型Fig.1 Model of fluid-saturated porous thermo-elastic half-plane problem

质量守恒方程为

动量和动量矩守恒方程为

质量守恒方程为

式中:ux,uz为固相位移;为固相有效应力;wx,wz为相对流速,p为有效孔隙压,θ为变温;物相φα(α=S,F,分别表示固相及液相)体积分数nα(nS+nF=1);宏观质量密度ρα(和微观质量密度ραR的关系为ρα=nαραR);流固耦合系数Sv=(nF)2γFR/κF,其中γFR为液相有效比重,κF为达西渗透系数;热交换因子β;热传导系数k;初始绝对温度Θ0;固相Lame系数μS,λS,体积模量KSλS+2μS/3;热膨胀系数αS;比热系数cα和热源强度rα.记ρc=ρScS+ρFcF,ρr=ρSrS+ρFrF.

对于各向同性线性热弹性固相,本构关系为

在有限变形的情况下,流体饱和多孔热弹性介质固相材料的几何关系可以表示为

假设流体饱和多孔弹性半平面的表面处于理想排水状态,同时承受变温φ(x,t)和垂直载荷q(x,t)的作用;底部为不排水刚性绝热底面;侧面不排水绝热且在x-方向位移被约束(见图1(b)).因为φ(x,t)和q(x,t)关于Oz轴对称,可以在区域0≤x≤L0,0≤z≤H0内考虑问题(见图1(b)),这时的边界条件为

假设t≤0时半平面处于静止和热力平衡状态,初始条件可以写成

因此,式(1)～(6),边界条件(7)和初始条件(8)构成了在热局部平衡条件下流体饱和多孔热弹性半平面动力学问题的一个数学模型,其中基本未知量有:固相位移ux、uz;固相有效应力相对流速wx、wz;有效孔隙压p和变温θ.

可以看出,上述动力学问题的解析或半解析解很难获得.在本工作中,上述的非线性数学模型在空间域中使用DQM离散,时间导数采用二阶后向差分格式处理,离散后的非线性代数方程组采用Newton-Raphson法迭代求解,可得到问题的数值模拟结果.

2 DQM和控制方程的离散化

Bellman等[11-12]提出了DQM,该方法具有公式简单、精度高、计算量少等优点,已有很多成功的应用[13-15].DQM的基本思想是将函数的偏导数近似为离散点处相应函数值的加权和.由于DQM权函数的选取与特定问题无关,微分方程能微分求积(differential quadrature,DQ)离散为相应的代数方程求解.下面运用DQM来对数学模型进行空间离散.

假设半平面占有区域Ω={〈x,z〉|0≤x≤L0,0≤z≤H0},分别在x、z方向布置Nx×Nz个节点(见图2),节点坐标为Chebyshev-Lobatto多项式的零点.

图2 半平面上的节点布置Fig.2 Nodes collocated on the half-plane

根据DQM,未知函数ψ(x,z)对x的n阶偏导数可近似表示为

式中:ψkη是在x=xk,z=zη节点处的函数值;是n阶偏导数的权系数,本工作中由Lagrange插值多项式来决定.同时,可给出ψ(x,z)对z的m阶偏导数,权系数为

在热局部平衡条件下,利用式(9),xz平面内的基本动力学方程(式(1)～(3))的DQ离散化方程为

式中:ξ=2,3,···,Nx-1;η=2,3,···,Nz-1.

本构关系(式(4))、几何关系(式(5))和固相的有效应力分量与总应力分量关系(式(6))的DQ离散化形式为

式中:ξ=1,2,···,Nx;η=1,2,···,Nz.

相应的,边界条件和对称性条件的DQ离散化形式为

式中:qξ1,φξ1是半平面表面处给定的外载荷q(x,t)和温度载荷φ(x,t)的值.

在空间域内对半平面问题的数学模型DQM离散化后,产生了关于时间t的微分代数方程组.使用二阶向后差分格式来处理关于微分代数方程组中的时间导数

式中:Δt是时间步长;φ(tn,x)是在时刻t=tn的函数值.

利用式(13),时间离散DQ离散化方程组(10)～(14).最后利用DQ离散的初始条件式(8),对离散化方程组进行Newton-Raphson迭代求解,从而可得到问题的数值结果.

3 动力学特性

3.1 验证方法

忽略体积力的影响,Boer等[3]分析了几何小变形和常温条件下不可压饱和多孔弹性介质的一维动力固结问题,并给出了解析解.为了验证本方法的正确性、有效性,在给定的数学模型中设平面表面变温为φ(x,t)=0,从而可以认为热局部平衡条件下的解近似等于常温状态时的解.本工作用一个高H0=10 m,宽L0=10 m的流体饱和多孔弹性体来模拟Boer等[6]所考虑的类似问题,并与解析解进行比较.在计算验证中,分别考虑了阶梯载荷q(x,t)=q0h(t)和周期循环载荷q(x,t)=q0[1-cos(ωt)]的影响,其中q0=3 kN/m2,ω=0.4 s-1,h(t)为Heaviside函数.表1是计算时选取的材料参数.

表1 材料参数[13]Table 1 Material parameters[13]

图3分别给出了阶梯载荷和周期载荷下对称轴(x=0)处不同深度的位移uz(沉降),其中时间步长Δt=1 s.从图3中可以看到:数值解与解析解吻合良好;当布置结点数为7×7时,便能得到令人满意的结果.这说明用本方法计算流体饱和多孔弹性介质的动力学特性是有效的.

图3 不同深度的位移时程曲线(x=0)Fig.3 Time-history curves of the displacement in differen depths(x=0)

3.2 参数和非线性影响

3.2.1 阶梯温度载荷下的动力学特性

考察高H0=10 m,宽L0=10 m的流体饱和多孔热弹性平面域,边界条件由式(7)给定.考察表面的阶梯温度载荷φ(x,t)=φ0h(t)对半平面动力学特性的影响,其中φ0=10°C,h(t)是Heaviside函数.计算中布点数为7×7,时间步长Δt=1 s,平面表面垂直外载q(x,t)=3 kN/m2,材料参数见表1.

图4阶梯温度载荷下κF对半平面时-程曲线的影响Fig.4 The effect ofκF on time-history curves of half-plane under the step temperature loading

图4 为在阶梯温度载荷作用下κF对半平面对称轴(x=0)处未知量uz,wz,p及θ的影响.图4(a)表明,随着时间的逐渐增加,即使在不同的κF下,沉降uz最终也会达到相同的稳定状态.κF较大时,由于在初始阶段固结速度比热传导快,导致热体积膨胀不明显,所以半平面的沉降uz变大.图4(b)表明了相对流速wz的变化趋势.当系统趋于稳定时,相对流速wz趋于0;但在初始阶段,平面内部的流速小于上表面附近,并且κF较大时流速的峰值较大.图4(c)表明孔隙压p随着时间的增加由初始孔压逐渐至0,并且半平面上表面附近孔隙压的消散速度较平面内部快,同时孔隙压的消散速度随κF增大而增大.图4(d)表明温度从上表面传播到半平面深处,并且逐渐到达至给定的表面温度,最终达到等温状态.同时,κF对变温θ的影响很小.

图5为热传导系数k对半平面对称轴(x=0)处未知量uz,wz,p及θ的影响.由图5可以看到,在给定的参数下,当k较大时,初始阶段时热传导比固结过程快,半平面的沉降uz较小,温度场达到等温状态所需的时间更短.

图5阶梯温度载荷下k对半平面时-程曲线的影响Fig.5 The effect of k on time-history curves of half-plane under the step temperature loading

图6 为附加热交换系数β对半平面对称轴(x=0)处未知量uz,wz,p及θ的影响.当β较小时,流相和固相之间的耦合作用减弱,初始阶段热体积膨胀不明显,半平面的沉降uz较大.

图6 阶梯温度载荷下β对半平面时-程曲线的影响Fig.6 The effect ofβon time-history curves of half-plane under the step temperature loading

图7为半平面对称轴不同深度处沉降uz的几何非线性结果与相应的线性解的比较.可以看到:当载荷较小时(q(x,t)=3 kN/m2),非线性解和线性解趋于一致;当载荷较大时(q(x,t)=30 kN/m2),非线性解略大于线性解.

图7 半平面的非线性解与线性解的对比Fig.7 The comparison between nonlinear solutions and linear solutions of half-plane

3.2.2 循环温度载荷下的动力学特性

限于篇幅,本工作只给出了循环温度载荷下达西渗透系数κF对系统的影响(见图8).由图8可以看出:在循环温度载荷和阶梯温度载荷作用下,半平面整体热动力学性能相似;各物理量呈周期性变化且周期相同,但在深度方向上存在相位差.

图8 循环温度载荷下κF对半平面时-程曲线的影响Fig.8 The effect ofκF on time-history curves of half-plane under the cycle temperature loading

4 结束语

在几何非线性条件和热局部平衡条件下,本工作对流体饱和不可压多孔热弹性半平面在受到表面温度载荷作用下的动力学特性进行了研究.首先基于PMT并考虑几何非线性的影响,给出了问题的数学模型;其次提出了一个系统的综合数值计算方法来模拟问题的数值结果,通过DQM和二阶后向差分格式分别在空间域和时间域离散数学模型,利用Newton-Raphson法求解非线性代数方程组,从而可得到问题的数值结果,分析半平面的动力学特性.

与现有解析解的对比性研究验证了本方法是正确可靠的,且具有精度高、计算量小、数值稳定等优点.在此基础上,分析了在表面温度载荷作用下非线性流体饱和多孔热弹性半平面的热力学特性,考察了系统的材料参数达西渗透系数κF、附加热交换系数β、热传导系数k以及几何非线性对半平面动力学特性的影响,获得了一些有益的结论.

(1)在阶梯温度载荷和循环温度载荷下,半平面整体热动力学性能相似:随着时间增加,半平面沉降趋于稳定,流速逐渐趋于0,孔隙压由初始值逐渐至零,温度逐渐上升并由上表面向纵深处传导和扩散.在循环温度载荷作用下,各物理量呈周期性变化且周期相同,但在深度方向上存在相位差.

(2)达西渗透系数较大或附加热交换系数较小时,固结作用过程快于热传导过程,初始阶半平面体沉降较大.

(3)热传导系数较大时,初始阶段时热传导比固结过程快,半平面的沉降较小,温度场达到等温状态所需的时间更短.

(4)当载荷较小时半平面的非线性解和线性解趋于一致,当载荷较大时,非线性解略大于线性解,因此当载荷较大时,有必要考虑非线性的影响.