多孔介质多相流hm全耦合模型的构建及试验验证

多孔介质多相流hm全耦合模型的构建及试验验证

1thm耦合模型高废物地质处置、石油和液化天然气的地下储存、热能源的开发、深层资源的开发、地下空间的开发和利用等建设项目通常与物质、能量和变量交换有关，并表现出复杂的热耦合过程。thm。经过近30a的发展,特别是在DECOVALEX等国际及欧盟合作项目的推动下,岩土多场耦合研究取得了长足进展,开展了大量室内及现场试验[3~7],建立了一系列数学模型及数值模拟方法[8~12],开发了数十套数值模拟程序,并在代码验证及成果比较等方面取得了丰富成果[13~16]。国内尽管尚未建立大型地下实验室,但在岩石及缓冲材料的物理力学特性、地下贮库的概念设计、THM耦合的数学模型及数值模拟等方面也开展了卓有成效的研究工作。理论上,岩土THM耦合数学模型的建立以连续介质力学中的动量、质量和能量三大守恒原理为基础,其方法包括基于推广到非等温过程的Biot固结理论的宏观唯像方法[23~25]和基于混合体理论的微观均匀化方法[26~29]。尽管后者在理论上可能突破经验性本构关系的束缚,但在当前的试验技术条件下,2种方法都受到有效应力原理、土水特性和热驱动效应等经验性物理机制模型的制约。上述经验性本构模型对于建立THM耦合数学模型是必须的,但由此建立的控制方程组往往出现非对称、非协调、强非线性及高度病态等特征。特别是将气体及蒸气传输过程纳入THM耦合体系之后,控制方程组的非线性及病态特征进一步恶化,使问题的求解难度增大。在DECOVALEX等国际合作项目的研究成果[13~16]中,气体及蒸气传输过程要么完全予以忽略,要么仅将蒸气作为水的组分之一加以考虑,而尚未实现在THM耦合的框架下统一求解气体及蒸气传输。然而,对于核废料贮库等工程而言,气体及蒸气传输过程对屏障系统的安全性能及设计优化起着重要的作用。它通过水分的传输和运移影响屏障系统中的湿度分布及吸湿饱和过程,并通过吸湿膨胀效应、水力特性及热传导特性对应力–应变过程、水体传输过程和热传输过程产生深刻影响。因此,工程屏障系统中,相对湿度的时空演化是THM耦合作用的综合反映,同时也是安全性能评价的重要指标。在过去的研究中,相对湿度多通过借鉴Kelvin关系式将其表示为基质吸力和绝对温度的负指数函数,其实质是在THM耦合体系中忽略气体及蒸气传输过程的情况下引入的补充关系,本身缺乏明确的物理意义。在物理学上,相对湿度定义为蒸气压力与饱和蒸气压力的百分比。采用基质吸力(严格地讲,应该是饱和度或含水量)定义相对湿度显然混淆了岩土中水分的存在形式。基于连续介质力学原理和混合体理论,导出了多孔介质多相流THM全耦合数学模型。该模型从固、液、气三相系统的动量、质量和能量守恒出发,考虑了应力–应变、水体流动、气体传输、蒸气传输、能量传输和孔隙率演化等6个过程的耦合作用,实现了对相变、溶解、热驱动、湿度传输和吸湿膨胀等物理现象的模拟,从而力图全面反映多相流THM耦合理论发展的丰富成果,确保控制方程组的封闭性和协调性,并使相对湿度的计算和评价建立在严格的物理学意义之上。其次,通过选取位移、孔隙水压、气压、蒸气压、温度和孔隙率为基本未知量,建立了有限元数值计算格式,自主研发了多孔介质三维八自由度多相流THM全耦合有限元程序THYME3D。该程序通过克服数值求解的困难和对程序架构的合理设计,能够实现任意有物理意义的单一过程、2个过程直至上述6个过程的全耦合分析。采用法国原子能委员会(CEA)开展的膨润土THM耦合Mock-up试验验证了多相流THM全耦合数值模型及其有限元程序THYME3D的正确性,从而深化了对多相流THM全耦合控制方程组特性以及耦合过程中涉及的大量本构关系及参数的作用效应的认识。2多相流thm全耦合数学模型2.1不含水层存在的不稳定性相流质控制模型将岩土材料视为包含固相(s)、液相(l)和气相(g)三相系统以及固体(s)、水(w)、干空气(a)和水蒸气(v)四种组分的等效多孔连续介质,其在当前构形中的表征单元体积及质量分别为式中:V,M分别为三相混合体的总体积和总质量。定义Vα和Mα=ραVα分别为α相(即s,l和g)的体积和质量,ρα为α相的本征密度。α相占据的体积比例为nα=Vα/V,并定义α相的表观密度为显然有式中:ρ为三相混合体的密度,ρ=M/V。三相混合体的孔隙体积为Vv=Vl+Vg,故其孔隙率可表示为定义饱和度为:则根据式(3)~(5),可得如下关系式:为了建立多相流THM全耦合数学模型,引入如下假设:(1)液相流体和气相流体在介质的骨架空隙中运动,且每一相流体均包含水和干空气2种组分。其中,气相被视为由干空气和水蒸气2种组分构成的理想气体,二者完全混合并占据相同的空间;而液相则包含水和溶解气2种组分。(2)力学响应满足小变形条件,且假定无论岩土介质处于饱和还是非饱和状态,修正的Terzaghi有效应力原理都适用。(3)流体运动满足广义Darcy定律,气体扩散符合广义Fick定律,混合气体中水蒸气和干空气的压力满足Dalton分压定律,且饱和水蒸气压力符合Kelvin湿度定律。(4)岩体固相密度是静水压力、基质吸力、温度和平均应力的函数,而液相和气相密度则是温度和压力的函数,且后者满足理想气体状态方程。(5)热传导符合广义Fourier定律,并假定固、液、气三相处于局部热平衡状态,即在空间中的每一点,固、液、气三相具有相同的温度,而相间的热量交换仅存在于相变过程中。2.2种基本理论对于三相混合体中的任一单元,其准静态过程的动量守恒方程为式中:σ为总应力张量,规定以拉应力为正;g为重力加速度向量;∇为梯度算子。将Terzaghi有效应力原理推广到非饱和情形,并假定孔隙气和孔隙水的效应可分开考虑,则岩土介质的总应力增量可表示为式中:σ′为有效应力张量;δ为KroneckerDelta张量;pl,pg分别为孔隙水压力和孔隙气压力,规定以压力为正;αl,αg为增量有效应力系数,一般要求满足αl+αg=1。考虑介质在非等温过程中的热膨胀效应和吸湿过程中的膨胀效应,则其弹性应力–应变关系可进一步表示为式中:D为介质的弹性张量;ε为应变张量;T为温度;K为介质的排水体积模量;βT为介质的体积热膨胀系数;εsw为吸湿膨胀应变,通常可表示为基质吸力s或饱和度S的函数,其中基质吸力s定义为引入几何方程则由式(7)~(9),可得动量守恒方程为式中:u为位移向量,t为时间。需要特别指出的是,上述动量守恒方程的建立仅涉及线弹性或非线性弹性本构关系。当考虑塑性或黏性本构关系时,还应使本构关系严格满足热力学第二定律的约束。2.3混合气体状态方程引入参考于固相骨架的物质导数并考虑到式(2),则多相混合体中质量守恒方程的普遍形式可表示为式中:vs为固相骨架的运动速度,vαr≡nαvαr为α相的相对表观流速,vαr=vα-vs为α相参考于固相骨架的相对本征流速;Qα为相变产生的质量交换项。(1)固相骨架质量守恒方程根据式(6),(14),固相骨架质量守恒方程(即孔隙率演化方程)可表示为式中:ρs为固相的本征密度。考虑温度、流体压力、基质吸力及体积应力对介质骨架的作用,介质的固相密度可表示为式中:εv为体积应变;α1为Biot系数,α1≡1-K/Ks,Ks为固相颗粒的体积模量。考虑到则由式(15),(16)可得孔隙率演化方程为需要指出的是,式(17)中的近似等式(即前文提到的小变形假设)并非必须。放弃该近似等式,并对式(12)采用Lagrangian描述,即可在完全相同的框架下求解大变形多相流THM全耦合问题。(2)水体质量守恒方程根据式(14),(6),水体质量守恒方程可表示为式中:vlr为液相水体的相对表观流速;ρl为水的密度;jlg为液相与气相之间的湿度传输速率(汽化为正,液化为负),可表示为式中:psv为饱和蒸气压力,pv为水蒸气压力,ω为液相迁移系数。根据广义Darcy定律,并考虑温度梯度对液相流场的热驱动效应,则水体的相对表观流速为式中:k为介质的本征渗透率张量,krl为水体的相对渗透率,klT为水体热驱动张量,µl为水的动力黏滞系数,vlr为液相参考于固相骨架的相对本征流速。对于正压流,其质量可表示为压力和温度的函数。因此有式中:clp,clT分别为水体的压缩性系数和热膨胀系数，。饱和度S通常是基质吸力s、温度T、孔隙率n等变量的函数,但以基质吸力和温度的影响最为显著。因此有其中,将式(21)~(23)代入式(19),并考虑到式(17),有(3)气体质量守恒方程如前所述,气体是干空气(a)和水蒸气(v)二者的混合体;同时,干空气既可存在于气相混合体中,也可溶解于水中。根据式(14),(6),气相混合气体的质量守恒方程可表示为式中:vgr为混合气体的相对表观流速;ρg为混合气体的密度。设干空气和水蒸气完全混合,则有式中:ρa,ρv分别为干空气和水蒸气的密度。根据广义Darcy定律,并考虑温度梯度对气相流场的热驱动效应,则混合气体的相对表观流速为式中:krg为气体的相对渗透率,kgT为气体热驱动张量,µg为气体的动力黏滞系数,vgr为气相参考于固相骨架的相对本征流速。视混合气体为理想气体,则根据理想气体状态方程,有式中:R为普适气体常数;Mg为混合气体的分子质量,满足ρg/Mg=ρa/Ma+ρv/Ml,Ma和M1分别为干空气和水蒸气的分子质量;Rg为混合气体参数,Rg=R/Mg。由式(28)可得式中:cgp,cgT分别为混合气体的压缩性系数和热膨胀系数,cgp=1/pg,cgT=1/T。需要说明的是,气体状态方程中涉及的温度T均应理解为绝对温度。另一方面,根据式(6),(14)和Henry定律,溶解于水中的干空气的质量守恒方程可表示为式中:H为Henry溶解系数。干空气的理想状态方程为式中:pa为干空气压力;Ra为干空气的气体常数,Ra=R/Ma。根据Dalton分压定律,有因此,有式中:cap,caT分别为干空气的压缩性系数和热膨胀系数,cap=1/pa,caT=1/T。将式(23),(27),(29)代入式(25),将式(21),(23),(33)代入式(30),然后对式(25)和(30)求和,并考虑到式(17),可得考虑溶解效应的气体质量守恒方程:(4)水蒸气质量守恒方程根据式(6),(14),水蒸气的质量守恒方程可表示为式中:vvr为水蒸气的相对表观流速,可分解为水蒸气在混合气体中的扩散流速和水蒸气随混合气体一起运动的对流流速之和,即式中:vvd,vva分别为水蒸气的扩散速度和对流速度,vvd≡n(1-S)(vvr-vgr),vva≡n(1-S)vgr;vvr为蒸气参考于固相骨架的相对本征流速。其中,对流速度由式(27)确定。根据广义Fick定律,水蒸气的扩散速度与其浓度梯度成正比,即式中:Dv为水蒸气在孔隙混合气体中的分子扩散张量。同理,水蒸气的理想状态方程为式中:Rv为水蒸气的气体常数,Rv=R/Ml。于是,有式中:cvp,cvT分别为水蒸气的压缩性系数和热膨胀系数,cvp=1/pv,cvT=1/T。将式(23),(23),(37),(39)代入式(35),并考虑到式(17),有2.4元热传导张量的计算引入质量守恒方程,则多相混合体中能量守恒方程的普遍形式可表示为式中:iα为α相热传导通量,σα为α相表观应力张量,Qα为α相热源;eα为α相比内能,可表示为式中:Cα为α相比热,T0为参考温度。设比热的导数项可忽略不计,则三相混合体中固、液、气三相的能量守恒方程可分别表示为式中:L为液相水的汽化潜热;〈⋅〉为Macaulay括号函数:当x>0,〈x〉=x;否则〈x〉=0。在局部热平衡假设条件下,对式(43)~(45)求和,即可得三相混合体的能量守恒方程。其中,根据广义Fourier定律,三相混合体的热传导通量可表示为式中:λ为三相混合体的热传导张量,通常是饱和度S的函数。注意到表观应力的定义式:以及如下恒等关系式:则三相混合体的能量守恒方程可最终表示为3计算thm全耦合值的格式3.1求解目标及边界条件式(12),(18),(24),(34),(40),(48)构成了岩土介质多相流THM全耦合控制方程组,其基本未知量为位移u(t)、孔隙水压力1p(t)、孔隙气压力pg(t)、水蒸气压力pv(t)、温度T(t)和孔隙率n(t)。对于三维问题,基本未知量的个数为8个。多相流THM全耦合过程的问题提法为:设岩土介质的空间体积为Ω,边界为Γ。问题的求解目标是在空间域Ω及时间域[0,τ]内寻找一组以{u(t),p1(t),pg(t),pv(t),T(t),n(t)}为基本未知量的解答,使之在Ω内满足由式(12),(18),(24),(34),(40),(48)组成的控制方程组,同时在t=0时刻满足如下初始条件:在Γ上满足如下边界条件:式(49)~(51)中:Γu,Γσ分别为已知位移和面力边界;Γl,Γlq分别为已知液相流体压力和流量边界;Γg,Γgq分别为已知气相混合流体压力和流量边界;Γv,Γvq分别为已知水蒸气压力和流量边界;ΓT,ΓTq分别为已知温度和热通量边界;分别为边界上的已知位移和面力;分别为边界上的已知水压和水体流量;分别为边界上的已知气压和气体流量;分别为边界上的已知蒸气压和蒸气流量;分别为边界上的已知温度和热通量;u0,pl0,pg0,pv0,0T和0n分别为初始位移、初始水压力、初始气压力、初始蒸气压、初始温度和初始孔隙率。需要指出的是,式(50),(51)考虑了各个过程之间的耦合作用对边界条件的影响,且均可表示为时间的函数。此外,孔隙率不需要给定边界条件,其演化完全受微分方程及初始条件控制。3.2thm全耦合过程的有限元迭代格式上述多相流THM耦合问题往往具有很强的非线性,因而只能用数值方法求解,这里采用有限单元法。(1)空间域离散引入形函数矩阵N,并对介质位移u(t)、孔隙水压力p1(t)、气压力pg(t)、蒸气压力pv(t)、温度T(t)和孔隙率n(t)等基本未知量采用相同的插值模式,则有式中:和分别为t时刻有限元节点的位移向量、孔隙水压力向量、气压力向量、蒸气压力向量、温度向量和孔隙率向量。显然,对于u(t),N为矩阵;而对于1p(t),gp(t),pv(t),T(t)和n(t),N为行向量。采用相同的符号可精简表达,但二者的区别很容易根据具体的表达式判断。对控制偏微分方程组式(12),(18),(24),(34),(40),(48)应用Galerkin法,并通过引入边界条件式(50),(51),可得多孔介质多相流THM全耦合过程的有限元支配方程组:其中,式(54a)~(54d)中,系数矩阵C,K以及右端项向量F各元素的计算式详见附录。(2)时间域差分现对式(53)进一步应用有限差分法,即可得随时步递进的迭代计算格式。令式中:i为时间步,∆t为i时步到i+1时步的时间步长,θ为积分常数。当θ=0时,为向前差分格式;当θ=1时,为向后差分格式;当θ=0.5时,为中心差分格式;当θ=2/3时,为Galerkin差分格式。为保证数值计算的稳定性,应选取θ≥0.5。将式(55)代入式(53),得多相流THM全耦合过程的有限元迭代格式(3)求解策略多相流THM全耦合过程涉及8个基本变量,且由于耦合作用关系的强非线性和非对称性,使有限元迭代方程组式(56)呈现强非线性、高度病态、非对称、非协调和计算规模随节点数迅速增长等突出特点。特别是气体和蒸气传输过程的引入,加剧了方程组的病态及非线性,使问题的求解难度增大。为确保多相流全耦合过程得到顺利求解,以下求解策略至关重要:①真实物理力学过程的耦合作用是协调和耗散的,这种协调性和耗散性确保了多相流THM全耦合问题解的存在性。因此,在求解过程中,应尽量减少对数学模型的简化,并应利用封闭方程组进行联立统一求解。对数学模型的简化尽管可以减少基本未知量,但往往是以丧失物理过程的协调性和封闭性为代价的;另一方面,过去常采用交错法对耦合过程进行逐一迭代,这种方法由于在迭代过程中未能充分利用各过程的最新更新,因而也损失了部分协调性。协调性的丧失往往是导致数值求解收敛缓慢或求解失败的重要原因。②在求解技术上,THM全耦合过程的求解涉及自动时间步进、耦合迭代和非对称病态方程组求解三大技术。其中,耦合迭代技术涉及各种耦合变量的更新,是保证耦合关系协调性的关键;代数方程组求解技术则直接决定了求解耦合问题的能力和效率。对各项技术进行严密设计,并在代数方程组求解层面上施加严格的收敛控制,避免其出现非稳定性及这种非稳定性向耦合迭代和时步迭代层面的传播,是确保耦合过程得到成功求解的关键。③不合理的经验本构关系及其参数的引入是导致耦合过程丧失协调性的又一重要原因。对经验本构关系的物理意义进行细致分析,正确把握其局限和适用条件,有助于为耦合求解指明正确的方向。④初始条件和边界条件的协调性也是确保多相流THM全耦合过程得到顺利求解的关键因素。初始条件的确定应以试验为依据,结合特定本构关系取值,避免耦合初始参数出现矛盾现象。边界条件的确定应考虑控制方程组边界项的物理意义及各过程在边界处的耦合作用效应。3.3tyne3-4耦合分析软件的创新点基于前述多相流THM全耦合数学模型、有限元格式及求解策略,自主研发了多孔介质三维八自由度多相流THM全耦合有限元程序THYME3D。其主要特性如下:(1)THYME3D以位移、孔隙水压、气压、蒸气压、温度和孔隙率等8个变量为基本未知量,涉及应力–应变、水体流动、气体传输、蒸气传输、能量传输和孔隙率演化等6个过程。THYME3D通过程序架构的合理设计,可求解任意有物理意义的单个过程、2个过程耦合或多个过程耦合问题,这既方便了程序调试,又增强了程序功能。若以应力–应变过程的3个位移自由度为基础,而将其他耦合过程的自由度数表示为X,此即为周创兵和龚玉锋提出的3+X维耦合。THYME3D不但实现了3+X维耦合,还赋以了更强的数值模拟能力。(2)THYME3D通过对多相流THM全耦合过程协调性的深入把握和对自动时间步进、耦合迭代及非对称病态方程组求解等3大核心技术的严密设计,确保了数值计算的稳定性,特别是气体传输过程和蒸气传输过程的数值稳定性。这是该程序区别于其他THM耦合分析软件的重要特色之一。(3)THYME3D包括八节点六面体、六节点棱柱体和四节点棱锥体等单元类型,并实现了THM耦合过程中可能涉及的复杂边界条件,能够在任意时刻施加边界、移除边界或进行边界类型转换。(4)依赖于特定试验成果的经验性本构关系完全独立于程序的其他主要模块,可通过编译成动态链接库与主程序链接,从而保证了主程序代码的稳定性和使用的便利性。(5)THYME3D的架构力图全面反映多相流THM耦合理论发展的最新成果,能够实现对相变、溶解、热驱动、湿度传输和吸湿膨胀等物理现象的模拟。同时,其良好的程序架构及求解能力也为将来扩展到THMC全耦合分析奠定了坚实的基础。4动力学过程的测试与分析为了确保耦合求解成果的正确性和可靠性,首先必须确保单一物理或力学过程得到正确求解。采用内压作用下的厚壁圆筒问题、外水压力作用下圆柱体的水体传输问题、外表面已知温度下圆柱体的热传导问题等对THYME3D的应力–应变过程、流体运动过程、能量传输过程及求解器进行了广泛的测试和考题,计算结果与理论解完全吻合,表明单一物理力学过程的计算成果准确、可靠。但限于篇幅,这方面的成果这里不予介绍。4.1材料和试件加工法国原子能委员会(CEA)针对核废料处置开展的膨润土THMMock-up试验是目前较为先进、且得到广泛引用的室内试验之一,其目的是研究核废料贮库黏土屏障在高温、大温度梯度及吸湿过程中的THM耦合特性。依托欧盟合作项目“应用于贮库安全性评价的THMC耦合过程”(THERESAProject),瑞典、芬兰、英国、法国、德国、西班牙、荷兰等的相关研究团队应用该试验成果开展了大量THM耦合分析代码验证及成果比较工作。试验材料采用MX–80膨润土,经33MPa高压密实后加工成2块不同含水量的试件Cell1和Cell2。下面仅考虑Cell1的试验成果,Cell2完全类似。试件为圆柱形,直径和高度均为203mm。试件依次置于PTFE套筒、不锈钢衬及绝热封套中,以确保试验过程中体积不变,且热能损失最小,如图1所示。试件安装后施以0.5MPa的轴压。试验包括2个阶段。第一阶段为2003年5月26日15:27~9月16日09:00。在此阶段,试样底部的温度逐步从室温(20℃)升高到150℃;第二阶段为2003年9月16日14:23~2004年5月25日,此阶段在维持试样底部温度为150℃的条件下,在试样顶部施加1MPa的水压力。在试验过程中,通过试件内预埋的15个温度传感器、7个相对湿度传感器、4个孔压传感器、8个径向应力传感器和1个轴压测力计,对温度、相对湿度、孔隙水压力、径向应力和轴向应力进行了全程监控。但受试验控制条件和传感器灵敏度的影响,仅有部分成果相对可靠。4.2维力学模型根据试验条件,试件内的水、气、热传输过程为一维问题,但力学过程为轴对称问题。因此,采用1/4柱体进行研究。数值模拟的有限元网格如图2所示,其中八节点等参单元432个,节点629个。(2)初始本构关系根据试验条件,试件的初始轴压为0.5MPa,初始温度为20℃,初始含水量为13.66%,初始饱和度为0.755,初始孔隙气压力为0.101325MPa,初始相对湿度为60%,初始孔隙率为0.3242,初始孔隙水压及蒸气压则根据上述初始条件由随后介绍的相关本构关系确定。对于力学过程,模型的上、下表面及侧面均取法向支座约束。对于能量传输过程,模型的上、下表面均取已知温度边界,上表面的温度值恒为20℃,下表面的温度值在第一阶段由20℃逐步升高到150℃,升温过程见表1,在第二阶段维持150℃不变;侧面取绝热边界。对于渗流过程,模型的上表面在第一阶段取隔水边界,在第二阶段取已知水压边界,水压值为1MPa;下表面及侧面均为隔水边界。对于气体及蒸气传输过程,模型的所有边界均为气密边界。然而,需要指出的是,试验条件没有保证气密,因此气体及蒸气传输过程的边界条件是近似的。(3)水体的本构关系前述多相流THM耦合数学模型涉及大量的本构关系和计算参数。下面结合Mock–up试验成果,逐一加以介绍。但需要特别指出的是,多相流THM耦合过程涉及的本构关系是相当复杂的,有些经验性本构关系至今尚未得到学术界的一致认可。限于篇幅,本文不拟对这些本构关系做深入的评述,仅结合试验验证对部分关键本构关系的作用效应加以讨论。①MX–80膨润土的本构关系及计算参数借鉴对数弹性模型,膨润土的体积模量K(MPa)采用如下经验公式:式中:为平均有效应力;κ为柔度系数,取为0.0025;K0为初始体积模量,取为16.7MPa;Kmax为体积模量上限,取为500MPa。根据MX–80膨润土的吸湿膨胀特性,膨胀应变采用如下经验公式计算:式中:s为基质吸力(MPa);ς为试验参数,取为0.004。膨润土的土水特性曲线采用如下经验关系式:式中:ρd为干密度,取为1.791g/cm3;a=10.96;b=41.89;p0=12.68MPa;m=0.211;η=7.97;n0=0.4;α=0.00647℃-1;T0=20℃。需要指出的是,式(59)是根据含水量–基质吸力关系曲线变换得到的,未必完全满足s=0时,S=0。因此,使用时应加以修正。有效应力系数αl和αg可分别表示为式中:χ为Bishop系数。为了考虑非饱和介质的进气延滞效应,并保证饱和条件下αl=α,采用下式计算:式中:sa为非饱和介质的进气吸力,取为1kPa。膨润土的热传导系数λ(W/(m·K))与饱和度S有关,可表示为式中:A1,A2,x0,dx为参数,分别为0.57,1.28,0.65,0.10。膨润土本征渗透率k取1×10-21m2;热膨胀系数βT取为1×10-5℃-1;比热容Cs取为920J/(kg·℃)。②液相水的本构关系及计算参数借鉴液体–蒸气界面之间的相变特征,液相迁移系数ω(m·s/kg)可表示为式中:C为调节系数,取为5.5×10-5。式(63)应理解为系数C通过单位体积上的质量迁移能力。水体的密度ρl(kg/m3)采用下式计算:式中:pl0为参考水压,可取为0或一个大气压;ρlT为常压下水体密度(kg/m3),随温度的变化,其表达式为在0.01℃~360℃温度范围内,式(65)拟合相关系数为R2=0.9997。水体的热膨胀系数clT(℃-1)采用如下多项式表达:在0.01℃~260℃范围内,式(66)拟合相关系数为R2=0.9975。水的动力黏滞系数µl(Pa·s)与绝对温度T有关,即水体在介质中的相对渗透率krl可表示为饱和度S和孔隙率n的函数:式中:Ss为最大饱和度,取为1.0;Sr为残余饱和度,取为0.01;0n为初始孔隙率,取为0.3242。考虑饱和度的影响,温度梯度对水体的热驱动系数klT(m2/(s·K))可表示为饱和度S的线性函数:式中:klT0为饱和条件下的热驱动系数,取为0。③气体和蒸气的本构关系及计算参数根据Sutherland公式,气体的动力黏滞系数µg(Pa·s)为式(70)各参数取值为:c=120K,0T=291.15K,µ0=1.827×10-5Pa·s。气体在介质中的相对渗透率krg可表示为饱和度S的函数:温度梯度对气体的热驱动系数kgT(m2/(s·K))可表示为饱和度S的线性函数:式中:kgT0为完全干燥状态下的热驱动系数,取为6.0×10-12m2/(s·K)。由于干空气主要由氧和氮组成,根据二者在水中的溶解能力,Henry溶解系数H可表示为气压pg、蒸气压pv和绝对温度T的函数:气体的比热容Cg可表示为pg和pv的函数:借鉴Kelvin关系式,孔隙气中的饱和蒸气压力psv(MPa)可表达为基质吸力s(或含水量)和绝对温度T的函数:式中:s0为初始基质吸力,根据初始条件由土水特性曲线计算得到;a,b为调节系数,取为a=2.3,b=0.35。式(75)的物理意义应理解为:在特定含水量和温度条件下,孔隙气中的水蒸气在平衡态所具有的蒸气压力。由于饱和蒸气压可能受土体结构和孔隙特征等因素的影响,因此需要用调节系数a,b加以修正。水蒸气在孔隙气中的分子扩散系数Dv取为3.0×10-10m2/s。液–汽蒸发潜热L取为2.46×106J/kg。4.3成果分析(1)计算结果的分析试件中从底至顶13个温度传感器T1~T13的实测温度及计算温度变化过程曲线如图3所示(图中的数据点为不同位置的实测值,实线为不同位置的计算值)。T0和T14传感器分别位于试件的底部和顶部边界,为边界温度,予以忽略。从图中可见,计算温度过程与实测温度基本吻合。在第一阶段,计算温度主要受饱和度变化(通过热传导系数)和升温过程的影响。由于在计算过程中,边界的升温过程采用每个升温阶段始末的温度值对时间进行内插确定,而未考虑其稳定时间,因此与实测温度存在一些偏差。在第二阶段,计算温度主要受饱和度的影响。在施加水压之后的一段时间内,试件上部的饱和度迅速增大,由式(62)可知热传导系数有所增大,因此试件中、上部的计算温度过程曲线略向下弯曲。随后,随着试件从上到下逐步吸湿饱和,热传导系数在试件内的差别逐步减小,计算温度分布也逐渐趋于均匀。但在整个计算过程中,固、液、气三相的应力功、相变能和对流传输对温度分布的影响很小。此外,埋设在试件顶部197.50mm处T13传感器的计算温度与实测温度差别较大。T13距顶部边界仅5.5mm,但实测温度却明显大于边界温度,且与位于其下6.25mm处T12传感器的温度值更为接近,这可能是由于T13的初始温度率定误差或埋设位置误差造成的。(2)试件内部的孔隙率和饱和度分布试件在不同位置处的孔隙率和饱和度变化过程曲线如图4,5所示。尽管孔隙率和饱和度不是试验的直接指标,但其变化过程是THM耦合过程的综合反映。因此,分析孔隙率和饱和度的变化过程有助于理解试件中的THM耦合机制。在第一阶段,试件内的THM耦合过程受底部边界的升温过程及其内的热能传输过程驱动。随着试件从下至上逐步升温,试件中、下部出现升温膨胀、向上挤压和孔隙收缩等现象,且水分蒸发逐步加剧,并通过温度梯度对气体及蒸气的热驱动效应使水分向试件中、上部迁移、凝结。从而表现出整个试件的孔隙率有所减小;同时,试件中、下部的饱和度明显减小,而顶部的饱和度有所增大。在第二阶段,试件顶部的水压对孔隙率和饱和度的分布起关键作用。施加水压之后,试件顶部迅速趋向饱和,且吸湿饱和过程逐步向试件的中、下部推进。但由于膨润土的本征渗透率很小,吸湿饱和过程推进较为缓慢。经过250d吸湿过程之后,整个试件还远未达到完全饱和状态。与此同时,在试件顶部施加水压之后,试件顶部附近的孔隙率先因水压作用而明显增大,而后则由于试件吸湿膨胀而导致孔隙率逐渐减小;但对于试件中、下部,在顶部水压作用下,试件的孔隙率先明显减小,而后随着水体的入渗,在有效应力、水压、气压及吸湿膨胀的共同制约下,孔隙率有缓慢增大趋势。显然,在定容试验条件下,第二阶段试件中、上部孔隙率的先剧增、后渐降以及试件中、下部孔隙率的先剧降、后缓增是固相骨架质量守恒的必然结果,这从一个侧面反映了孔隙率演化方程对于维护THM耦合系统协调性所具有的重要作用。(3)湿度传输特性如前所述,孔隙介质中的相对湿度Hr常表示为基质吸力s和绝对温度T的负指数函数,即然而,式(76)混淆了孔隙介质中水分的存在形式,且与湿度传感器的物理工作机制不符。在物理学上,相对湿度定义为蒸气压pv与饱和蒸气压psv的百分比,即根据式(77),试件中7个相对湿度传感器HR1~HR7的实测和计算相对湿度变化过程曲线如图6所示,相应的计算蒸气压变化过程曲线如图7所示。由图6可知,计算给出的相对湿度过程曲线总体上较好地反映了试件中的湿度传输过程,表明给出的多相流THM耦合模型及相对湿度定义式(77)具备描述多孔介质中湿度传输特性的能力。相对湿度及蒸气压变化过程的物理机制如下:在试验的初始阶段,试件内的温度梯度较小,其对气体及蒸气的热驱动效应相对较弱,试件中的THM耦合作用以水分的缓慢蒸发为主,从而表现为相对湿度及蒸气压成整体增大趋势。在第一阶段随后的试验过程中,试件中的温度梯度逐步增大。相应地,温度梯度对气体及蒸气的热驱动效应也逐步增强,促使试件中、下部高温区的气体及蒸气向试件中、上部的低温区迁移,从而表现出试件中、下部的相对湿度及蒸气压减小,而中、上部的相对湿度及蒸气压持续增大。在第二阶段,顶部边界的水压及试件内的水体传输导致试件内的饱和度不断增大,气体及蒸气占据的空间减小,从而导致蒸气压及相对湿度增大。在此阶段,相变和热驱动将影响曲线的斜率和形态。当蒸气压达到饱和蒸气压时,液–汽相变由蒸发转变为凝结,从而使蒸气压逐步稳定到饱和蒸气压上,相对湿度也相应地稳定到100%。需要说明的是,孔隙介质中湿度传感器的监测数据可能同时受孔隙气中的蒸气浓度和介质中其他水分的影响。因此,在这个意义上,式(77)也是近似的,但其物理意义较式(76)明确,特别是对于含水量较小、基质吸力较大的干土或极干土尤为如此。对于Mock-up试验,若以式(76)定义相对湿度,则根据图5给出的饱和度变化过程曲线,不可能得到类似于图6的相对湿度变化过程曲线,也不可能对控制相对湿度变化过程的物理机制作出合理的解释。(4)轴压含量在初始和测定轴压时轴压在形态、温度和模量等基试件的实测及计算轴向压力变化过程如图8所示。但需要指出的是,根据试验报告,试验的应力成果可靠性可能较差,仅有参考意义。在第一阶段,试件所有边界处于隔水状态,尽管升温过程将影响试件内水分的分布,但试件内总的含水量不变。因此,试件的体积热膨胀是导致轴压变化的主导因素,而水分迁移导致的体积膨胀及收缩则处于次要地位。在此阶段,计算轴压逐步从初始的0.5MPa增大到1.2MPa;但实测轴压却迅速增长到5.9MPa,超过了理论分析可能的轴压范围。在第二阶段,试件内的温度分布已基本趋于稳定,而水体传输过程导致的吸湿膨胀成为轴压增大的主要驱动力。顶部边界施加水压之后,在水压、气压及饱和度变化的共同作用下,轴压迅速增大,并通过吸湿膨胀和体积模量调整,使轴压增长过程与实测值基本吻合。需要指出的是,施加水压前后实测轴压的跌落与对试验设备的操作有关。4.4土水特性分析如前所述,多相流THM耦合过程涉及大量的本构关系及计算参数。为了模拟试验过