基于虚拟加载上限法和黏土小应变特性的桩基<i>p</i>-<i>y</i>曲线

俞剑; 朱俊霖; 黄茂松; 沈侃敏

doi:10.11779/CJGE202111009

基于虚拟加载上限法和黏土小应变特性的桩基p-y曲线 English Version

俞剑^{1, 2,},
朱俊霖^{1, 2},
黄茂松^{1, 2, ,},
沈侃敏³

1.
同济大学地下建筑与工程系,上海　200092
2.
同济大学岩土及地下工程教育部重点实验室,上海　200092
3.
中国电建集团华东勘测设计研究院有限公司,浙江　杭州　310058

基金项目:

国家自然科学基金项目 51908420

国家自然科学基金项目 51579177

详细信息

作者简介:
俞剑(1987— ),男,助理研究员,博士,主要从事岩土工程方面的科研和教学工作。E-mail：002yujian@tongji.edu.cn

通讯作者:
黄茂松, E-mail：mshuang@tongji.edu.cn

中图分类号: TU470
计量
- 文章访问数: 340
- HTML全文浏览量: 52
- PDF下载量: 226
出版历程
- 收稿日期: 2021-03-08
- 网络出版日期: 2022-12-01
- 刊出日期: 2021-10-31

T-EMSD-based p-y curve of laterally loaded piles in clay considering small-strain behavior

YU Jian^{1, 2,},
ZHU Jun-lin^{1, 2},
HUANG Mao-song^{1, 2, ,},
SHEN Kan-min³

1.
Department of Geotechnical Engineering, Tongji University, Shanghai 200092, China
2.
Key Laboratory of Geotechnical and Underground Engineering of Ministry of Education, Tongji University, Shanghai 200092, China
3.
Power-China Huadong Engineering Corporation Limited, Hangzhou 310058, China

摘要

摘要: 国际主流规范API推荐p-y曲线分析海上风电钢管桩基础的水平非线性响应。该p-y曲线通过单参数刻画了土体应力应变的发展规律从而间接控制桩土非线性响应。但因过度简化,在长期使用中已暴露出低估桩侧承载力,且无法准确预测桩土初始刚度等问题。为此,首先定义了含土体小应变特性的应力应变曲线,再基于虚拟加载上限法获得考虑该土体应力应变关系的数值p-y骨干曲线。结合参数分析,拟合得到二维p-y骨干曲线表达式。进一步结合三维极限承载力系数、桩土初始刚度和剪应变系数,将p-y骨干曲线拟三维化,并分别与三维有限元和工程实例对比验证其合理性。与API规范相比,该p-y骨干曲线不但准确预测桩侧承载力,对变形控制极严的海上风电钢管桩基础而言,其更显著优点在于通过考虑土体小应变特性实现准确预测桩土初始刚度。
- 水平受荷桩 /
- 软黏土 /
- p-y曲线 /
- 虚拟加载上限法 /
- 小应变
Abstract: The international popular API code recommends the p-y curve method to analyze the nonlinear behavior of offshore wind turbine (OWT) steel pipe piles. The p-y curve controls the pile-soil nonlinear response only through one parameter regarding the development of stress-strain relation of soils. This over-simplification results in the inaccurate evaluation of the lateral initial stiffness of pile-soil and the underestimation of the bearing capacity. Therefore, the stress-strain curve with soil small-strain behavior is first introduced to achieve a numerical p-y backbone curve by using the total-displacement-loading extended mobilized strength design method (T-EMSD). The expression for the two-dimensional p-y backbone curve is then fitted from the numerical results. The three-dimensional effect of the proposed p-y curve is further considered by incorporating the three-dimensional ultimate capacity factor, the initial subgrade modulus and the compatibility factor. The rationality of the proposed p-y curve is verified against the results from the three-dimensional finite-element analysis and field tests. Compared with API code, the proposed p-y curve can provide a more reasonable prediction for both the bearing capacity and the initial stiffness of pile-soil by considering the soil small-strain behavior, which is a significant advantage for the OWT pile foundation with strict deformation control.
- laterally loaded pile /
- soft clay /
- p-y curve /
- T-EMSD method /
- small strain

HTML全文

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 T-EMSD迭代示意图

Figure 1. Iteration diagram of T-EMSD method

下载: 全尺寸图片幻灯片

图 2 桩身荷载-位移曲线

Figure 2. Load-displacement curves of fully bounded pile

下载: 全尺寸图片幻灯片

图 3 分段示意图

Figure 3. Piecewise diagram of p-y curve

下载: 全尺寸图片幻灯片

图 4 拟合p-y曲线与T-EMSD方法对比

