电动法修复镉污染黏性土的试验研究

王艳; 史鑫锋; 薛逸翔; 郑佳挺; 陈冠年

doi:10.11779/CJGE2023S10004

摘要: 使用自行设计电动修复装置，针对镉污染黏性土开展了室内土柱试验，分析了土柱中镉的去除效果、修复前后污染土的微观结构变化和镉的形态特征，探究了土壤中镉的去除机制。结果表明：镉初始浓度较小时，增大电压可以显著提高镉的去除率，电压的增大可加快土中不稳定形态镉的迁移。当镉初始浓度较大时，残余在土柱中的镉含量呈现出从阳极侧到阴极侧逐渐增大的现象。经过电动修复后，土颗粒从球形呈现的黏聚体结构变为了光滑、松散的结构。电压为40V时，约95%的弱酸提取态及可还原态镉从土壤中去除，修复后土柱中的残渣态等稳定形态的镉约为总镉的60%，有效降低了镉污染土的毒性。

关键词:

电动修复 /
镉 /
污染黏性土

Abstract: A self-designed electrokinetic remediation device is used to carry out the soil column tests on cadmium- contaminated clayey soil. The removal mechanism of cadmium in soil is explored by analyzing the removal effects of cadmium in a soil column, the microstructural change of contaminated soil before and after remediation and the species analysis of cadmium. The results show that when the initial concentration of cadmium is relatively small, the removal rate of cadmium increases remarkably when increasing the voltage, since the unstable cadmium migrates faster at greater voltage. When the initial concentration of cadmium is relatively great, the residual content of cadmium exhibits a gradual increase trend from the anode to cathode. After the electrokinetic remediation, the appearance of soil particles changes from a spherical cohesive structure to a smooth and loose one. When the voltage is 40 V, about 95 % of the weak acid extractable and reducible cadmium is removed from the soil, and the residual cadmium in the soil column after remediation is about 60 % of the total cadmium, which effectively reduces the toxicity of the cadmium-contaminated soil.

Keywords:

0. 引言

可降解土体是指土体中的固相颗粒会在外部环境作用下出现质量损失、颗粒破碎等现象^[1]，进而导致土体的物理力学特性发生改变。广义上讲，自然界中的可降解土体包含城市生活垃圾、可燃冰、核废料、石灰性土壤、季节性冻土等。相比传统不可降解土体，可降解土体经历着生物、化学或物理降解，伴随生化反应、物理变化和机械运动等过程^[2]，其内部作用相互影响，往往存在物质、能量及动量交换，表现出复杂的热–水–力–化耦合过程。因此，开展可降解土体的多场耦合行为分析至关重要，对于研究与可降解土体相关的工程问题具有重要意义。

目前，国内外学者对可降解土体内部的多场耦合过程开展了试验和数值模拟研究。室内试验和现场试验的结果可为多场耦合数值模型提供参数依据；反之，数值模型可以优化并提升多场耦合试验的设计。相比传统土力学试验，开展可降解土体的多场耦合试验相对困难，且试验成本高，而进行针对可降解土体耦合行为的数值模拟则优势明显，其成本低、边界条件明确，有助于揭示可降解土体多场耦合作用的机理。现阶段，相关学者已建立了针对生活垃圾体^[2-10]、季节性冻土^[11-12]、可燃冰^[13-14]等可降解土体的多场耦合模型，加深了对可降解土体的降解过程、骨架变形、孔隙水运移、溶质迁移和孔隙气运移等过程的理解。然而，上述的多场耦合模型侧重点往往不一，对可降解土体的某些过程进行了简化^{[9, 15]}，现阶段仍有必要进一步开展可降解土体水–力–热–化等过程的耦合分析及相关模型的开发。

本文对已建立的生活垃圾多场耦合模型^{[7, 16]}进行拓展，从多孔介质质量、动量及能量守恒出发，推导了可降解土体的多场耦合模型控制方程，并给出了相关的本构关系。该模型考虑了物质相变、热量传递、应力–应变、液相流动、气体传输、组分及溶质扩散、孔隙率演化等过程。选取液相压力、气相压力、气体组分质量分数、液相溶质质量浓度、温度、骨架位移和孔隙率作为基本未知量，采用有限体积法对建立的耦合模型控制方程进行数值离散，以顺序耦合的求解思路对耦合模型进行求解。在开源场操作软件OpenFOAM中编写了相应的数值求解器，开展Liakopoulos砂柱排水试验、多孔介质中溶质迁移和热量传递、南安普顿生活垃圾CAR1试验的模拟，验证了耦合模型及数值求解器的正确性。

1. 多场耦合理论模型

骨架可降解土体多场耦合模型建立在连续介质力学的前提上。建模主要以均匀化的表征单元体（REV）为对象，引入孔隙率、饱和度、组分分数等物理量，如图 1所示，采用相应的数学方程对可降解土体的行为进行合理表征。

图 1 可降解土体表征单元体

Figure 1. Representative element volume for degradable soils

下载: 全尺寸图片幻灯片

本文提出的多场耦合模型采用的假设如下：①可降解土体固体颗粒是不可压缩的，且满足小变形假定，修正的Terzaghi有效应力原理都适用；②孔隙水和孔隙气在孔隙中连通，其运移满足达西定律，降解过程中化学溶质的质量浓度变化不影响孔隙水运移特性，不考虑溶解的气相在液相中的对流扩散运移；③假定骨架力学性质为各向同性，且不考虑温度、基质吸力对骨架变形的影响；④热传导符合广义Fourier定律，并假定固、液、气三相处于局部热平衡状态。

1.1 耦合模型基本物理力学关系

可降解土体有效孔隙率可以表示为

$n = \frac{{{V_{\text{v}}}}}{{{V_{\text{v}}} + {V_{\text{s}}}}} = \frac{{{V_{\text{l}}} + {V_{\text{g}}}}}{{{V_{\text{l}}} + {V_{\text{g}}} + {V_{\text{s}}}}} \text{，}$

(1)

式中， ${V_{\text{v}}}$ 和 ${V_{\text{s}}}$ 分别为孔隙体积和固相体积， ${V_{\text{l}}}$ 和 ${V_{\text{g}}}$ 分别为液相和气相的体积。

可降解土体液气两相的饱和度可以表示为

$\left. {\begin{array}{*{20}{l}} {{S_1} = \frac{{{V_1}}}{{{V_1} + {V_{\rm{g}}}}}}&,\\ {{S_{\rm{g}}} = \frac{{{V_{\rm{g}}}}}{{{V_1} + {V_{\rm{g}}}}}}&。{} \end{array}} \right\}$

(2)

式中， ${S_l}$ 和 ${S_{\text{g}}}$ 分别为可降解土体液相和气相饱和度。

对于可降解土体中多组分混合气体，根据道尔顿分压定律和理想气体状态方程：

$\left. {\begin{array}{*{20}{l}} {{\rho _{\rm{g}}} = \sum\limits_k {\rho _{\rm{g}}^k} ,}\\ {{p_{\rm{g}}} = \sum\limits_k {p_{\rm{g}}^k} ,}\\ {x_{\rm{g}}^k = \frac{{p_{\rm{g}}^k}}{{{p_{\rm{g}}}}} = {{\frac{{V_{\rm{g}}^k}}{{{V_{\rm{g}}}}}}^{}。}} \end{array}} \right\}$

(3)

式中： ${\rho _{\text{g}}}$ 为混合气体的密度； $\rho _{\text{g}}^k$ 为组分k的密度； ${p_{g} }$ 为混合气体的压强； $p_{\text{g}}^k$ 为组分k的分压； $x_{\text{g}}^k$ 为组分i的物质的量分数，等于其体积分数和压强分数。

各组分气体的物质的量分数和质量分数满足如下关系：

$\sum\limits_k {x_{\text{g}}^k} = \sum\limits_k {y_{\text{g}}^k} = 1 \text{，}$

(4)

式中， $y_{\text{g}}^k$ 为混合气体组分k的质量分数。

定义可降解土体的平均密度为

$\rho = {\rho _{\text{l}}}n{S_{\text{l}}} + {\rho _{\text{g}}}n{S_{\text{g}}} + {\rho _{\text{s}}}(1 - n) \text{，}$

(5)

式中， $\rho$ 为可降解土体的平均密度， ${\rho _{\text{l}}}$ 为液相的密度， ${\rho _{\text{s}}}$ 为土体骨架的表观密度。

则可降解土体的质量含水率（湿含水率）和质量含气量可以表示为

$\left. {\begin{array}{*{20}{l}} {{\theta _1} = \frac{{{\rho _1}}}{\rho }n{S_1}\quad ,}\\ {{\theta _{\rm{g}}} = \frac{{{\rho _{\rm{g}}}}}{\rho }n{S_{\rm{g}}}\quad。} \end{array}} \right\}$

(6)

对于可降解土体骨架，其运动速度 ${v_{\text{s}}}$ 可以表示为

$\nabla \cdot {v_{\text{s}}} = \frac{{\partial {\varepsilon _{\text{v}}}}}{{\partial t}} 。$

(7)

式中， ${\varepsilon _{\text{v}}}$ 为可降解土体的骨架体应变，t为时间。

1.2 多场耦合模型单元体守恒方程

准静态条件下可降解土体基质的动量守恒为

$\nabla \cdot ({\text{d}}\sigma ) + {\text{d}}(\rho g) = {\text{0}}。$

(8)

式中： $\sigma$ 为施加在可降解土体上的总应力；ρ为可降解土体的平均密度；g为重力加速度；d为增量符号。

采用Darcy定律描述流动相的动量守恒，液相和气相表观速度为

$\boldsymbol{v}_{\mathrm{lr}} \equiv\left(\boldsymbol{v}_1-\boldsymbol{v}_{\mathrm{s}}\right) \cdot\left(n S_1\right)=-\boldsymbol{M}_1 \cdot\left(\nabla p_1-\rho_1 \boldsymbol{g}\right),$

(9)

$\boldsymbol{v}_{\mathrm{gr}} \equiv\left(\boldsymbol{v}_{\mathrm{g}}-\boldsymbol{v}_{\mathrm{s}}\right) \cdot\left(n S_{\mathrm{g}}\right)=-\boldsymbol{M}_{\mathrm{g}} \cdot\left(\nabla p_{\mathrm{g}}-\rho_{\mathrm{g}} \boldsymbol{g}\right) \text { 。 }$

