渗流、温度场、应力场的耦合作用

渗流、温度场、应力场的耦合作用

1岩石体多孔介质的耦合地下水和地应力是地质环境的主要因素。它们之间的联系、相互作用和制约性是相互联系、相互制约的，因此在地下水、应力场和温度场之间形成了多阶段耦合效应。该多场耦合问题广泛存在于水利水电工程、环境工程、石油开发、地热开发、地震预报与控制、地下工程以及核废料处置等诸多领域,是人类活动和地质环境相互作用研究的重要课题。在水利水电工程和岩土工程实际应用中,考虑岩土类多孔介质的多场耦合作用,可以更加准确地评价高坝复杂基础、高地应力和高外水压力下的深埋洞室、高边坡岩体的稳定性,并为采取有效、合理的渗控措施提供科学依据。岩土类多孔介质的多场耦合问题实际上是由固体力学、流体力学、传热学、传质学、物理化学与反应理论等基础学科与众多工程学科相互交叉而形成的新兴边缘学科。研究岩土体多场耦合问题的关键在于耦合数学模型的建立,目前国内外研究者已经提出了很多数学模型,经历了由初期的研究多孔介质饱和状态下的多场耦合情形,发展到现在主要研究多孔介质在非饱和-饱和状态下的多场耦合情形,国外如Nguyen等在Biot固结理论的基础上提出了适用于饱和孔隙介质的非等温固结的数学模型;Bower等针对裂隙岩体提出了饱和双重介质的多场耦合数学模型;Rutqvist等提出了一种非饱和介质的多场耦合数学模型,该模型中较为充分地考虑了多孔介质的多相效应以及各场间的相互作用关系。国内方面,如赖远明、李宁等针对冻土问题提出了带相变的多场耦合数学模型;刘亚晨等对核废料贮存库围岩裂隙岩体介质多场耦合问题进行了探讨,等等。这些数学模型假设有些是结合具体的问题提出,这样得到的数学模型应用范围较窄;有些提出的假设虽然可以使模型大大简化,便于数值计算,但是与多场的实际耦合过程并不一定完全符合。本文在总结国内外学者对多孔介质多场耦合的机理、数学模型建立及求解方法的研究成果基础上,以岩土工程中的岩土类多孔介质为研究对象,综合运用固体力学、水力学、热学等基本理论,结合多孔介质的热本构关系以及孔隙流体的热运动规律,建立了多孔介质多相、多场全耦合数学模型方程组,在此基础上采用加权余量法和有限差分法对其分别进行了空间和时域上的离散,并讨论了数值求解方法的稳定性,最终得到了多孔介质的多相、多场耦合的有限元求解模型。2密度vl将非饱和状态下岩土工程中的岩土类多孔介质视为连续性介质,其概念模型为:假设一个非饱和多孔介质的特征单元体主要由三相构成:固相、液相和气相,其体积V可以用下式表示:其中Vs、Vl、Vg分别是固体介质、液体介质和气体介质的体积。假设ρs、ρl、ρg分别代表固体介质、液体介质和气体介质的密度,ms、ml、mg分别表示固体介质、液体介质和气体介质的质量,且满足m=ρV,ms=ρsVs,ml=ρlVl,mg=ρgVg。定义ns=Vs/V,nl=Vl/V,ng=Vg/V分别代表固相、液相、气相的体积分数;定义表征单元体的孔隙体积为Vv,则有Vv=Vl+Vg,相应的,定义n为孔隙率,有:定义单元体中液体的饱和度为S,且满足S=Vl/Vv。3多孔热耦合结构的关系模型3.1根据本构关系求解热应力应变关系当温度发生变化时,多孔介质将由于受热发生膨胀、受冷而发生收缩,从而产生应变。如果多孔介质热变形受到边界的约束,就会产生应力。这种由于温度变化而引起的应力可以称为“热应力”或“温度应力”。考虑等温条件下的岩土体本构关系比较常用,如反映非饱和土中力学行为的Camclay模型、Alonso-Gens模型和Bolzon-Schrefler模型等;反映材料塑性屈服行为的Mohr-Coulomb模型、Drucker-Prager模型等,为了利用这些成熟的本构关系,同时考虑温度变化对多孔介质的影响,可以给出以增量形式表述的热本构关系。另外,在非饱和多孔介质中,对变形起主要作用的是有效应力,于是可以将非饱和多孔介质的热应力应变关系写为:式中:σ′为有效应力张量;D为刚度张量,可以为弹性张量,可以为弹塑性张量;ε为总应变张量;εT=βT(T-T0)δ可以为受温度影响而产生的应变张量;βT为固体骨架的热膨胀系数。3.2相对渗透率系数式中:k为介质的本征饱和渗透张量;krl为液体的相对渗透性系数,非饱和区域krl∈(0,1),饱和区域krl=1.0;klT为液相的热耦合张量。3.3相热运动的广义darcy规则式中:krg为气体的相对渗透性系数;kgT为气相的热耦合张量;μg为气相的动力粘滞系数;vgr为气体相对于固体骨架的速度向量。3.4热通量矢量t式中:q为热流密度矢量即单位时间内通过单位面积的热通量矢量;T表示温度分布矢量;λ为导热系数,其影响因素包括介质的种类、材料组分、湿度、压力等。4广义darcy定律推导多孔介质非饱和状态下多场耦合的控制方程将主要基于上述非饱和多孔介质的概念模型,同时将遵循以下基本假设:(1)假设多孔介质为连续、均匀介质,多孔介质中的液体和气体只在连通孔隙中流动,气体为理想气体,且液体和气体的运动满足式(4)和式(5)所示的广义Darcy定律;(2)假设多孔介质中的固相介质的密度是应力、吸力、温度以及水压力的函数;流体介质(包括液体和气体)的密度是温度和气压的函数;(3)假设多孔介质的变形满足小变形假设,而且其应力状态满足Biot有效应力原理;无论是饱和状态,还是非饱和状态,应力应变关系满足式(3)所示的热应力应变关系;应力以拉应力为正,压应力为负;(4)假设多孔介质满足局部热平衡,即同一局部区域内固、液、气三相的温度相同;液相与气相不相混溶,且不考虑液相与气相的相互转化,热量在多孔介质中的传导过程满足式(6)所示的Fourier定律。多孔介质多场耦合的控制方程主要包括变形场控制方程、连续性方程以及能量守恒方程构成。4.1应力增量算子pg式中σ为总应力张量。对于多孔介质非饱和状态下,总应力张量的增量形式可以用Biot有效应力原理来表示,即:式中σ′为有效应力张量;δ为Kronecker张量算子;αl、αg为有效应力增量系数;pl,pg分别为孔隙水压力和孔隙气压力。将介质热应力应变关系代入上式,并将总应力表示的平衡方程替换成有效应力表示的平衡方程,于是得到增量形式的变形场控制方程有:4.2物质导数的及其应用考虑到液体的密度是压力和温度的函数1/ρldρl=clpdpl+clTdT,其中clp为液体的压缩性系数、clT为液体的热膨胀系数,二者满足clp=∂ρl/ρl∂pl,clT=-∂ρl/ρl∂pl。考虑到饱和度为温度、孔隙率以及吸力等的函数,而其中吸力和温度是影响饱和度的最重要因素,因此有:dS=cSpd(pg-pl)+cSTdT,其中cSp=∂S/∂s为吸力对饱和度的影响系数;cST=∂S/∂T为温度对饱和度的影响系数。引入物质导数,式中vs为固相介质的绝对速度向量。对于物性导数来说,在小变形条件下vs一般较小,可以认为,因此物性导数可以认为是D(×)/Dt≈∂(×)/∂t,于是。由式(10)可以得到液相介质的连续性方程:式中us为固体骨架的位移向量;εv为固体骨架的体积应变。同样地,多孔介质中的气相介质的连续性方程可以写成:式中vgr为气体平均相对速度向量。考虑到气相为理想状态下的气体,因此其状态可以用以下方程来描述:式中:mg为气体的分子量;R为常数。一般地,R/mg为常数,于是可以得到:1/ρgdρg=1/pgdpg-1/TdT,且定义cgp=1/pg,cgT=1/T分别为气体的压缩系数和温度膨胀系数,结合式(10)—式(13),可以得到多孔介质中理想气体的连续性方程:同样地,将固相介质的连续性方程写为:固相介质的密度可以写为:式中Ks为固体颗粒的压缩模量。如果将多孔介质视为弹塑性材料,近似地认为,其中K1、K2为系数,如果与多孔介质的体积变形模量K结合起来,有:,式中λ为Lame常数。于是可以得到考虑材料弹塑性的固相介质连续性方程:4.3多孔介质的热重分析式中:qs为流入单元体中固相的热通量;Cs是固相介质的比热;ρs是固相介质的密度;t是时间;σ′s为作用于单元体中固相的应力;Qs为单元体内部热源,如水泥水化热等,表示在单位时间内单位体积中发出的热量。非饱和介质单元体中液相的能量守恒方程:非饱和多孔介质单元体中气相的能量守恒方程:可以看出,σ′s、σ′l、σ′g满足下式:且有:根据广义Fourier定律,多孔介质的总热通量可以写成:式中:q、qs、ql、qg分别多孔介质总热通量、固相的热通量、液相的热通量和气相的热通量;λ为多孔介质的导热向量,是由固相、液相、气相的导热通量集合而成的,通常认为是多孔介质饱和度S的函数。联合式(18)—式(21)可以得到非饱和多孔介质的能量守恒方程:综上所述,非饱和多孔介质THM耦合的控制方程主要由式(9)、式(11)、式(12)、式(17)和式(22)构成,其中包含位移向量u、孔隙水压力pl、孔隙气压力pg、温度T和孔隙率n等七个基本未知量。4.4面力向量及边界液体流量(1)非稳定非饱和状态下的初值条件,即t=0时刻在Ω内有:式中,t、分别为面力向量和已知面力向量;ql、分别为边界液体流量和已知边界液体流量;qg、分别为边界气体流量和已知边界气体流量;qT、分别为边界热通量和已知边界热通量;n为边界面的外法向向量。5多孔介质多孔耦合模型的有限离散5.1求解方程对非饱和多孔介质水-热-力全耦合数学模型采用加权余量法进行有限元计算公式的推导:将计算域Ω离散成有限个单元(如NE个),并将基本未知量写成节点上的插值形式:,记应变矩阵B=LNu,并经过适当的整理,可以得到非饱和多孔介质水-热-力耦合问题的有限元求解方程组:其中:5.2耦合方程差分法非饱和多孔介质的三相多场耦合方程组(36)是一个以时间t为变量的常微分线性代数方程组,有限元分析中求解该方程组通常有差分法:将求解的时间域划分为若干个时间步:t0、t1、t2、…、tn、…,假设在步长Dt时间内未知量列阵X是线性变化的,将未知量列阵X按泰勒级数在ti附近展开并取线性项,并取X、F的差分格式,于是将耦合方程组(36)写成如下增量形式:式中含上标i、i+1量分别代表了时步i、i+1对应的量;α为差分系数,在0到1之间变化,α=0表示前差分格式(Euler差分格式),α=1表示后差分格式,α=0.5表示中间差分格式(Crank-Nicholson差分格式),α=0.667表示Galerkin差分格式;Dt为迭代的时间步长。5.3差分解稳定收敛的条件考虑到微分方程的稳定性与非齐次项无关上式成为:易知该方程的解析解为X=ae-ωt,其中α为任意常数;ω=K/C且总是正值,相当于上式的特征值。容易看出该解析解随时间单调趋于零。假设α为正的常数,则该解析解单调恒大于零。将式(39)写成如下差分格式为:将ΔXi+1=Xi+1-Xi代入可以得到:可以看出,为了使差分解稳定收敛,则λ必须满足两个基本要求:(1)|λ|<1,这是因为如果|λ|>1,则解是发散的;(2)λ>0,这是因为如果λ<0,则解会在正负之间来回震荡,不满足解解析解单调恒正或恒负的要求。由于ω恒正且0ue08dαue08d1,则λ的表达式中1-ωΔti+1(1-α)/(1+αωΔti+1)总是小于1,因此为了满足|λ|<1,必须有:可以看出,当αue08e0.5时,解是无条件稳定的;当0<α<0.5时,解是有条件稳定的,此时要求时间步长Δti+1满足Δti+1<2/ω(1-2α)。另一方面,还需要保证解不发生震荡即λ>0,此时时间步长Dt还需要满足Δti+1<1/(ω(1-α))。可以看出当α=1时,Δti+1<∞。因此,为了保证耦合有限元方程解的稳定收敛,将取α=1即后差分格式,此时解是无条件稳定的,时间步长Δti+1的选择从理论上讲无任何限制条件此时耦合方程组(37)的实际迭代格式为:式中6基于广义darcy定律的多相耦合非饱和岩土类多孔介质多场耦合作用的数学模型,研究的主要是某一物理场方程中因变量或源汇项受其他物理场作用的变化的数学描述,也包括本构规律的影响在控制方程中的反映,因此该类数学模型较为复杂,必定包含多个控制方程:(1)本文根据连续介质方法给出了非饱和多孔介质的概念模型,归纳总结了非饱和多孔介质的热本构关系模型,包括变形场的热本构关系,流体运动的广义Darcy定律以及温度场热传导的广义Fourier;在给出基本假设的基础上,基于质量守恒定律、能量守恒定律详细推导了非饱和多孔介质多相多场耦合的主要控制方程,包括变形场控制方程;固相、液相、气相的连续性方程;能量守恒方程等,建立了以位移、孔隙水压力、孔隙气压力、温度、孔隙率为未知量的非饱和多孔介质多相多场全耦合研究的数学模型。(2)本文给出了非饱和多孔介质多相多场耦合数学模型的相关定解条件,包括初值条件和水力、应力、温度边值条件,采用加权余量法对非饱和多孔介质多场耦合的数学模型进行了空间和时间域的离散,获得了有限元格式的多场耦合求解方程组,并详细讨论了该方程组解的稳定性,为进一步开展工程应用分析奠定了基础。在考虑温度变化的条件下,温度梯度同样是引起渗流运动的因素,为此可以将Dar