Figure 4. Comparison of p-y curves by fitting and T-EMSD

下载: 全尺寸图片幻灯片

图 5 水平受荷桩拟三维分析

Figure 5. Construction of p-y curve considering three-dimensional effects

下载: 全尺寸图片幻灯片

图 6 土体参数沿深度变化

Figure 6. Variation of soil parameters along depth

下载: 全尺寸图片幻灯片

图 7 p-y曲线对比（Winkler方法与有限元）

Figure 7. Comparison of p-y curves between Winkler and FEM

下载: 全尺寸图片幻灯片

图 8 现场试验1的p-y曲线对比

Figure 8. Comparison of p-y curves of field test Case 1

下载: 全尺寸图片幻灯片

图 9 现场试验2的p-y曲线对比

Figure 9. Comparison of p-y curves of field test Case 2

下载: 全尺寸图片幻灯片

表 1 剪应变系数 $M_{c}$

Table 1 Values of compatibility factor $M_{c}$

上限解	完全胶结	完全光滑
Klar^[18]（MSD：单层）	1.30	1.35
黄茂松等^[6]（T-EMSD：100层）	0.80	0.90

下载: 导出CSV

表 2 现场试验1土体参数

Table 2 Soil parameters of field test Case 1

深度z/D	$s_{u}$ /kPa	$G_{0}$ /MPa	$γ_{0.7}$ /(10^-4)	$E_{50}$ /MPa
2	18	10	1	1.8
4	20	14	1	2

下载: 导出CSV

表 3 现场试验2土体参数

Table 3 Soil parameters of field test Case 2

深度z/m	$s_{u}$ /kPa	$G_{0}$ /MPa	$γ_{0.7}$ /(10^-4)	$E_{50}$ /MPa
3.6	13	12	4	0.87
6.3	15	15	4	1.0

下载: 导出CSV

参考文献(29)