(10)

式中： ${v_{\text{l}}}$ 和 ${v_{\text{g}}}$ 分别为液相和气相的孔隙实际迁移速度，M_l(=Kk_lr/µ_l)和M_g(=Kk_gr/µ_g)为简化控制方程引入的液相和气相相动量，其中K为可降解土体骨架的固有渗透率或固有渗透矩阵，k_lr和k_gr分别为液相和气相在固相中的相对渗透系数，µ_l和µ_g分别为液相和气相的动力黏滞系数， ${p_{\text{l}}}$ 为液相压力。

可降解土体骨架的质量守恒可表示为

$\frac{{\partial ((1 - n){\rho _{\text{s}}})}}{{\partial t}} + \nabla \cdot ((1 - n){\rho _{\text{s}}}{v_{\text{s}}}) = {\text{B}}{{\text{Q}}_{\text{s}}} \text{，}$

(11)

式中，t为时间， ${\text{B}}{{\text{Q}}_{\text{s}}}$ 为由降解引起的源项。

将可降解土体内气体作为混合物，液相和气相的质量守恒如下：

$\frac{{\partial (n{S_{\text{l}}}{\rho _{\text{l}}})}}{{\partial t}} + \nabla \cdot (n{S_{\text{l}}}{\rho _{\text{l}}}{v_{\text{l}}}) = {\text{B}}{{\text{Q}}_{\text{l}}} + {\text{EQ}}_{\text{l}}^{\text{w}} \text{，}$

(12)

$\frac{{\partial (n{S_{\text{g}}}{\rho _{\text{g}}})}}{{\partial t}} + \nabla \cdot (n{S_{\text{g}}}{\rho _{\text{g}}}{v_{\text{g}}}) = {\text{B}}{{\text{Q}}_{\text{g}}} + {\text{D}}{{\text{Q}}_{\text{g}}} + {\text{EQ}}_{\text{g}}^{\text{w}} 。$

(13)

式中： ${\text{B}}{{\text{Q}}_{\text{l}}}$ 和 ${\text{EQ}}_{\text{l}}^{\text{w}}$ 分别表示因生化降解和蒸发造成的液相源项； ${\text{B}}{{\text{Q}}_{\text{g}}}$ 和 ${\text{D}}{{\text{Q}}_{\text{g}}}$ 分别为表征单元体内因降解和气相溶解造成的源项； ${\text{EQ}}_{\text{g}}^{\text{w}}$ 为气体中水蒸气的蒸发引起的源项。

气体成分的质量守恒为

$\begin{gathered} \frac{\partial\left(n S_{\mathrm{g}} \rho_{\mathrm{g}}^k\right)}{\partial t}+\nabla \cdot\left(n S_{\mathrm{g}} \rho_{\mathrm{g}}^k \boldsymbol{v}_{\mathrm{g}}\right)=\nabla \cdot\left(n S_{\mathrm{g}} \boldsymbol{J}_{\mathrm{g}}^k\right)+ \\ \mathrm{BQ}_{\mathrm{g}}^k+\left.\mathrm{DQ}_{\mathrm{g}}^k\right|_{k \neq w}+\left.\mathrm{EQ}_{\mathrm{g}}^k\right|_{k=w} \quad \circ \end{gathered}$

(14)

式中： $\rho _{\text{g}}^k$ 代表组分k的密度； $\boldsymbol{J}_{\mathrm{g}}^k$ 为气相中组分k的扩散通量； ${\text{BQ}}_{\text{g}}^k$ ， ${\text{DQ}}_{\text{g}}^k$ 和 ${\text{EQ}}_{\text{g}}^k$ 分别为表征单元体中因生化降解、溶解和蒸发造成的气相组分k的源项。

忽略溶解气体组分的输送，液相中溶质的质量守恒方程可以写成

$\frac{\partial\left(C_{\mathrm{v}}^s\right)}{\partial t}+\nabla \cdot\left(C_{\mathrm{v}}^s v_1\right)=\nabla \cdot \boldsymbol{J}_1^s+\mathrm{BQ}_1^s 。$

(15)

式中： $C_{v} ^s$ 为溶质s的在单位体积土体中的含量，且 $C_{\text{v}}^s$ = ${S_{\text{l}}}$ n $c_{\text{l}}^s$ ， $c_{\text{l}}^s$ 为液相中溶质质量浓度， $\boldsymbol{J}_1^s$ 为液相中溶质s的扩散通量， ${\text{BQ}}_{\text{l}}^s$ 为表征单元体中生化降解造成的溶质s的源项。

基于局部热平衡假设，可降解土体的能量守恒为

$\frac{{\partial (C(T - {T_{\rm{r}}}))}}{{\partial t}} + \nabla \cdot {\mathit{\boldsymbol{Q}}_{{\rm{conv}}}} = \nabla \cdot (\mathit{\Gamma }\nabla T) + {\rm{B}}{{\rm{Q}}_{\rm{T}}} + {\rm{E}}{{\rm{Q}}_{\rm{T}}}。$

(16)

式中：T和 ${T_{\text{r}}}$ 分别为温度和参考温度；Г和C分别为与可降解土体参考温度对应的导热系数和比热容； ${\mathit{\boldsymbol{Q}}_{{\rm{conv}}}}$ 为液相和气相流动引起的热对流； ${\text{B}}{{\text{Q}}_{\text{T}}}$ 和 ${\text{E}}{{\text{Q}}_{\text{T}}}$ 分别为由降解和水分蒸发引起的热源项。

1.3 封闭关系

液气两相之间的饱和度满足

${S_{\text{l}}} + {S_{\text{g}}} = 1 。$

(17)

可降解土体孔隙中始终存在气相和液相，将其当作存在液气两相的多孔介质，引入基质吸力概念：

