2.3 热流固耦合理论

流固之间的耦合现象在水利水电工程、岩土工程和石油开采工程中广泛存在。在某些情况下（如核废料地质存储、地下热源开采利用等），多孔介质之间除了要考虑孔隙流体和固体之间的相互影响外，还必须考虑孔隙介质温度变化与流体和固体之间的相互影响。

2.3.1 达西方程

已知孔隙介质中流体的渗流速度φν_r，φ为孔隙度，常温下流体的动力黏滞系数为常数，达西方程可写成

实践表明，流体的黏滞系数随温度的变化而变化，且差别较大，因此对于温度变化较大的情况，需要考虑流体黏滞性对渗流速度的影响。水的动力黏滞系数与温度的关系为

不同温度情况下水的动力黏滞系数见表2.1。水的温度越低，其动力黏滞系数越大；0℃时的动力黏滞系数是100℃的动力黏滞系数的6.34倍。由此可见，水的动力黏滞系数对孔隙介质渗透系数取值影响不容忽视。

考虑动力黏滞系数后的达西定律为

考虑温度影响情况，孔隙介质水流连续性方程形式与流固耦合条件下的连续性方程形式相同，但其中的渗流速度与温度非线性相关。

式（2.71）的张量表达形式为

2.3.2 状态方程

考虑温度影响条件下，流体密度ρ_f、固体密度ρ_s、孔隙度φ和渗透率张量κ与状态变量p、T有关。

1.密度和孔隙度

基于小变形假定，孔隙介质的密度和孔隙度可表示为

式中:c_f、c_φ分别为流体和孔隙的压缩系数；K_r为固体基质体积模量;β_f、β_φ和β_Tm分别为流体、孔隙和固体基质的线性热膨胀系数。

各系数定义为

在式（2.75）中，由于假定孔隙介质为小变形，所以忽略了应力和变形对孔隙度φ的影响。其实，大多数情况下，应力或应变对孔隙度的影响应该予以考虑。

将式（2.75）对时间t求导，整理得

式中:K_b为孔隙介质整体的体积模量;β_Tb为整体线性膨胀系数;其余符号意义同前。

2.渗透率

对各向异性孔隙介质，其在各个方向上的渗透性均不相同。各向异性的结果是孔隙介质渗流在不同方向上的渗流速度也不相同。这种特性对渗流分析结果影响是很重大的。对于各向同性均质孔隙介质，渗透率在各个方向是一致的，渗透率取值大小与孔隙介质当前孔压及温度状态密切相关，其关系式为

式中:κ₀为参考状态下的孔隙介质渗透率。

式（2.77）给出了孔压、温度和固体骨架变形对渗透率影响的定量修正表达式。式（2.77）只是孔隙介质渗透率变化众多修正表达式中的一种。

2.3.3 连续方程

将达西定律描述的渗流速度式（2.70）代入连续性方程（2.71），消去速度项v_r，得

将式（2.76）代入式（2.78），整理得

式中:α为比奥系数;其余符号意义同前。

式（2.79）左边第一项描述了固体骨架孔隙体积变形，第二项则描述了流体压力变量引起的孔隙体积改变，第三项为温度变化引起的孔隙体积变化。式（2.79）右边则孔隙介质中流体体积的改变量。式（2.79）实际上描述了孔隙在骨架应力、孔隙流体压力以及变温作用下孔隙介质的孔隙体积改变量与孔隙介质体中的流量改变量相等。

2.3.4 控制方程

基于孔隙介质线弹性假设，孔隙介质总的应变是应力引起的应变、孔隙压力引起的应变和温度变化引起的应变之代数和，即

式(2.80)逆变换得

注意到比奥系数α=1-K/K_m，并应用有效应力原理，式（2.81）改写为

式中：σ_ij为总应力。

式（2.82）就是多孔连续介质热流固耦合本构方程。

根据几何方程ε_ij=(u_i,j+u_j,i)/2可得

将式（2.83）对x_j求导，并利用平衡方程σ_ij,j+f_i=0，可得力学平衡方程为

上述3个平衡方程中包含了3个位移分量以及1个孔压分量p和1个温度分量T共5个未知量。加上1个渗流微分方程，还需要补充一个方程，方能求解。

2.3.5 能量方程

根据热力学第一定律，单位体积单位时间内由外界传入系统内的能量与内部热源产生的能量之和等于物质内能的增量与力对外做功之和。考虑固体弹性变形，单相液体饱和多孔介质情形的能量方程为

式中：e、h分别为单位质量的内能和热含量（即比内能和比焓）；k_t是总的热导率。

式（2.85）左边第一项为系统内能的非稳态变化率，称累积项；第二项为热应变能的变化率；第三项是导热项；第四项为单位时间内传入与传出的能量差，是对流项；第五项是流体压力对外做的功。等式右边是总的热源强度。

式中:c_v和c_p分别为定容比热和定压比热；下标f和s分别对应液体和固体。

将式（2.85）左边第一项展开，并将第四项进行改写，整理得能量方程

2.3.6 耦合方程的求解

式（2.79）、式（2.84）和式（2.86）组合构成求解p、T、u_i的完整方程组。方程组的求解需要结合相应的边界条件和初始条件进行。

1.边界条件

热力学边界：给定域Ω界面S上已知温度边界T⁰ （x_i，t）或热通量边界），其中x_i=（x，y，z）。

渗流边界：给定域Ω界面S上已知孔压边界p⁰(x_i,t)或流量边界Q⁰(x_i,t)，其中x_i=(x,y,z)。

力学边界：给定域Ω界面S上已知应力边界σ⁰(x_i,t)或已知位移u⁰(x_i,t)，其中x_i=(x,y,z)。

2.初始条件

给定域Ω上的温度、孔压和应力的初始值T(x,0)、p(x,0)和σ(x,0)。

3.求解方法

式（2.79）、式（2.84）、式（2.86）构成的方程组具有高度的非线性特性，除极少数情况下可以通过简化获得位移、孔压和温度的解析解外，绝大多数情况都需要借助有限元等数值分析方法来求解。