[1]	DNV. Design of Offshore Wind Turbine Structures: DNV—OS—J101[S]. 2014.
[2]	American Petroleum Institute (API). Recommended practice for planning, designing and constructing fixed offshore platforms—working stress design[C]//API 2A-WSD, twenty-second ed. Washington, D.C, 2014.
[3]	MATLOCK H S. Correlation for design of laterally loaded piles in soft clay[C]//Proceedings Second Annual Offshore Technology Conference. Houston, 1970: 1204.
[4]	朱斌, 熊根, 刘晋超, 等. 砂土中大直径单桩水平受荷离心模型试验[J]. 岩土工程学报, 2013, 35(10): 1807-1815. https://www.cnki.com.cn/Article/CJFDTOTAL-YTGC201310007.htm ZHU Bin, XIONG Gen, LIU Jin-chao, et al. Centrifuge modelling of a large-diameter single pile under lateral loads in sand[J]. Chinese Journal of Geotechnical Engineering, 2013, 35(10): 1807-1815. (in Chinese) https://www.cnki.com.cn/Article/CJFDTOTAL-YTGC201310007.htm
[5]	王卫, 闫俊义, 刘建平. 基于海上风电试桩数据的大直径桩p-y模型研究[J]. 岩土工程学报, 2021, 43(6): 1131-1138. https://www.cnki.com.cn/Article/CJFDTOTAL-YTGC202106023.htm WANG Wei, YANG Jun-yi, LIU Jian-ping. Study on p-y models of large-diameter pile foundation based on in-situ tests of offshore wind power[J]. Chinese Journal of Geotechnical Engineering, 2021, 43(6): 1131-1138. (in Chinese) https://www.cnki.com.cn/Article/CJFDTOTAL-YTGC202106023.htm
[6]	黄茂松, 俞剑, 张陈蓉. 基于应变路径法的黏土中水平受荷桩p-y曲线[J]. 岩土工程学报, 2015, 37(3): 400-409. https://www.cnki.com.cn/Article/CJFDTOTAL-YTGC201503003.htm HUANG Mao-song, YU Jian, ZHANG Chen-rong. p-y curves of laterally loaded piles in clay based on strain path approach[J]. Chinese Journal of Geotechnical Engineering, 2015, 37(3): 400-409. (in Chinese) https://www.cnki.com.cn/Article/CJFDTOTAL-YTGC201503003.htm
[7]	ASHOUR M, NORRIS G. Modeling lateral soil-pile response based on soil-pile interaction[J]. Journal of Geotechnical and Geoenvironmental Engineering, 2000, 126(5): 420-428. doi: 10.1061/(ASCE)1090-0241(2000)126:5(420)
[8]	OSMAN A S, BOLTON M D. Simple plasticity-based prediction of the undrained settlement of shallow circular foundations on clay[J]. Géotechnique, 2005, 55(6): 435-447. doi: 10.1680/geot.2005.55.6.435
[9]	YU J, HUANG M S, LI S, et al. Load-displacement and upper-bound solutions of a loaded laterally pile in clay based on a total-displacement-loading EMSD method[J]. Computers and Geotechnics, 2017, 83: 64-76. doi: 10.1016/j.compgeo.2016.10.025
[10]	YU J, ZHU J L, SHEN K M, et al. Bounding-surface-based p-y model for laterally loaded piles in undrained clay[J]. Ocean Engineering, 2020, 216: 107997. doi: 10.1016/j.oceaneng.2020.107997
[11]	BENZ T, VERMEER P A, SCHWAB R. A small-strain overlay model[J]. International Journal for Numerical and Analytical Methods in Geomechanics, 2009, 33(1): 25-44. doi: 10.1002/nag.701
[12]	KLAR A, OSMAN A S. Load-displacement solutions for piles and shallow foundations based on deformation fields and energy conservation[J]. Géotechnique, 2008, 58(7): 581-589. doi: 10.1680/geot.2008.58.7.581
[13]	黄茂松, 李森, 俞剑. 水平受荷桩的弹性有限元虚拟加载上限分析[J]. 岩土力学, 2016, 37(8): 2399-2403, 2410. https://www.cnki.com.cn/Article/CJFDTOTAL-YTLX201608036.htm HUANG Mao-song, LI Sen, YU Jian. Analysis of laterally loaded pile by elastic finite element based EMSD method[J]. Rock and Soil Mechanics, 2016, 37(8): 2399-2403, 2410. (in Chinese) https://www.cnki.com.cn/Article/CJFDTOTAL-YTLX201608036.htm
[14]	LI S, YU J, HUANG M S, et al. Application of T-EMSD based p-y curves in the three-dimensional analysis of laterally loaded pile in undrained clay[J]. Ocean Engineering, 2020, 206: 107256. doi: 10.1016/j.oceaneng.2020.107256
[15]	LI S, HUANG M S, YU J. Continuous field based upper-bound analysis for the undrained bearing capacity of strip footings resting near clay slopes with linearly increased strength[J]. Computers and Geotechnics, 2019, 105: 168-182. doi: 10.1016/j.compgeo.2018.10.002
[16]	LI S, YU J, HUANG M S, et al. Upper bound analysis of rectangular surface footings on clay with linearly increasing strength[J]. Computers and Geotechnics, 2021, 129: 103896. doi: 10.1016/j.compgeo.2020.103896
[17]	HUANG M S, LI S, YU J, et al. Continuous field based upper bound analysis for three-dimensional tunnel face stability in undrained clay[J]. Computers and Geotechnics, 2018, 94: 207-213. doi: 10.1016/j.compgeo.2017.09.014
[18]	KLAR A. Upper bound for cylinder movement using “elastic” fields and its possible application to pile deformation analysis[J]. International Journal of Geomechanics, 2008, 8(2): 162-167. doi: 10.1061/(ASCE)1532-3641(2008)8:2(162)
[19]	USACE. Settlement analysis, Engineer Manual EM 1110-1-1904. U.S. Army Corps of Engineers: Washington, D.C. 1990.
[20]	GEORGIADIS M, ANAGNOSTOPOULOS C, SAFLEKOU S. Cyclic lateral loading of piles in soft clay[J]. Geotechnical Engineering, 1992, 23(1): 47-60.
[21]	JEANJEAN P. Re-assessment of p-y curves for soft clays from centrifuge testing and finite element modeling[C]//Offshore Technology Conference. Houston, 2009.
[22]	DUNNAVANT T W, O'NEILL M W. Experimental p-y model for submerged, stiff clay[J]. Journal of Geotechnical Engineering, 1989, 115(1): 95-114. doi: 10.1061/(ASCE)0733-9410(1989)115:1(95)
[23]	BYRNE B W, HOULSBY G T, BURD H J, et al. PISA design model for monopiles for offshore wind turbines: application to a stiff glacial clay till[J]. Géotechnique, 2020, 70(11): 1030-1047. doi: 10.1680/jgeot.18.P.255
[24]	ZHANG C R, YU J, HUANG M S. Winkler load-transfer analysis for laterally loaded piles[J]. Canadian Geotechnical Journal, 2016, 53(7): 1110-1124. doi: 10.1139/cgj-2015-0394
[25]	YU J, HUANG M S, ZHANG C R. Three-dimensional upper-bound analysis for ultimate bearing capacity of laterally loaded rigid pile in undrained clay[J]. Canadian Geotechnical Journal, 2015, 52(11): 1775-1790. doi: 10.1139/cgj-2014-0390
[26]	ZHANG Y H, ANDERSEN K H, TEDESCO G. Ultimate bearing capacity of laterally loaded piles in clay-Some practical considerations[J]. Marine Structures, 2016, 50: 260-275. doi: 10.1016/j.marstruc.2016.09.002
[27]	FAN C C, LONG J H. Assessment of existing methods for predicting soil response of laterally loaded piles in sand[J]. Computers and Geotechnics, 2005, 32(4): 274-289. doi: 10.1016/j.compgeo.2005.02.004
[28]	KIM Y, JEONG S, LEE S. Wedge failure analysis of soil resistance on laterally loaded piles in clay[J]. Journal of Geotechnical and Geoenvironmental Engineering, 2011, 137(7): 678-694. doi: 10.1061/(ASCE)GT.1943-5606.0000481
[29]	SUN C G, CHO C S, SON M, et al. Correlations between shear wave velocity and in-situ penetration test results for Korean soil deposits[J]. Pure and Applied Geophysics, 2013, 170(3): 271-281. doi: 10.1007/s00024-012-0516-2

