2.4 非饱和热流固耦合理论
基于线性化的湿-热弹性理论，建立热流固耦合的本构方程、水和气体的渗流微分方程以及能量方程。其中非饱和渗流部分考虑到水的蒸发-凝结相变过程，气液两相流中液相包含溶解的空气，气相中包含空气和水蒸气。
含相变过程的热流固耦合非饱和渗流在现代岩土力学中有非常重要的应用价值，特别是核废料存储库缓冲区的分析计算、冻土带路基融化的耦合过程分析中意义重大。
非饱和情况下，整体研究区域由固、液和气三相构成。固体骨架之间的孔隙空间由水和气共同占据，如果水的饱和度为s，则气相饱和度为sg=1-s。对于地下水非饱和带，气体一般是空气。
固体介质可以是孔隙岩石、裂隙岩石或土体。控制方程包括：混合物的本构方程（由多孔湿-热弹性本构关系以及修正后的有效应力公式联立导出）；水流和气体渗流微分方程；能量方程。
为建立上述控制方程，需要采用以下假定：
（1）固体介质为各向同性介质。
（2）研究区域瞬间处于热平衡状态。
（3）对地下水饱和带，孔隙压力pg为1个大气压，且其宏观速度vg=0。
（4）液体的流动遵循达西定律及其在多相流中的推广。
（5）空气遵循理想气体的状态方程，蒸汽由Kelvin关系表示。
（6）对于土壤，其湿应变大于热应变，必须予以考虑；对坚硬的岩石，湿应变忽略不计。
（7）固体速度vs与φ、s、ρ等标量梯度的乘积为高阶小量，可忽略。
（8）应力符号规定：拉正压负。
2.4.1 本构方程
基于湿-热弹性假设，孔隙介质总应变为宏观应力导致的应变、湿应变、热应变与孔压引起的应变之和。
1.湿应变与热应变
对于膨胀性的土体，当多孔介质中含水饱和度变化时，会产生明显的湿胀现象。在线弹性假设下，湿应变为
式中：βM为湿胀系数。
孔隙介质的热应变表达式为
式中：βT为热膨胀系数。
2.孔压引起的应变
对各向同性的多孔介质，孔隙压力的变化理论上只会引起固体颗粒沿3个轴线方向的应变，不会产生剪切应变。体积应变量由孔隙压力大小和固体介质的体积模量联合决定，即
，单向线性应变为
式中:
为加权平均孔隙压力。
有以下两种表达式
式中：p和pg分别为孔隙中水压力和气体压力；pc为毛细管吸力；χ为Bishop因子，在0~1之间取值。
3.本构方程
非饱和孔隙介质的变形由外应力引起的变形、孔隙压力变化引起的变形、湿度和温度变化引起的变形共同构成，总应变表达式为
对式(2.91)进行逆变换，得到应力表达式为
式中：Kb为多孔介质骨架的体积模量；βTb 和βMb分别为固体骨架的体积热膨胀系数和体积湿胀系数。
应用有效应力原理
，总应力形式的本构方程为
式中：α为比奥耦合系数。
2.4.2 平衡方程
根据固体力学平衡方程▽·σij+fi=0，若体积力只考虑重力，取铅直向上为z轴正方向，混合物密度为
式(2.94)表明多孔介质整体密度取各相分量的体积加权平均。于是，体积力表达式可写作
式中:gi为孔隙介质在各方向上的重力加速度分量。
将应力和体力表达式代入平衡方程，并对时间t求导，整理得
将几何方程εij=(ui,j+uj,i)/2代入，可得含位移分量的平衡方程
分量表达形式如下，
2.4.3 渗流方程
渗流微分方程由连续方程、达西定律、菲克定律和状态方程联立求出。
1.连续方程
非饱和多孔介质三相（水、气体、固体）连续方程分别为
式中:ρfw、ρfa、ρv、ρa和ρs分别为液相中的水和溶解空气的密度、气相中的水蒸气和空气的密度以及固体密度；
分别为水、气相中的水蒸气、气相中的空气和固体颗粒的速度，对非饱和带，
；qmw和qma分别为水和空气质量源（汇）强度，由于平衡条件下蒸发和液化的质量总相等，所以互相抵消，即qmw=0和qma=0。
上述速度均指相对于固定坐标系而言，引入液相和蒸汽相对于固体的相对速度vrf和vrv，即vrf=vf-vs，vrv=vv-vs
将上述表达式代入连续方程，展开得
将式（2.102）乘以（sρfw+ρvsg）/ρs，然后与式（2.100）相加，消去∂φ/∂t，得到孔隙介质中水的连续性方程为
孔隙介质中气体连续性方程为
式（2.103）和式（2.104）就是非饱和孔隙介质中的两个流体连续方程。
2.达西定律和菲克定律
连续性方程中包含4个速度项，为了使渗流方程中只出现p、T、s和ui 4个未知量，需要对这些速度项进行处理。对于非饱和带中的空气，认为
。液体相对速度和蒸汽相对速度分别用达西定律和菲克定律表示。固体速度项的散度可用位移分量的导数表示。据达西定律有
据菲克定律有
式中：κrw为气水两相流中水的相对渗透率；Dv为蒸汽在空气中的扩散系数；μw为温度相关的黏滞系数。
3.状态方程及相关导数
连续性方程中，不同物质的密度与孔隙压力、温度以及饱和度相关。各种密度的状态方程为
其中 R=8314J/（kmol·K）、 Rv（R/Mv）=641.5J/（kg·K）, Ra=287J/（kg·K）
式中:βw、βfa、βTm、βMm分别为液相水、溶解孔隙的体积热膨胀系数、固体基质的热膨胀系数和湿胀系数；R、Rv(R/Mv)、Ra分别为普适气体常数、水蒸气的气体常数和空气的气体常数；Mv为蒸汽的分子量，kg/kmol；RH为相对湿度；ρvs为饱和水蒸气密度。
溶解空气的压缩系数cfa可以认为与液体水的压缩系数cw相等，溶解空气的热膨胀系数βfa可以认为与水的热膨胀系数βw相等。
另外，上述变量之间还存在以下关系
式中：pc为毛细管压力。
式（2.107）中相关导数关系为
4.渗流微分方程
为建立水相和空气的渗流微分方程，将渗流本构方程和状态方程代入连续方程，并略去较小量∂pg/∂t项，整理得
这就是热流固非饱和渗流控制方程中的第二组方程。其中系数定义为
2.4.4 能量方程
非饱和渗流中需要考虑饱和度、面积力所做的功以及湿应变能的变化率。热流固非饱和渗流的能量方程可写为
式中：e、h分别为单位质量的内能和热含量（即比内能和比焓）；ht为总的热导率。
式（2.112）左边第一项为系统内能的非稳态变化率，称累积项；第二项为热应变能和湿应变能的变化率；第三项是单位时间内传入与传出的能量差，为对流项；第四项为流体压力对外做功之和；第五项是导热项。等式右边是总的热源强度。
这里有3个部分因其量值很小可以忽略，即质量对外力做的功率、动能变化率、黏性耗散率。
同时，式（2.112）中的系数
式中:cfp、cgp和csp和cp分别为定容比热和定压比热；ht是总的热导率；hf、hg、hs分别为液体、气体和固体的热导率。
与渗流微分方程类似，能量方程最终形式整理为
其中
2.4.5 耦合方程的求解
鉴于热流固非饱和渗流方程组的复杂性，即使对一些简单的问题，也不易获得其解析解。因此，对于非饱和热流固问题结合已知边界和初始条件，一般采用数值分析方法求解包含孔压、温度、饱和度以及位移矢量在内的6个未知量，进而求得研究域内的应力、流量和热量等物理量。数值求解具体过程略。