$p_{\mathrm{c}}=\left\{\begin{array}{ll} p_{\mathrm{g}}-p_1 & \left(p_{\mathrm{g}}>p_1\right) \\ 0 & \left(p_{\mathrm{g}} \leqslant p_1\right) \end{array},\right.$

(18)

式中， ${p_{\text{c}}}$ 为基质吸力。

对于可降解土体固相骨架，其物质导数为^[17]

$\frac{1-n}{{\rho }_{\text{s}}}\frac{\text{D}{\rho }_{\text{s}}}{\text{D}t}=(b-n)\frac{1}{{K}_{\text{s}}}\frac{\text{D}p}{\text{D}t}-(1-b)\nabla \cdot ({v}_{\text{s}})-(b-n){\beta }_{\text{s}}\frac{\text{D}T}{\text{D}t}。$

(19)

式中：p为可降解土体骨架孔隙平均压力；β_s为固相骨架的体积热膨胀系数；K_s为固相颗粒的体积模量；b为比奥系数，b = 1–K_T/K_s，其中K_T为骨架体积模量。

同理，对于液相和气相以及气相组分，其密度的时间导数可以表示为^[17]

$\frac{1}{{{\rho _{\text{l}}}}}\frac{{{\text{D}}{\rho _{\text{l}}}}}{{{\text{D}}t}} = \frac{1}{{{K_{\text{l}}}}}\frac{{{\text{D}}{p_{\text{l}}}}}{{{\text{D}}t}} - {\beta _{\text{l}}}\frac{{{\text{D}}T}}{{{\text{D}}t}} \text{，}$

(20)

$\frac{1}{{{\rho _{\text{g}}}}}\frac{{{\text{D}}{\rho _{\text{g}}}}}{{{\text{D}}t}} = \frac{1}{{{K_{\text{g}}}}}\frac{{{\text{D}}{p_{\text{g}}}}}{{{\text{D}}t}} + \frac{1}{{{M_{\text{g}}}}}\frac{{{\text{D}}{M_{\text{g}}}}}{{{\text{D}}t}} - \frac{1}{T}\frac{{{\text{D}}T}}{{{\text{D}}t}} \text{，}$

(21)

$\begin{array}{*{20}{c}} {\frac{1}{{\rho _{\rm{g}}^k}}\frac{{{\rm{D}}\rho _{\rm{g}}^k}}{{{\rm{D}}t}} = \frac{1}{{y_{\rm{g}}^k{\rho _{\rm{g}}}}}\frac{{{\rm{D}}\left( {y_{\rm{g}}^k{\rho _{\rm{g}}}} \right)}}{{{\rm{D}}t}} = \frac{1}{{y_{\rm{g}}^k}}\frac{{{\rm{D}}y_{\rm{g}}^k}}{{{\rm{D}}t}} + \frac{1}{{{\rho _{\rm{g}}}}}\frac{{{\rm{D}}{\rho _{\rm{g}}}}}{{{\rm{D}}t}}}\\ { = \frac{1}{{y_{\rm{g}}^k}}\frac{{{\rm{D}}y_{\rm{g}}^k}}{{{\rm{D}}t}} + \frac{1}{{{p_{{\rm{g}}R}}}}\frac{{{\rm{D}}{p_{\rm{g}}}}}{{{\rm{D}}t}} + \frac{1}{{{M_{\rm{g}}}}}\frac{{{\rm{D}}{M_{\rm{g}}}}}{{{\rm{D}}t}} - \frac{1}{T}\frac{{{\rm{D}}T}}{{{\rm{D}}t}}{\rm{ }}} \end{array}$

(22)

式中：K_l为液相的压缩模量；β_l为液相的热膨胀系数； $K_{\text{g}}^{}$ 为气相的体积模量，近似等于气相绝对压力； $M_{\text{g}}^{}$ 为气相的相对分子质量； $K_{\text{g}}^k$ 为气相组分k的体积模量，近似等于气相组分k的分压； $M_{\text{g}}^k$ 为气相组分k的相对分子质量。

2. 耦合模型本构关系

2.1 可降解土体降解相变动力学

固相可降解土体的降解与土体自身特性及所处环境密切相关，可采用相关的降解动力学模型表征，其降解速率可与所处降解环境因子（如温度、含水率、酸碱度等）建立联系。

通常，可以用一阶动力学方程描述土体的降解过程，其降解速率可以表述为^[18]

${R_x} = \frac{{{\text{d}}{C_x}}}{{{\text{d}}t}} = - {f_{\text{a}}}{k_{\text{a}}}{C_x} 。$

(23)

式中：R_x为物质x的降解反应速率；C_x为土体中可降解物质x的含量； ${k_{\text{a}}}$ 为最大降解反应速率常数； ${f_{\text{a}}}$ 为降解速率环境影响因子，与含水率、温度、底物质量浓度等相关。

而对于生物降解，还可采用Monod方程描述土体的生物降解过程^[18-19]。降解速率可以表示为

${R_x} = \frac{{{\text{d}}{C_x}}}{{{\text{d}}t}} = - {f_{\text{a}}}{k_{\text{a}}}{C_{\text{B}}}\frac{{{C_x}}}{{{K_x} + {C_x}}} 。$

(24)

式中： ${C_{\text{B}}}$ 为降解微生物的含量；K_x为可降解物质x的半饱和速率常数。

针对具体的降解过程，基于上述降解动力学方程，根据反应质量守恒，可确定降解引起的 $\alpha$ 相中组分k的源项，即不同反应过程i中该物质的生成速率总和减去不同反应过程j中该物质消耗速率的总和：

${\text{BQ}}_\alpha ^k = \sum\limits_i {R_\alpha ^{kg}} - \sum\limits_j {R_\alpha ^{kd}} 。$

(25)

此外，非平衡状态下，液气交界面的水的蒸发源项可根据道尔顿蒸发定律计算，气相的溶解源项可用亨利定律表征。

2.2 可降解土体应力应变关系

定义拉应力为正、压应力为负，基于有效应力原理，则土体的总应力增量可以表示为

$\mathrm{d} \sigma=\mathrm{d} \sigma^{\prime}-b \boldsymbol{I} \mathrm{d} p \quad \text { 。 }$

(26)

式中： ${\rm{d}}{\mathit{\boldsymbol{\sigma }}^\prime }$ 和dp分别为有效应力增量和平均孔压增量；p为孔隙平均压力，可以表示为

$p = {p_{\text{l}}}{S_{\text{l}}} + {p_{\text{g}}}{S_{\text{g}}} 。$

(27)

假定土体的应变符合小应变假设，则可降解土体骨架应变可以表示为

$\mathrm{d} \boldsymbol{\varepsilon}=\frac{1}{2}\left(\nabla(\mathrm{d} \boldsymbol{u})+\nabla(\mathrm{d} \boldsymbol{u})^{\mathrm{T}}\right) \quad \text { 。 }$

(28)

式中： ${\rm{d}}\mathit{\boldsymbol{u}}$ 为骨架位移增量向量； ${\rm{d}}\mathit{\boldsymbol{\varepsilon }}$ 为骨架应变张量增量。

可降解土体骨架总应变主要包括力学压缩引起的弹性应变 ${\mathit{\boldsymbol{\varepsilon }}^e}$ 和非弹性应变 ${\mathit{\boldsymbol{\varepsilon }}^{{\rm{ne}}}}$ ，非弹性应变主要包括塑性应变 ${\mathit{\boldsymbol{\varepsilon }}^p}$ 和降解引起的应变 ${\mathit{\boldsymbol{\varepsilon }}^{\rm{b}}}$ ，可表示为

$\mathrm{d} \boldsymbol{\varepsilon}=\mathrm{d} \boldsymbol{\varepsilon}^{\mathrm{e}}+\mathrm{d} \boldsymbol{\varepsilon}^{\mathrm{ne}}=\mathrm{d} \boldsymbol{\varepsilon}^{\mathrm{e}}+\mathrm{d} \boldsymbol{\varepsilon}^{\mathrm{p}}+\mathrm{d} \boldsymbol{\varepsilon}^{\mathrm{b}} \text { 。 }$

(29)

弹性应变–应力可用拉梅常数μ和 $\lambda$ 表示为

$\mathrm{d} \boldsymbol{\sigma}^{\prime}=2 \mu\left(\mathrm{d} \boldsymbol{\varepsilon}-\mathrm{d} \boldsymbol{\varepsilon}^{\mathrm{ne}}\right)+\lambda \operatorname{tr}\left(\mathrm{d} \boldsymbol{\varepsilon}-\mathrm{d} \boldsymbol{\varepsilon}^{\mathrm{ne}}\right) 。$

(30)

对于降解，固相体积变化和孔隙体积变化的关系可表示为^[20]

${\rm{d}}{V_{\rm{v}}} = \Lambda \cdot {\rm{d}}{V_{\rm{s}}} = \mathit{\Lambda } \cdot \left( {{\rm{B}}{{\rm{Q}}_{\rm{s}}}/{\rho _{\rm{s}}}} \right) \cdot {\rm{d}}t 。$

(31)

式中： ${\text{d}}{V_{\text{v}}}$ 和 ${\text{d}}{V_{\text{s}}}$ 分别为降解造成的孔隙体积变化和固相体积变化； $\mathit{\Lambda }$ 为降解引起的孔隙体积变化与骨架体积变化之比。根据体应变表达式，可得

${\text{d}}{\varepsilon ^{\text{b}}} = \frac{{{\text{d}}{V_{\text{v}}} + {\text{d}}{V_{\text{s}}}}}{V} = \frac{{1 + \Lambda }}{{{V_0}({\varepsilon _{\text{v}}} + 1)}}({\text{B}}{{\text{Q}}_{\text{s}}}/{\rho _{\text{s}}}) \cdot {\text{d}}t 。$

(32)

2.3 可降解土体的热力参数

土体降解产热可用化学反应焓变公式计算，即标准反应焓为产物的标准生成焓值和减去反应物的标准生成焓值和：

$\Delta H_{\text{R}}^{\text{0}} = \sum {\Delta H_{{\text{fp}}}^{\text{0}}} - \sum {\Delta H_{{\text{fe}}}^{\text{0}}} 。$

(33)

式中： $\Delta H_{{\text{fp}}}^{\text{0}}$ 为生成物的标准生成焓值； $\Delta H_{{\text{fe}}}^{\text{0}}$ 为反应物的标准生成焓值。

降解热源为单位时间内各反应物的单位反应热与其反应物质量乘积的和，可以表示为

${\text{B}}{{\text{Q}}_{\text{T}}} = \sum\limits_x {\frac{{\partial {m_x}}}{{\partial t}}\Delta H_{\text{R}}^x} 。$

(34)

可降解土体的热容量用如下的公式计算^[10]：

$C = (1 - n){\rho _{\text{s}}}{H_{\text{s}}} + n{S_{\text{l}}}{\rho _{\text{l}}}{H_{\text{l}}} + n{S_{\text{g}}}{\rho _{\text{g}}}\sum\limits_k {y_{\text{g}}^kH_{\text{g}}^k} 。$

(35)

式中， ${H_{\text{s}}}$ ， ${H_{\text{l}}}$ ， $H_{\text{g}}^k$ 分别为生活土体骨架、液相和气相组分k的单位热容量。

可降解土体的有效热传导系数可根据其最大和最小值进行估算^[21]。认为可降解土体的对流传热主要由液相和气相运移引起，忽略骨架运动导致的液气运移热对流，可表示为

$\boldsymbol{Q}_{\mathrm{conv}}=\left(H_1 \rho_1 \boldsymbol{v}_{\mathrm{lr}}+\rho_{\mathrm{g}} \sum\limits_k H_{\mathrm{g}}^k y_{\mathrm{g}}^k \boldsymbol{v}_{\mathrm{gr}}\right)\left(T-T_{\mathrm{r}}\right) 。$

(36)

2.4 可降解土体的持水及扩散特性

孔隙介质基质吸力与饱和度之间的关系可以用持水曲线描述。忽略温度对饱和度的影响，常见的描述持水曲线的模型有VG模型^[11]和BC模型^[12]。土体液相和气相的相对渗透系数与液相和气相和含水率密切相关，可以采用VG-M^[22]描述。但降解过程中，土体持水特性会随之改变，后续应进一步深入研究可降解土体的持水特性。

液相中溶质主要通过分子扩散和机械弥散的形式进行主动迁移。分子扩散造成的溶质迁移与溶质质量浓度梯度成正比例关系，比例系数称为扩散系数。机械弥散是指流体通过多孔介质孔隙的平均速度和质量浓度与多孔介质断面上平均流动速度和质量浓度的不一致导致的分散现象。忽略机械弥散，液相内的溶质扩散以采用Fick定律描述：

$\boldsymbol{J}_1^s=\tau_1 D_1^s \nabla c_1^s \quad 。$

(37)

式中： ${\tau _{\text{l}}}$ 为液相内扩散弯曲因子，是介质孔隙特性参数，与溶质无关； $D_{\text{l}}^s$ 为液相内溶质s的扩散系数，在0~20℃下，液相内溶质的扩散系数可以由相关文献得到^[23]。

混合气体中各组分的扩散包括因质量浓度差异造成的扩散和因分子运动造成的扩散。其相对于混合气体平均速度的每一组分质量通量可以表示为

$\boldsymbol{J}_{\mathrm{g}}^k=\tau_{\mathrm{g}} D_{\mathrm{g}}^k \rho_{\mathrm{g}} \nabla y_{\mathrm{g}}^k \quad 。$

(38)

式中： $D_{\text{g}}^k$ 为组分k在自由混合气体中的分子扩散系数； ${\tau _{\text{g}}}$ 为气相内扩散弯曲因子。

3. 耦合模型控制方程

本文耦合模型忽略溶解的气体随液相的迁移，只考虑其界面质量传递。通过整理控制方程，将上述的部分本构关系和封闭关系带入控制方程，引入物质导数概念D(*)/Dt = ∂(*)/∂t +v_s‧∇(*)，可得到多场耦合模型的主要控制方程为

（1）液相运移方程

$\begin{gathered} \left(\frac{b-n}{K_{\mathrm{s}}} S_1^2+\frac{n S_1}{K_1}\right) \frac{\partial p_1}{\partial t}+\frac{1}{\rho_1} \nabla \cdot\left[-\rho_1 \boldsymbol{M}_1 \cdot\left(\nabla p_1-\rho_1 \boldsymbol{g}\right)\right]+ \\ \left(\frac{b-n}{K_{\mathrm{s}}} S_1 S_{\mathrm{g}}\right) \frac{\partial p_{\mathrm{g}}}{\partial t}-\left(\frac{b-n}{K_{\mathrm{s}}} S_1 p_{\mathrm{c}}-n\right) \frac{\partial S_1}{\partial t}+b S_1 \frac{\partial \varepsilon_{\mathrm{v}}}{\partial t}- \\ \quad\left[(b-n) S_1 \beta_{\mathrm{s}}+n S_1 \beta_1\right] \frac{\partial T}{\partial t}=\frac{\mathrm{BQ}_1+\mathrm{EQ}_1^{\mathrm{w}}}{\rho_1}+S_1 \frac{\mathrm{BQ}_{\mathrm{s}}}{\rho_{\mathrm{s}}} 。 \end{gathered}$

(39)

（2）气相运移方程

$\begin{array}{l} \left( {\frac{{b - n}}{{{K_{\rm{s}}}}}S_{\rm{g}}^2 + \frac{{n{S_{\rm{g}}}}}{{{p_{{\rm{g}}R}}}}} \right)\frac{{\partial {p_{\rm{g}}}}}{{\partial t}} + \frac{1}{{{\rho _{\rm{g}}}}}\nabla \cdot \left[ { - {\rho _{\rm{g}}}{{\bf{M}}_{\rm{g}}} \cdot \left( {\nabla {p_{\rm{g}}} - {\rho _{\rm{g}}}{\bf{g}}} \right)} \right] + \\ \left( {\frac{{b - n}}{{{K_{\rm{s}}}}}{S_1}{S_{\rm{g}}}} \right)\frac{{\partial {p_1}}}{{\partial t}} - \left[ {\frac{{n{S_{\rm{g}}}}}{T} + (b - n){S_{\rm{g}}}{\beta _{\rm{s}}}} \right]\frac{{\partial T}}{{\partial t}} - \\ \left( {\frac{{b - n}}{{{K_{\rm{s}}}}}{S_{\rm{g}}}{p_{\rm{c}}} + n} \right)\frac{{\partial {S_1}}}{{\partial t}} + b{S_{\rm{g}}}\frac{{\partial {\varepsilon _{\rm{v}}}}}{{\partial t}} + \frac{{n{S_{\rm{g}}}}}{{{M_{\rm{g}}}}}\frac{{\partial {M_{\rm{g}}}}}{{\partial t}}\\ = \frac{{{\rm{B}}{{\rm{Q}}_{\rm{g}}} + {\rm{D}}{{\rm{Q}}_{\rm{g}}} + {\rm{EQ}}_{\rm{g}}^w}}{{{\rho _{\rm{g}}}}} + {S_{\rm{g}}}\frac{{{\rm{B}}{{\rm{Q}}_{\rm{s}}}}}{{{\rho _{\rm{s}}}}}。 \end{array}$

(40)

（3）气相组分运移方程

$\begin{aligned} & \frac{n}{y_{\mathrm{g}}^k} \frac{\partial y_{\mathrm{g}}^k}{\partial t}+\frac{1}{\rho_{\mathrm{g}}^k} \nabla \cdot\left(y_{\mathrm{g}}^k \rho_{\mathrm{g}} M_{\mathrm{g}}\left(-\nabla p_{\mathrm{g}}+\rho_{\mathrm{g}} \mathrm{g}\right)\right)- \\ & \frac{1}{\rho_{\mathrm{g}}^k} \nabla \cdot\left(n S_{\mathrm{g}} \tau_{\mathrm{g}} \rho_{\mathrm{g}} D_{\mathrm{g}}^k \nabla y_{\mathrm{g}}^k\right)+\left[\frac{(b-n) S_{\mathrm{g}}^2}{K_{\mathrm{s}}}+\frac{n S_{\mathrm{g}}}{p_{\mathrm{g} R}}\right] \frac{\partial p_{\mathrm{g}}}{\partial t}+ \\ & \frac{(b-n) S_{\mathrm{g}} S_1}{K_{\mathrm{s}}} \frac{\partial p_{\mathrm{l}}}{\partial t}-\left[\frac{n S_{\mathrm{g}}}{T}+(b-n) S_{\mathrm{g}} \beta_{\mathrm{s}}\right] \frac{\partial T}{\partial t}- \\ & {\left[\frac{(b-n)}{K_{\mathrm{s}}} S_{\mathrm{g}} p_{\mathrm{c}}+n\right] \frac{\partial S_1}{\partial t}+b S_{\mathrm{g}} \frac{\partial \varepsilon_{\mathrm{v}}}{\partial t}+\frac{n S_{\mathrm{g}}}{M_{\mathrm{g}}} \frac{\partial M_{\mathrm{g}}}{\partial t} } \\ = & \frac{1}{\rho_{\mathrm{g}}^k}\left(\mathrm{BQ}_{\mathrm{g}}^k+\left.\mathrm{DQ}_{\mathrm{g}}^k\right|_{k \neq w}+\left.\mathrm{EQ}_{\mathrm{g}}^k\right|_{k=w}\right)+S_{\mathrm{g}} \frac{\mathrm{BQ}}{\rho_{\mathrm{s}}}。 \end{aligned}$

(41)

（4）溶质运移

$\begin{aligned} \frac{\partial\left(C_1^s\right)}{\partial t}+\nabla \cdot & {\left[-\boldsymbol{M}_1 \cdot\left(\nabla p_1-\rho_1 \boldsymbol{g}\right) /\left(S_1 n\right) \cdot C_1^s\right]+\frac{\partial \varepsilon_{\mathrm{v}}}{\partial t} C_1^s } \\ & =\nabla \cdot\left(\tau_1 D_1^s \nabla C_1^s\right)+\mathrm{BQ}_1^s。 \end{aligned}$

(42)

（5）热量传递

$\begin{gathered} \frac{\partial\left(C\left(T-T_r\right)\right)}{\partial t}+\nabla \cdot\left[\left(H_{\mathrm{l}} \rho_{\mathrm{l}} v_{l r}+\rho_{\mathrm{g}} \sum\limits_k H_{\mathrm{g}}^k y_{\mathrm{g}}^k v_{\mathrm{gr}}\right)\left(T-T_r\right)\right]- \\ \nabla \cdot(I \nabla \mathrm{T})=\mathrm{BQ}_{\mathrm{T}}+\mathrm{EQ}_{\mathrm{T}}。 \end{gathered}$

(43)

（6）骨架变形

$\begin{aligned} & \nabla \cdot {\left[\mu \nabla(\mathrm{d} \boldsymbol{u})+\mu \nabla(\mathrm{d} \boldsymbol{u})^{\mathrm{T}}+\lambda \boldsymbol{I} \operatorname{tr}(\nabla(\mathrm{d} \boldsymbol{u}))\right]-} \\ & \quad \nabla \cdot\left[2 \mu\left(\mathrm{d} \boldsymbol{\varepsilon}^{n e}\right)+\lambda \boldsymbol{I} \operatorname{tr}\left(\mathrm{d} \boldsymbol{\varepsilon}^{n e}\right)\right]=\nabla \cdot(b \boldsymbol{I} \mathrm{d} p)-\mathrm{d} \rho \boldsymbol{g} 。 \end{aligned}$

(44)

（7）孔隙率方程

$\frac{{\partial n}}{{\partial t}} = (b - n)\left( {\frac{{{S_{\text{l}}}}}{{{K_{\text{s}}}}}\frac{{\partial {p_{\text{l}}}}}{{\partial t}} + \frac{{{S_{\text{g}}}}}{{{K_{\text{s}}}}}\frac{{\partial {p_{\text{g}}}}}{{\partial t}} + \frac{{\partial {\varepsilon _{\text{v}}}}}{{\partial t}} - {\beta _{\text{s}}}\frac{{\partial T}}{{\partial t}}} \right) - \frac{{{\text{B}}{{\text{Q}}_{\text{s}}}}}{{{\rho _{\text{s}}}}} 。$

(45)

4. 耦合模型有限体积离散及程序实现

式（39）~（45）构成了可降解土体多组分耦合模型控制方程组，其基本未知量为液相压力、气相压力、气体组分质量分数、液相溶质质量浓度、温度、骨架位移增量和孔隙率。设可降解土体的空间体积为Ω，边界为Γ，则该耦合问题可以描述为：在空间域Ω及时间域[0，t]内寻找一组以[ ${p_{\text{l}}}$ ， ${p_{\text{g}}}$ ， $y_{\text{g}}^k$ ， $C_{\text{l}}^s$ ，T，du，n]为基本未知量的解，使之在空间域Ω内满足组成的控制方程组，同时在t = 0时刻满足对应的初边值条件。在t = 0时刻满足如下的初始条件：

$\phi (0) = {\phi _0}\quad ({\rm{ }}在\mathit{\Omega }内{\rm{ }}){\rm{ }}。$

(46)

在边界Γ上满足如下边界条件：

$\phi = {\phi _{\rm{b}}}\ \ \ \ {\rm{ (}}在边界{\mathit{\Gamma }_\phi }上{\rm{), }}$

(47)

$\mathit{\Phi } = {\mathit{\Phi }_{\rm{b}}}\ \ \ \ \ \ {\rm{ (}}在边界\mathit{ }{\mathit{\Gamma }_\phi 上}{\rm{ ) }}。$

(48)

式中：ϕ为基本未知量； ${\phi _{\text{0}}}$ 为基本未知量初始值； ${\phi _{\text{b}}}$ 为基本未知量在边界值；Φ为基本未知量的等效形式（流量、应力等）， ${\mathit{\Phi }_{\rm{b}}}$ 为基本未知量等效形式在边界上的值。

4.1 控制方程顺序迭代思路

考虑到建立的耦合模型涉及的场较多，采用全耦合求解方法较为昂贵，本文采用顺序耦合求解方法对耦合模型进行数值求解。对于液气两相渗流计算，模型采用的压力–压力计算形式可以保证较好的收敛性，但在进行不动点迭代求解时，易造成质量不守恒的数值失真现象，这是相饱和度的时间导数计算造成的，需要对其改进。在t+1时间步，i+1迭代步，采用如下的饱和度时间导数可以在保证较好收敛性的前提下满足质量守恒要求^[24]：

$\begin{aligned} & \left(\frac{\partial S_1}{\partial t}\right)^{t+1, i+1}=\frac{S_1^{t+1, i+1}-S_1^t}{\Delta t} \\ & =\frac{S_1^{t+1, i}+\left(\frac{\partial S_1}{\partial p_{\mathrm{c}}}\right)^{t+1, i}\left(p_{\mathrm{c}}^{t+1, i+1}-p_{\mathrm{c}}^{t+1, i}\right)-S_1^t}{\Delta t} 。 \end{aligned}$

(49)

因液气运移方程的求解先于骨架位移求解方程，为保证数值求解稳定性和收敛性，采用固定应力求解技巧，即在求解液气运移方程时，固定骨架的体应力不变^[25]。故在t+1时间步，i+1迭代步，体应变时间导数可以表示为

$\frac{{\partial {\varepsilon _{\text{v}}}}}{{\partial t}} = \frac{{\varepsilon _{\text{v}}^{t + 1,i + 1} - \varepsilon _{\text{v}}^t}}{{\Delta t}} = \frac{{\varepsilon _{\text{v}}^{t + 1,i} + \frac{b}{{{K_{{\text{dr}}}}}}({p^{t + 1,i + 1}} - {p^{t + 1,i}}) - \varepsilon _{\text{v}}^t}}{{\Delta t}} \text{，}$

(50)

式中， ${K_{{\text{dr}}}}$ 为骨架的排水体积模量，b为比奥系数。

迭代求解未知量中，将式（49），（50）代入耦合模型控制方程中，则可得耦合模型最终控制方程。在t+1时间步，i+1迭代步，方程可表示为如下形式。

对于热量传递过程：

$\frac{{\partial \left( {{A_{\rm{T}}}{T^{t + 1,i + 1}}} \right)}}{{\partial t}} + \nabla \cdot \left[ {{\mathit{\boldsymbol{B}}_{\rm{T}}}{T^{t + 1,i + 1}}} \right] - \nabla \cdot \left( {{C_{\rm{T}}}\nabla {T^{t + 1,i + 1}}} \right) = {R_{\rm{T}}},$

(51)

式中，

$\left. {\begin{array}{*{20}{l}} {{A_{\rm{T}}} = {C^{t + 1,i}},{\mathit{\boldsymbol{B}}_{\rm{T}}} = {{\left( {{H_1}{\rho _{\rm{l}}}{\mathit{\boldsymbol{v}}_{{\rm{lr}}}} + {\rho _{\rm{g}}}\sum\limits_k {H_{\rm{g}}^k} y_{\rm{g}}^k{\mathit{\boldsymbol{v}}_{{\rm{g}}r}}} \right)}^{t + 1,i}},}\\ {{C_{\rm{T}}} = {\mathit{\Gamma }^{t + 1,i}},}\\ {{R_{\rm{T}}} = {{\left( {B{Q_{\rm{T}}} + E{Q_{\rm{T}}}} \right)}^{t + 1,i}} + \frac{{\partial {{\left( {C{T_{\rm{r}}}} \right)}^{t + 1,i}}}}{{\partial t}} + }\\ {\nabla \cdot {{\left[ {\left( {{H_1}{\rho _{\rm{1}}}{\mathit{\boldsymbol{v}}_{{\rm{lr}}}} + {\rho _{\rm{g}}}\sum\limits_k {H_{\rm{g}}^k} y_{\rm{g}}^k{\mathit{\boldsymbol{v}}_{{\rm{gr}}}}} \right){T_r}} \right]}^{t + 1,i}}{\rm{ }}。} \end{array}} \right\}$