施引文献

资源附件(0)

图(9) / 表(3)

计量

文章访问数: 340
HTML全文浏览量: 52
PDF下载量: 226
被引次数: 0

0. 引言
1. 多场耦合理论模型
1.1 耦合模型基本物理力学关系
1.2 多场耦合模型单元体守恒方程
1.3 封闭关系
2. 耦合模型本构关系
2.1 可降解土体降解相变动力学
2.2 可降解土体应力应变关系
2.3 可降解土体的热力参数
2.4 可降解土体的持水及扩散特性
3. 耦合模型控制方程
4. 耦合模型有限体积离散及程序实现
4.1 控制方程顺序迭代思路
4.2 有限体积离散
4.3 求解步骤
4.4 数值计算程序
5. 模型及程序验证
5.1 Liakopoulos砂柱排水固结试验
5.2 溶质迁移及热量传递COMSOL验证
5.3 生活垃圾CAR1单元模型试验验证
6. 结论

0. 引言
1. 多场耦合理论模型
1.1 耦合模型基本物理力学关系
1.2 多场耦合模型单元体守恒方程
1.3 封闭关系
2. 耦合模型本构关系
2.1 可降解土体降解相变动力学
2.2 可降解土体应力应变关系
2.3 可降解土体的热力参数
2.4 可降解土体的持水及扩散特性
3. 耦合模型控制方程
4. 耦合模型有限体积离散及程序实现
4.1 控制方程顺序迭代思路
4.2 有限体积离散
4.3 求解步骤
4.4 数值计算程序
5. 模型及程序验证
5.1 Liakopoulos砂柱排水固结试验
5.2 溶质迁移及热量传递COMSOL验证
5.3 生活垃圾CAR1单元模型试验验证
6. 结论

参考文献(29)

施引文献

资源附件(0)

基于虚拟加载上限法和黏土小应变特性的桩基p-y曲线 English Version

作者简介: 俞剑(1987— ),男,助理研究员,博士,主要从事岩土工程方面的科研和教学工作。E-mail：002yujian@tongji.edu.cn

通讯作者: 黄茂松, E-mail：mshuang@tongji.edu.cn

计量

出版历程

T-EMSD-based p-y curve of laterally loaded piles in clay considering small-strain behavior

0. 引言

1. 多场耦合理论模型

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

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

1.3 封闭关系

2. 耦合模型本构关系

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

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

2.3 可降解土体的热力参数

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

3. 耦合模型控制方程

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

4.1 控制方程顺序迭代思路

4.2 有限体积离散

4.3 求解步骤

4.4 数值计算程序

5. 模型及程序验证

5.1 Liakopoulos砂柱排水固结试验

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

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

6. 结论

计量

出版历程

目录

0. 引言

1. 多场耦合理论模型

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

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

1.3 封闭关系

2. 耦合模型本构关系

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

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

2.3 可降解土体的热力参数

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

3. 耦合模型控制方程

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

4.1 控制方程顺序迭代思路

4.2 有限体积离散

4.3 求解步骤

4.4 数值计算程序

5. 模型及程序验证

5.1 Liakopoulos砂柱排水固结试验

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

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

6. 结论

作者简介:
俞剑(1987— ),男,助理研究员,博士,主要从事岩土工程方面的科研和教学工作。E-mail：002yujian@tongji.edu.cn

通讯作者:
黄茂松, E-mail：mshuang@tongji.edu.cn