对于液相运移过程：

${A_1}{\left( {\frac{{\partial {p_1}}}{{\partial t}}} \right)^{t + 1,i + 1}} - \nabla \cdot {\left( {{\mathit{\boldsymbol{B}}_1}\nabla {p_1}} \right)^{t + 1,i + 1}} + {C_1}{\left( {\frac{{\partial {p_{\rm{g}}}}}{{\partial t}}} \right)^{t + 1,i + 1}} = {R_1},$

(52)

式中，

$\begin{aligned} & A_1=\left[\left(\frac{b-n}{K_{\mathrm{s}}} S_1^2+\frac{n S_1}{K_1}+\left(\frac{b-n}{K_{\mathrm{s}}} S_1 p_{\mathrm{c}}-n\right) \frac{\partial S_1}{\partial p_{\mathrm{c}}}\right)+\frac{\left(b S_1\right)^2}{K_{\mathrm{dr}}}\right]^{t+1, i}, \\ & \boldsymbol{B}_1=\boldsymbol{M}_1, \quad C_1=\left[\left(\frac{b-n}{K_{\mathrm{s}}} S_1 S_{\mathrm{g}}-\left(\frac{b-n}{K_{\mathrm{s}}} S_1 p_{\mathrm{c}}-n\right) \frac{\partial S_1}{\partial p_{\mathrm{c}}}\right)+\frac{b^2 S_1 S_{\mathrm{g}}}{K_{\mathrm{dr}}}\right]^{t+1, i}, \\ & R_1=\left[\frac{\mathrm{BQ}_1+\mathrm{EQ}_1^w}{\rho_1}+S_1 \frac{\mathrm{BQ}}{Q_{\mathrm{s}}}\right]^{t+1, i+1}-\left[\nabla \cdot\left(\rho_{\mathrm{l}} \boldsymbol{M}_1 \boldsymbol{g}\right)\right]^{t+1, i}- \\ & \left(b S_1\right)^{t+1, i} \frac{\varepsilon_{\mathrm{v}}^{t+1, i}-\varepsilon_{\mathrm{v}}^t}{\Delta t}+\left[(b-n) S_1 \beta_{\mathrm{s}}+n S_1 \beta_1\right]^{t+1, i} \frac{T^{t+1, i+1}-T^t}{\Delta t}- \\ & {\left[-\left(\frac{b-n}{K_{\mathrm{s}}} S_1 p_{\mathrm{c}}-n\right) \frac{\partial S_1}{\partial p_{\mathrm{c}}}-\frac{\left(b S_1\right)^2}{K_{\mathrm{dr}}}\right]^{t+1, i} \frac{p_1^{t+1, i}-p_1^t}{\Delta t}-} \\ & \left.\left[\left(\frac{b-n}{K_{\mathrm{s}}} S_1 p_{\mathrm{c}}-n\right) \frac{\partial S_1}{\partial p_{\mathrm{c}}}\right)-\frac{b^2 S_1 S_{\mathrm{g}}}{K_{\mathrm{dr}}}\right]^{t+1, i} \frac{p_{\mathrm{g}}^{t+1, i}-p_{\mathrm{g}}^t}{\Delta t}- \\ & {\left[-\left(\frac{b-n}{K_{\mathrm{s}}} S_1 p_{\mathrm{c}}-n\right)\right]^{t+1, i} \frac{S_1^{t+1, i}-S_1^t}{\Delta t} 。} \end{aligned}$

对于气相运移过程：

$A_{\mathrm{g}}\left(\frac{\partial p_{\mathrm{g}}}{\partial t}\right)^{t+1, i+1}-\nabla \cdot\left(\boldsymbol{B}_{\mathrm{g}} \nabla p_{\mathrm{g}}\right)^{t+1, i+1}+C_{\mathrm{g}}\left(\frac{\partial p_1}{\partial t}\right)^{t+1, i+1}=R_{\mathrm{g}},$

(53)

式中，

$\begin{aligned} & A_{\mathrm{g}}=\left[\frac{b-n}{K_{\mathrm{s}}} S_{\mathrm{g}}^2+\frac{n S_{\mathrm{g}}}{p_{\mathrm{g}}}-\left(\frac{b-n}{K_{\mathrm{s}}} S_{\mathrm{g}} p_{\mathrm{c}}+n\right) \frac{\partial S_{\mathrm{w}}}{\partial p_{\mathrm{c}}}+\frac{\left(b S_{\mathrm{g}}\right)^2}{K_{\mathrm{dr}}}\right]^{t+1, i} \\ & \boldsymbol{B}_{\mathrm{g}}=\boldsymbol{M}_{\mathrm{g}}, C_{\mathrm{g}}=\left[\frac{b-n}{K_{\mathrm{s}}} S_{\mathrm{w}} S_{\mathrm{g}}+\left(\frac{b-n}{K_{\mathrm{s}}} S_{\mathrm{g}} p_{\mathrm{c}}+n\right) \frac{\partial S_{\mathrm{w}}}{\partial p_{\mathrm{c}}}+\frac{b^2 S_{\mathrm{w}} S_{\mathrm{g}}}{K_{\mathrm{dr}}}\right]^{t+1, i} \\ & R_{\mathrm{g}}=\left[\frac{\mathrm{BQ}_{\mathrm{g}}+\mathrm{DQ}_{\mathrm{g}}+\mathrm{EQ}_{\mathrm{g}}^w}{S_{\mathrm{g}} \rho_{\mathrm{g}}}+\frac{\mathrm{BQ}}{\rho_{\mathrm{s}}}\right]^{t+1, i+1}-\left[\rho_{\mathrm{g}} \nabla \cdot\left(M_{\mathrm{g}} g\right)\right]^{t+1, i}- \\ & \left(b S_{\mathrm{g}}\right)^{t+1, i} \frac{\varepsilon_{\mathrm{v}}^{t+1, i}-\varepsilon_{\mathrm{v}}^t}{\Delta t}+\left(\frac{n S_{\mathrm{g}}}{T}+(b-n) S_{\mathrm{g}} \beta\right)^{t+1, i} \frac{T^{t+1, i+1}-T^t}{\Delta t}- \\ & \left(\frac{n S_{\mathrm{g}}}{M_{\mathrm{g}}}\right)^{t+1, i} \frac{M_{\mathrm{g}}^{t+1, i}-M_{\mathrm{g}}^t}{\Delta t}-\left[-\left(\frac{b-n}{K_{\mathrm{s}}} S_{\mathrm{g}} p_{\mathrm{c}}+n\right) \frac{\partial S_{\mathrm{w}}}{\partial p_{\mathrm{c}}}-\right. \\ & \left.\frac{b^2 S_{\mathrm{w}} S_{\mathrm{g}}}{K_{\mathrm{dt}}}\right]^{t+1, i} \frac{p_{\mathrm{w}}^{t+1, i}-p_{\mathrm{w}}^t}{\Delta t}-\left[\left(\frac{b-n}{K_{\mathrm{s}}} S_{\mathrm{g}} p_{\mathrm{c}}+n\right) \frac{\partial S_{\mathrm{w}}}{\partial p_{\mathrm{c}}}-\frac{\left(b S_{\mathrm{g}}\right)^2}{K_{\mathrm{dr}}}\right]^{t+1, i} \\ & \left.\frac{p_{\mathrm{g}}^{t+1, i}-p_{\mathrm{g}}^t}{\Delta t}-\left[\frac{b-n}{K_{\mathrm{s}}} S_{\mathrm{g}} p_{\mathrm{c}}+n\right)\right]^{t+1, i} \frac{S_{\mathrm{w}}^{t+1, i}-S_{\mathrm{w}}^t}{\Delta t} \cdot \end{aligned}$

对于气相组分运移过程：

$\begin{array}{*{20}{c}} {{A_k}{{\left( {\frac{{\partial y_{\rm{g}}^k}}{{\partial t}}} \right)}^{t + 1,i + 1}} + \nabla \cdot {{\left( {{\mathit{\boldsymbol{B}}_k}y_{\rm{g}}^k} \right)}^{t + 1,i + 1}} - \nabla \cdot {{\left( {{C_k}\nabla y_{\rm{g}}^k} \right)}^{t + 1,i + 1}} + }\\ {{D_k}{{\left( {y_g^k} \right)}^{t + 1,i + 1}} = {R_k},} \end{array}$

(54)

式中，

$\begin{aligned} & A_k=(n)^{t+1, i}, B_k=\left(v_{\mathrm{gr}}\right)^{t+1, i+1}, C_k=\left(n S_{\mathrm{g}} \tau_{\mathrm{g}} D_{\mathrm{g}}^k\right)^{t+1, i+1}, D_k= \\ & \left\{\left[\frac{(b-n) S_{\mathrm{g}}^2}{K_{\mathrm{s}}}+\frac{n S_{\mathrm{g}}}{p_{\mathrm{g}}}\left[\frac{\partial p_{\mathrm{g}}}{\partial t}+\frac{(b-n) S_{\mathrm{g}} S_1}{K_{\mathrm{s}}} \frac{\partial p_1}{\partial t}-\left[\frac{n S_{\mathrm{g}}}{T}+\right.\right.\right.\right. \\ & \left.\left.(b-n) S_{\mathrm{g}} \beta_{\mathrm{s}}\right] \frac{\partial T}{\partial t}-\left[\frac{(b-n)}{K_{\mathrm{s}}} S_{\mathrm{g}} p_{\mathrm{c}}+n\right] \frac{\partial S_l}{\partial t}\right\}^{t+1, i+1}+ \\ & \left(b S_{\mathrm{g}}\right)^{t+1, i+1} \frac{\varepsilon_{\mathrm{v}}^{t+1, i}-\varepsilon_{\mathrm{v}}^t}{\Delta t}+\left(\frac{b^2 S_{\mathrm{g}}}{K_{\mathrm{dr}}}\right)^{t+1, i+1} \frac{(p)^{t+1, i+1}-(p)^{t+1, i}}{\Delta t}+ \\ & \left(\frac{n S_{\mathrm{g}}}{M_{\mathrm{g}}}\right)^{t+1, i+1} \frac{M_{\mathrm{g}}^{t+1, i}-M_{\mathrm{g}}^t}{\Delta t}-\left(S_{\mathrm{g}} \frac{\mathrm{BQ}_{\mathrm{s}}}{\rho_{\mathrm{s}}}\right)^{t+1, i+1}, R_k=\left[\left(\mathrm{BQ}_{\mathrm{g}}^k+\right.\right. \\ & \left.\left.\left.\mathrm{DQ}_{\mathrm{g}}^k\right|_{k \neq w}+\left.\mathrm{EQ}_{\mathrm{g}}^k\right|_{k=w}\right) / \rho_{\mathrm{g}}\right]^{t+1, i+1} . \end{aligned}$

对于液相溶质运移过程：

$\begin{gathered} \left(\frac{\partial C_1^{\mathrm{s}}}{\partial t}\right)^{t+1, i+1}+\nabla \cdot\left(\boldsymbol{B}_{\mathrm{s}} C_1^{\mathrm{s}}\right)^{t+1, i+1}-\nabla \cdot\left(C_{\mathrm{s}} \nabla C_1^{\mathrm{s}}\right)^{t+1, i+1}+ \\ D_{\mathrm{s}}\left(C_1^{\mathrm{s}}\right)^{t+1, i+1}=R_{\mathrm{s}}, \end{gathered}$

(55)

式中，

$\begin{aligned} & \boldsymbol{B}_{\mathrm{s}}=\left(v_{\mathrm{lr}} / n S_1\right)^{t+1, i+1}, C_{\mathrm{s}}=\left(\tau_{\mathrm{l}} D_1^s\right)^{t+1, i+1}, \quad D_{\mathrm{s}}=\frac{b S_1^{+t i, t}}{K_{\mathrm{dr}}} \\ & \frac{p_1^{t+1, i+1}-p_1^t}{\Delta t}-\frac{b S_1^{t+1, i}}{K_{\mathrm{dr}}} \frac{p_1^{t+1, i}-p_1^t}{\Delta t}+\frac{b S_{\mathrm{g}}^{t+1, i}}{K_{\mathrm{dr}}} \frac{p_{\mathrm{g}}^{t+1, i+1}-p_{\mathrm{g}}^t}{\Delta t}- \\ & \frac{b S_{\mathrm{g}}^{t+1, i}}{K_{\mathrm{d}}} \frac{p_{\mathrm{g}}^{t+1, i}-p_{\mathrm{g}}^t}{\Delta t}+\frac{\varepsilon_{\mathrm{v}}^{t+1, i}-\varepsilon_{\mathrm{v}}^t}{\Delta t}, R_{\mathrm{s}}=\left(\mathrm{BQ}_1^{\mathrm{s}}\right)^{t+1, i+1} \mathrm{。} \end{aligned}$

对于骨架变形过程：

$\nabla \cdot \left[ {{A_u}\nabla \left( {{\rm{d}}{\mathit{\boldsymbol{u}}^{t + 1,i + 1}}} \right)} \right] = {R_{\rm{u}}},$

(56)

式中，

$\begin{aligned} & A_{\mathrm{u}}=2 \mu+\lambda, R_{\mathrm{u}}=-\nabla \cdot\left[-(\mu+\lambda) \nabla\left(\mathrm{d} \boldsymbol{u}^{t+1, i}\right)+\mu \nabla\left(\mathrm{d} \boldsymbol{u}^{t+1, i}\right)^T+\right. \\ & \left.\lambda \boldsymbol{I} \operatorname{tr}\left(\nabla\left(\mathrm{d} \boldsymbol{u}^{t+1, i}\right)\right)\right], \quad \nabla \cdot\left[\left(2 \mu \mathrm{d} \boldsymbol{\varepsilon}^{\mathrm{n}}+\lambda \boldsymbol{I} \operatorname{tr}\left(\mathrm{d} \boldsymbol{\varepsilon}^{\mathrm{we}}\right)\right)^{t+1, i}\right]+\nabla \cdot \\ & \left(b \boldsymbol{I} p^{t+1, i+1}\right)-\left(\mathrm{d} \rho^{t+1, i+1} \boldsymbol{g}\right)_。 \end{aligned}$

对于孔隙率更新过程：

${\text{ }}{\left( {\frac{{\partial n}}{{\partial t}}} \right)^{t + 1,i + 1}} + {A_n}{n^{t + 1,i + 1}} = {R_n} \text{，}$

(57)

式中，

$\begin{array}{l}{A}_{n}={\left(\frac{{S}_{l}}{{K}_{\text{s}}}\frac{\partial {p}_{l}}{\partial t}+\frac{{S}_{\text{g}}}{{K}_{\text{s}}}\frac{\partial {p}_{\text{g}}}{\partial t}+\frac{\partial {\epsilon }_{\text{v}}}{\partial t}-{\beta }_{\text{s}}\frac{\partial T}{\partial t}\right)}^{t+1,i+1}，{R}_{n}=\\ {\left(b{A}_{n}-\frac{{\text{BQ}}_{\text{s}}}{{\rho }_{\text{s}}}\right)}^{t+1,i+1}。\end{array}$

上述方程中，带下标的A，B，C，D，R分别为相应控制方程的系数，包括非线性项和耦合项，其下标代表待求解场量，在迭代求解过程中均可基于当前场量对其显式更新。

4.2 有限体积离散

采用单元中心有限体积法^[26]进行耦合模型数值离散求解，需要对上述问题进行求解域离散、方程离散和边界条件离散。

（1）求解域离散

求解域的离散包括空间域离散和时间域离散。时间域离散将时间划分为一系列时间步长Δt，然后采用一阶向后差分方法进行方程的求解。空间域的离散将待求解问题空间区域离散为一系列有限控制体。所示为一个典型的控制单元 ${V_P}$ ，其中心为P，包括不同的面f和面的面积向量S_f（大小为面的面积，指向为面的外法向向量方向）。其邻近控制单元V_N和其共有面f，两个网格中心的指向向量为d_f（由P指向N），向量与面f的交点为f_i。

图 2 有限体积法单元体示意图^[26]

Figure 2. Illustration of element for finite volume method^[26]

下载: 全尺寸图片幻灯片

在进行有限体积离散的时候，每个偏微分方程的项首先需要在有限控制体 ${V_P}$ 内进行积分。大部分空间导数项都使用广义高斯定律来将体积分转化为面积分，并对其近似线性化^[27]：

${\displaystyle \underset{{V}_{P}}{\int }\nabla \oplus \varphi \text{d}V}={\displaystyle \underset{\partial {V}_{P}}{\oint }\text{d}S\oplus \varphi }={\displaystyle \sum _{f}{\displaystyle \underset{f}{\oint }\text{d}S\oplus {(\varphi )}_{f}}}\approx {\displaystyle \sum _{f}{S}_{f}\oplus {(\varphi )}_{f}}。$

(58)

式中：𝜙为任意张量场；⊕为任意张量积，包含内积、外积、叉乘等，相对应的就是散度、梯度、旋度等。在积分之后，需要使用合适的格式来进行线性化。

（2）控制方程离散

对模型控制方程进行有限体积离散，在每一时间步，在控制体 ${V_P}$ 上对方程（51）~（57）进行体积分，可获得针对各未知场量的单元体积分方程。如热量传递方程可以表述为

$\int_{V_p} \frac{\partial\left(A_T T\right)}{\partial t} \mathrm{~d} V+\oint_{\partial V_p}\left(\mathrm{~d} \boldsymbol{S} \cdot \boldsymbol{B}_T T\right)-\oint_{\partial V_p}\left(\mathrm{~d} \boldsymbol{S} \cdot C_T \nabla T\right)=\int_{V_p} R_{\mathrm{T}} \mathrm{d} V_{\circ}$

(59)

液相和气相运移、组分和溶质运移、骨架变形和孔隙率更新对应的积分方程可类似整理。采用有限体积法对上述积分方程中的时间导数项、散度项、拉普拉斯项、源项进行进一步离散，即可得到关于待求变量在控制体上的线性方程。其中的方程系数及方程源项可通过上个时间步和上个迭代步的结果进行显式计算求得。

（3）边界条件离散

模型求解采用压力、位移等变量作为待求未知量，无法直接设置流量边界、应力边界，需对此类边界条件进行改进。

对于流量边界，可采用显式隐式迭代计算方法，将其转变为压力的纽曼边界条件^[28]：

$\underbrace{\nabla p_\alpha \cdot \boldsymbol{n}}_{\text {隐式 }}=\underbrace{\frac{\boldsymbol{M}_\alpha \rho_\alpha \boldsymbol{g} \cdot \boldsymbol{n}}{M_{\alpha, n}}-\frac{q}{M_{\alpha, n}}}_{\text {显式 }} 。$

(60)

式中：q为固定流量值（单位时间内的入渗速率，单位为m/s）；n为边界面外法向向量；M_{α, n}为相动量M_α的法向分量。

对于应力边界条件，采用显式隐式迭代计算，可将其转变为位移增量的纽曼边界^[29]：

$\begin{aligned} & \underbrace{(2 \mu+\lambda) \nabla(\mathrm{d} \boldsymbol{u}) \cdot \boldsymbol{n}}_{\text {隐式 }} \\ & =\mathrm{d} \boldsymbol{F}-\underbrace{(\mu \nabla(\mathrm{d} \boldsymbol{u}) \cdot \boldsymbol{n}+\lambda \operatorname{tr}(\nabla(\mathrm{d} \boldsymbol{u})) \boldsymbol{n}-(\mu+\lambda) \nabla(\mathrm{d} \boldsymbol{u}) \cdot \boldsymbol{n})}_{\text {显式耦合 }}+ \\ & \underbrace{\left(2 \mu \mathrm{d} \boldsymbol{\varepsilon}^{\mathrm{ne}} \cdot \boldsymbol{n}+\lambda \operatorname{tr}\left(\mathrm{d} \boldsymbol{\varepsilon}^{\mathrm{ne}}\right) \boldsymbol{n}\right)+b \mathrm{~d} p \boldsymbol{d} \rho \boldsymbol{g}}_{\text {显式非线性耦合 +显式压力耦合 }} \text { 。 } \\ & \end{aligned}$

(61)

式中， ${\rm{d}}\mathit{\boldsymbol{F}}$ 为边界应力增量。

4.3 求解步骤

（1）单一场量求解步骤

引入边界条件，每一控制方程最终被离散为单个控制体中心M上的线性方程：

$\alpha _{\text{M}}^{(F)}{(F)_{\text{M}}} + \sum\limits_N {\alpha _N^{(F)}{{(F)}_N}} = r_{\text{M}}^{(F)} \text{，}$

(62)

式中， $\alpha _{\text{M}}^{}$ ， $\alpha _N^{}$ 和 $r_{\text{M}}^{}$ 分别为未知场量F在空间和时间离散后对角系数、相邻系数和源项。在空间所有控制体上对各未知场量方程进行有限体积离散，则可组合得到针对该未知场量F的线性方程组：

$\boldsymbol{A F}=\boldsymbol{b}。$

(63)

式中：A为对称稀疏矩阵；F为未知场量F的列向量；b为离散后方程组右边列向量。

（2）模型总体求解步骤

耦合模型涉及较多未知变量，数值方程具有强非线性、高度病态、非对称等特点，故采用顺序迭代方法对其进行数值求解。为减小采用该方法造成的数值计算误差，引入CFL条件，在计算过程中限制最大时间步长，并尽可能采用较小的时间步长。此外，顺序耦合迭代求解使该模型易于扩展，其他复杂或非线性的本构关系可以在不改变全局求解程序的前提下嵌入模型中。耦合模型的总体求解步骤如图 3所示。

图 3 耦合模型顺序迭代耦合求解示意图

Figure 3. Schematic diagram of sequential iterative coupling solution of coupled model

下载: 全尺寸图片幻灯片

（3）收敛准则

为判别迭代计算是否达到收敛，模型中引入了绝对收敛准则和相对收敛准则。

$\left.\begin{array}{l} \max \left|\phi^{t+1, i+1}-\phi^{t+1, i}\right| \leqslant \varepsilon_\phi, \\ \max \left|\frac{\phi^{t+1, i+1}-\phi^{t+1, i}}{\phi^{t+1, i+1}}\right| \leqslant \zeta_\phi \circ \end{array}\right\}$

(64)

式中， ${\varepsilon _\phi }$ 和 ${\zeta _\phi }$ 分别为变量ϕ的收敛判别绝对精度和相对精度。

4.4 数值计算程序

依托于OpenFOAM^[30]编写相应的数值求解程序，将上述的离散过程在foam-extend-4.0版本中编程实现。程序核心代码分为3个部分，即求解器部分、边界条件部分、材料模型部分，如图 4所示。程序具有良好的拓展性，可在此基础上自定义材料本构模型。

图 4 耦合模型数值程序框架

Figure 4. Framework of numerical code for coupled model

下载: 全尺寸图片幻灯片

5. 模型及程序验证

5.1 Liakopoulos砂柱排水固结试验

采用上述模型及程序模拟经典Liakopoulos砂柱排水试验^[31]，该试验常被用于变形多孔介质中三相流动问题验证的基准案例。该砂柱由Del Monte砂组成，高度为1.0 m，直径为0.1 m。首先，从砂柱顶端注入稳定的水流，直到砂柱底部完全排水，确保整个砂柱处于完全饱和状态。试验开始时，停止注水，在重力作用下，孔隙水从底部排除，砂柱同时产生固结沉降变形。Del Monte砂的持水曲线及液相相对渗透率采用如下的关系式描述：

$\left.\begin{array}{l} S_{\mathrm{w}, \mathrm{r}}=1-1.9722 \times 10^{-11} p_{\mathrm{c}}^{2.4279} \\ k_{\mathrm{rw}}=1-2.207\left(1-S_{\mathrm{w}, \mathrm{r}}\right)^{1.0121} \end{array},\right\}$

(65)

该试验中没有测定砂土的力学参数和气相相对渗透系数，后人研究均假设砂柱为弹性介质，并假定气相的相对渗透系数用Brooks Corcy模型描述：

${k_{{\text{rg}}}} = {(1 - {S_{{\text{w,eff}}}})^2}(1 - S_{{\text{w,eff}}}^{5/3}) 。$

(66)

实际条件下，当砂柱处于完全饱和状态下时，气相可以以气泡的形式迁移，为考虑这一现象，取砂柱的气相相对渗透系数下限为1.0×10^{-4 [32-33]}。数值模型底部边界为自由排水、自由排气、位移完全约束边界： ${q_{\text{l}}}$ =0， ${p_{\text{g}}}$ = ${p_{{\text{atm}}}}$ ；顶部边界设为隔水、透气、自由位移边界： ${p_{\text{l}}}$ =0， ${p_{\text{g}}}$ = ${p_{{\text{atm}}}}$ ，du=0；四周边界设为隔水、不透气、法向位移约束边界： ${q_{\text{l}}}$ =0， ${q_{\text{g}}}$ =0，du_x=0，du_y=0。初始状态下，液气压力均设为0以表征砂柱处于完全饱和状态。数值模拟中材料参数和模型参数均取自前人研究结果^[32]。

图 5所示为砂柱不同时刻竖向位移、气相压力、液相压力随高度以及流速流量变化曲线图。可以发现，重力作用下砂柱逐渐排水，砂柱从饱和状态逐渐变为非饱和状态，同时伴随着固结沉降的发生。本文模型模拟的液相压力结果与试验结果存在差异，但与已有学者模拟的液相压力结果吻合。本文模拟结果表明在120 min后，砂柱顶部位移为1.5 mm，与已有学者模拟结果基本吻合。砂柱初始时刻的排水速率实测值为0.0291 cm/min，本文提出的模型模拟的结果与试验所测结果基本一致，证明了本文提出的模型及相应的程序在模拟非饱和固结问题的适用性。

图 5 Liakopoulos砂柱试验结果与模拟结果

Figure 5. Laboratory and simulated results of Liakopoulos sand column tests

下载: 全尺寸图片幻灯片

5.2 溶质迁移及热量传递COMSOL验证

溶质迁移针对一个长度为10 m的土柱，其左侧溶质质量浓度为1000 mg/L，右侧溶质质量浓度为0 mg/L，在固定压力梯度下溶质会在土柱内迁移扩散。初始时刻，土柱内溶质质量浓度为0 mg/L，对流速度v_lr采用固定流速，忽略降解对溶质质量浓度的影响，模拟溶质迁移行为。计算参数如表 1所示。

表 1 污染物迁移案例模型参数

Table 1. Model parameters for contaminant transportation case

孔隙率n	饱和度S_l	液相流速v_lr/ (m·s^-1)	扩散系数D^s _l/ (m²·s^-1)	孔隙弯曲因子τ_l	质量源 BQ^s _l/ (kg·m^-3)
0.6	0.4	4.8×10^-7	1.2×10^-9	0.1576	0

下载: 导出CSV

| 显示表格

图 6为溶质质量浓度在不同时刻的空间分布曲线。从图 6中可以发现，在对流和扩散的作用下溶质向右侧迁移，使内部液相溶质质量浓度从左至右逐渐升高，在第20天时，溶质已运移至8 m处，通过对比分析，Comsol有限单元法计算结果与计算程序结果一致，证明了溶质迁移模块及有限体积法求解的准确性。

图 6 污染物质量浓度随时间变化关系图

Figure 6. Variation of contaminant concentration with time

下载: 全尺寸图片幻灯片

热量传递基准案例模型长度为10 m，初始时刻其温度为30℃，模型左右两侧边界采用固定温度边界，左侧边界处温度为50℃，右侧为30℃，其传热过程采用多孔介质传热方程描述。考虑固相骨架的热传导和孔隙中液气的热对流作用，其内部气相组分假设为CH₄和CO₂等体积混合，模型参数如表 2所示。

表 2 热量传递案例模型参数

Table 2. Model parameters for heat transfer case

孔隙率n	饱和度S_l	固相密度ρ_s/ (kg·m^-3)	液相密度ρ_l/ (kg·m^-3)	气相密度ρ_g/ (kg·m^-3)	固相单位热容量H_s/ (J·kg^-1·K^-1)	液相单位热容量H_w/ (J·kg^-1·K^-1)	气相单位热容量H_g/ (J·kg^-1·K^-1)	液相流速v_lr/ (m·s^-1)	气相流速v_gr/ (m·s^-1)	热传导系数Г/ (m·s^-1·K^-1)	参考温度T_r/℃
0.6	0.4	800	1000	1.29	1300	4200	1514	10^-6	10^-6	0.436	20

下载: 导出CSV

| 显示表格

图 7给出了土体温度随时间及空间的变化图，左侧边界处温度较高，在骨架热传导和内部流体作用下热量从左至右传递。在第20天时，左侧3 m范围内温度与边界相同，处于等温状态，热量持续向右侧传递。多场耦合数值模型中的传热模块与数值软件Comsol计算结果一致，以上分析过程保证了多场耦合模型及计算程序的计算精度。

图 7 温度随时间变化关系图

Figure 7. Variation of temperature with time

下载: 全尺寸图片幻灯片

5.3 生活垃圾CAR1单元模型试验验证

南安普顿环境土工研究组在2008年进行了垃圾体填埋柱试验（CAR1试验^[34-35]），可作为生活垃圾多场耦合模型验证的基准试验。基于上述模型，引入生活垃圾两阶段降解模型^[36]，模拟CAR1试验的结果。南安普顿CAR1试验装置及示意图如图 8所示，其直径为480 mm，高度为900 mm。试验开始前，先在反应柱内装填粉碎的干生活垃圾，装填后生活垃圾高度约为0.8 m，质量约为27 kg，其后将90 L混合厌氧污泥（体积分数为10%）的渗滤液加入至反应柱中促使生活垃圾降解，混合后固液总体积为114.2 L。未加压条件下，垃圾体密度为237 kg/m³，此时试样的孔隙率为0.788。试验过程中将收集的渗滤液不断回灌至反应柱中，试验过程中保持室温为30℃。试验历时919 d，过程中记录堆体位移、产气速率、气相成分、渗滤液成分、渗滤液水位、pH等参数。

图 8 CAR1试验装置及示意图^[34]

Figure 8. Device and schematic diagram for CAR1 model test^[34]

下载: 全尺寸图片幻灯片

引入两阶段生化降解模型^[7]并嵌入数值程序，采用表 3所示的生活垃圾组成，对CAR1试验进行模拟。

表 3 CAR1试验中生活垃圾成分表

Table 3. Composition of municipal solid waste in CAR1 experiment

成分	纤维素^a	纤维素^b	糖类	蛋白质	脂肪	其他
质量/(kg·m^-3)	37.37	44.19	8.49	4.59	3.41	339.45
质量分数/%	8.54	10.12	1.94	1.05	0.78	77.57
注：^a代表快速降解纤维素，^b代表慢速降解纤维素。

下载: 导出CSV

| 显示表格

在顶部施加一大小为150 kPa的压力，达到平衡后获得模型的初始应力场，如图 9所示。

图 9 CAR1试验模型初始竖向应力分布图

Figure 9. Distribution of initial vertical stress of model for CAR1 tests

下载: 全尺寸图片幻灯片

初始条件下，渗滤液中VFA（挥发性脂肪酸）质量浓度假设为0，MB（甲烷菌质量浓度）质量浓度假设为0.1 g/L。试验测得的NH₄⁺离子质量浓度为0.189 g/L，快速降解和慢速降解组分的反应速率常数分别取0.02，0.0057 d，酸化和甲烷化抑制常数和指数均取2 kg/m³和8，初始pH值为7。对于模型顶部，采用固定液相水头回灌，气相压力为大气压力边界，位移边界为自由位移边界；对于模型底部，液相和气相压力均为零，以模拟自由排水边界，位移边界为零位移边界；对于四周，设为不渗透、固定法向位移边界，模型模拟参数见文献[7]。

图 10所示为CAR1试验结果与模型计算结果的对比。可以看出，本文计算的CAR1试验的产气量低于试验实测值，差异主要是模型中垃圾体可降解组分的估算不确定造成的，且试验过程中气体的密度未能测得，对于气体体积的估算也存在影响。

图 10 CAR1试验结果与模拟结果

Figure 10. Test and simulated results of CAR1 tests

下载: 全尺寸图片幻灯片

总体上，本文的模型及本文采取的垃圾体组分比例可以用于估算CAR1试验的产气量。模型的可计算pH值对生化降解的影响，同时可考虑铵根离子对溶液pH的影响，计算的pH与试验结果基本一致。对于CAR1试验，在150 kPa的压力下，反应柱高度从0.80 m降至0.341 m，之后降低至0.250 m。从图 10中可以看出，模型计算沉降值略低于实测堆体变形值，但总体上与试验结果一致。此外，慢速降解的纤维素质量分数的随时间变化的速率最小，符合模型实际情况，约在800 d左右，慢速降解的纤维素基本完成降解。

6. 结论

（1）建立了固相可降解土体降解–渗流–相变–传热–变形多场耦合理论模型。模型主要控制方程建立在多孔介质物质质量守恒、动量守恒、能量守恒基础上，模型液气两相运移严格建立在固结理论基础上。

（2）基于有限体积法进行模型的空间和时间离散，采用顺序耦合迭代方法对耦合模型进行数值求解。基于开场操作软件OpenFOAM编写了相应的数值求解程序，包括求解器、模型库和边界条件，可在其基础上实现不同降解模型、本构模型的嵌入。

（3）采用耦合模型模拟Liakopoulos砂柱排水试验，验证了模型计算土体渗流固结问题的正确性。通过模拟污染物对流-扩散和热量传递算例，将模拟结果与Comsol计算结果对比，验证了溶质迁移和热量传递模块的正确性。采用耦合模型模拟垃圾体降解CAR1试验，验证了耦合模型在模拟生活垃圾体多场耦合行为方面的合理性。

需要指出的是，本文模型进行了一定的简化和假设。后续应深入研究不同可降解土体的有效应力原理、本构关系、降解源项、固相质量损失与孔隙比关系等。

图 1 电动修复装置示意图

Figure 1. Schematic of electrokinetic remediation device

下载: 全尺寸图片幻灯片

图 2 阴极室中镉浓度随时间的变化

Figure 2. Variation of cadmium concentration with time in cathode chambers

下载: 全尺寸图片幻灯片

图 3 电压和镉初始浓度对电动修复效率的交互影响及电动修复后镉在土柱中的分布情况

Figure 3. Interactive effects of voltage and initial cadmium concentration on efficiency of electrokinetic remediation and distribution of cadmium in soil columns after electric remediation

下载: 全尺寸图片幻灯片

图 4 镉污染土修复前后SEM图像

Figure 4. SEM images of cadmium-contaminated soil before and after remediation

下载: 全尺寸图片幻灯片

图 5 各种形态镉含量及各种形态镉所占比例

Figure 5. Concentrations and proportions of various cadmium species

下载: 全尺寸图片幻灯片

[1]	YAO W K, CAI Z P, SUN S Y, et al. Electrokinetic-enhanced remediation of actual arsenic-contaminated soils with approaching cathode and Fe⁰ permeable reactive barrier[J]. Journal of Soils and Sediments, 2020, 20(3): 1526-1533. doi: 10.1007/s11368-019-02459-4
[2]	BARSOVA N, YAKIMENKO O, TOLPESHTA I, et al. Current state and dynamics of heavy metal soil pollution in Russian Federation: a review[J]. Environmental Pollution, 2019, 249: 200-207. doi: 10.1016/j.envpol.2019.03.020
[3]	ACAR Y B, GALE R J, ALSHAWABKEH A N, et al. Electrokinetic remediation: basics and technology status[J]. Journal of Hazardous Materials, 1995, 40(2): 117-137. doi: 10.1016/0304-3894(94)00066-P
[4]	GOMES H I, DIAS-FERREIRA C, RIBEIRO A B. Electrokinetic remediation of organochlorines in soil: enhancement techniques and integration with other remediation technologies[J]. Chemosphere, 2012, 87(10): 1077-1090. doi: 10.1016/j.chemosphere.2012.02.037
[5]	付志平, 冼子良, 雷畅, 等. 三种连续提取法对矿区土壤和尾砂中Pb形态分析研究[J]. 广东化工, 2022, 49(6): 80-83. https://www.cnki.com.cn/Article/CJFDTOTAL-GDHG202206018.htm FU Zhiping, XIAN Ziliang, LEI Chang, et al. Speciation analysis of Pb in soils and tailings surrounding mine area by three sequential extraction procedures[J]. Guangdong Chemical Industry, 2022, 49(6): 80-83. (in Chinese) https://www.cnki.com.cn/Article/CJFDTOTAL-GDHG202206018.htm
[6]	薛浩, 孟凡生, 王业耀, 等. 酸化-电动强化修复铬渣场地污染土壤[J]. 环境科学研究, 2015, 28(8): 1317-1323. https://www.cnki.com.cn/Article/CJFDTOTAL-HJKX201508021.htm XUE Hao, MENG Fansheng, WANG Yeyao, et al. Remediation of chromium residue-contaminated soil by preacidification electrokinetic remediation[J]. Research of Environmental Sciences, 2015, 28(8): 1317-1323. https://www.cnki.com.cn/Article/CJFDTOTAL-HJKX201508021.htm
[7]	李亚林, 刘蕾, 倪明, 等. 电动修复法对土壤中镉的去除及形态转化研究[J]. 工业安全与环保, 2016, 42(7): 76-79, 90. https://www.cnki.com.cn/Article/CJFDTOTAL-GYAF201607023.htm LI Yalin, LIU Lei, NI Ming, et al. Study on cadmium removal and speciation transformation from soil by electrokinetic remediation[J]. Industrial Safety and Environmental Protection, 2016, 42(7): 76-79, 90. (in Chinese) https://www.cnki.com.cn/Article/CJFDTOTAL-GYAF201607023.htm