气水两相流现象及其影响广泛存在于各类岩土工程中. 19世纪初期，King^[1]阐述了孔隙气压力相对于大气压变化的响应现象. Richards^[2]将Darcy定理应用到饱和-非饱和渗流研究中. 传统的水气两相流理论认为非饱和土壤中只有水的运移，假设孔隙气不影响水的流动，此种假设也被称为Richards假设^[3]，相应的方程称为Richards方程. 该理论模型为之后的渗流分析奠定了基础.

在水气二相流问题上，国内外学者大多从Richards方程出发，多是集中于单相流模型^[4-6]，对气水两相流以及气液固三相耦合的研究较少^[7-9]. 研究发现，当介质不均匀或在极短时间内渗流量较大时，则不可以忽视气相的存在，气相的存在会对水相的流动产生影响. 彭胜等^[10]进行了相关的室内试验，分析认为气体的运动方向主要沿着进水口和排气口之间的直线方向，在其他方向上气相的运动微乎其微.

在渗流计算领域，由于TOUGH2强大的多相流多组分计算能力，越来越多的模型跟TOUGH2耦合，以满足不同领域的气液固耦合计算需求. Rutqvist等^[11]将带有ECO2N模块的TOUGH与FLAC3D进行耦合，用在二氧化碳地质封存领域，开创了基于TOUGH2开展流固耦合模拟的先河. Gosavi 等^[12]将TOUGH2与GeoCrack3D耦合，以研究离散裂缝中的水力热耦合问题；Hurwitz等^[13]将TOUGH2与开源有限元软件Biot耦合，主要用于水热系统研究，包括与火山相关的模拟；Taron等^[14]将FLAC3D与TOUGHREACT耦合，以解决增强性地热系统开发领域的问题；Rohmer等^[15]将TOUGH2与开源有限元程序CodeAster耦合，在法国二氧化碳地质封存的地质调查过程中，用于分析力学相关的部分；Kim等^[16]将TOUGH2与ROCMECH有限元程序耦合，用于页岩气压裂研究；Pan等^[17]将TOUGH2与RDCA耦合，以研究多孔介质中的裂纹扩展问题；Lee等^[18]将TOUGH2与UDEC耦合，以研究裂隙岩体的力学问题；Lei等^[19]将TOUGH2与多孔弹性介质有限元程序耦合，以分析增强型地热系统与二氧化碳地质封存领域的问题；Park等^[20]通过TOUGH-FLAC模拟器研究在低渗透性岩石中注入流体引发断层活化的问题，研究表明该方法具备捕捉断层中流体力学过程的能力. An等^[21]将TOUGH2的EOS4模块中蒸气压变量用温度变量代替，开展典型油藏流动模拟研究，最后通过两相流模拟进行了验证. 可以看出，TOUGH2出色的渗流计算能力，得到了众多研究者的青睐，将TOUGH2与其他软件联合应用的案例越来越多，TOUGH2与其他软件的耦合程序的开发受到越来越多的研究者的关注.

Richards方程忽略了孔隙气压对水的流动的影响，在具体的应用中将孔隙气压力假设为恒定的大气压，即假设孔隙气压是与外界大气压连通的. 这不符合实际情况，也没有考虑结构体中各相介质之间的相互作用，也就是说，在物理机制上，无法探究气水渗流与固体颗粒之间的相互作用关系，其分析结果具有一定的局限性. 另外，越来越多的研究表明，非饱和土中孔隙气能够明显阻滞和加速土壤中水的入渗过程^[22-24]. 因此，有必要对渗流场中气水两相渗流进行详细研究，以加深对气水两相渗流机理的认识，准确分析渗流场对实际工程的影响.

本研究根据非饱和土降雨入渗过程中气水两相渗流-应力弱耦合理论，基于TOUGH2（EOS3和EOS9模块）与FLAC3D，开发了基于Python语言的气水两相渗流-应力弱耦合计算搭接程序. 该计算程序能够模拟真实的气水两相流，能够探讨渗流过程中气水的相互作用及其对渗流过程的影响. 程序考虑了气水两相渗流与土体骨架变形直接的相互作用，反映了这一过程中孔隙度、渗透率、毛管压力和土体物理力学参数的变化，实现了更为完善的气-水两相渗流与应力弱耦合分析.

1）质量守恒方程. 在TOUGH2中，通过流体中各组分的相对质量分数计算各自的质量守恒方程. 对应的质量守恒方程为

(1) $ \frac{\partial }{\partial t}{M}^{k }-{Q}^{k }=-\nabla \cdot ({\boldsymbol{q}}_{\mathrm{l}}^{k}+{\boldsymbol{q}}_{\mathrm{g}}^{k} ) \text{，} $

(2) $ {M}^{k }=\phi {S}_{\mathrm{l}}{\rho }_{\mathrm{l}}{w }_{{\rm{l}}}^{k}+\phi {S}_{\mathrm{g}}{\rho }_{\mathrm{g}}{w }_{{\rm{g}}}^{k} \text{，} $

(3) $ {\boldsymbol{q}}_{\psi }^{k}=-{\rho }_{\psi }{w}_{{\rm{g}}}^{k}\frac{\boldsymbol{k}{k}_{{\rm{r}}\psi }}{{\eta_\psi }}\left(\nabla {{p}}_{\psi }-{\rho }_{\psi }{g}\nabla z\right)+{{{\boldsymbol{i}}}}_{\mathrm{\psi }}^{{k}} \text{，} $

(4) $ {{{\boldsymbol{i}}}}_{\psi }^{k}=-{\rho }_{\psi }{D}_{\mathrm{v}}{{\boldsymbol{I}}}\nabla {w}_{\psi }^{k} . $

式中： $ \psi $为相， $ \psi =\left\{\mathrm{g},\mathrm{l}\right\}，\mathrm{g} $表示气相， $ \mathrm{l} $表示液相； $k$为组分， $ {M}^{k} $为组分 $ k $的单位体积质量， $ {Q}^{k} $为组分 $ k $的源汇相； ${\boldsymbol{q}}_{\psi }^{k} $为质量通量密度； ${\phi } $为孔隙度； $ S $为饱和度； $ \;{\rho }_{\psi } $为 $ \psi $相的密度；w为质量分数；k为固有渗透率张量； ${k}_{{\rm{r}}\psi }$为相对渗透率； $ {\eta }_{\psi } $为动力黏度； ${{p}}_{\psi }$为相的压力； ${g}$为重力加速度； $ z $为高程； ${{{\boldsymbol{i}}}}_{\psi }^{k}$为扩散质量通量； $ {D}_{\mathrm{v}} $为有效分子扩散系数； $ \boldsymbol{I} $为单位张量.

2）能量守恒方程.

(5) $ \frac{\partial }{\partial t}\left(\phi {S}_{{\rm{l}}}{\rho }_{{\rm{l}}}{e}_{{\rm{l}}}+\phi {S}_{{\rm{g}}}{\rho }_{{\rm{g}}}{e}_{{\rm{g}}}+\left(1-\phi \right){\rho }_{{\rm{s}}}{C}_{{\rm{s}}}T_0\right) -{Q}^{\mathrm{h}}=-\nabla \cdot \left({{{\boldsymbol{q}}}}^{\mathrm{h}}\right) , $

(6) $ {{\boldsymbol{q}}}^{\mathrm{h}}={\sum} _{\psi }{\sum} _{k}{h}_{\psi }^{k}{{\boldsymbol{q}}}_{\psi }^{k}+{{\boldsymbol{i}}}_{\mathrm{m}}^{\mathrm{h}} , $

(7) $ {{\boldsymbol{i}}}_{\mathrm{m}}^{\mathrm{h}}=-{\lambda }_{\mathrm{m}}{{\boldsymbol{I}}}\nabla T_0 . $

式中： ${e}$为内能， ${C}_{{\rm{s}}}$为固体的比热容， $T_0$为绝对温度， $ {Q}^{\mathrm{h}} $为热量源汇相， $ {\boldsymbol{q}}^{\mathrm{h}} $为能量通量密度， $ {h}_{\psi }^{k} $为比焓， $ {\boldsymbol{i}}_{\mathrm{m}}^{\mathrm{h}} $为表观热传导相， ${\lambda }_{{\rm{m}}}$为表观导热系数.

3）非线性求解的迭代方程.

(8) $ \begin{split} {R}_{{n}}^{\left(\psi \right){s}+1}=\;&{M}_{{n}}^{\left(\psi \right){s}+1}-{M}_{{n}}^{\left(\psi \right){s}}-\\ \;&\frac{\Delta t}{{{V}}_{{n}}}\left\{\sum _{{m}}{A}_{{n}{m}}{q}_{{n}{m}}^{\left(\psi \right){s}+1}+{V}_{{n}}{Q}_{{n}}^{\left(\psi \right){s}+1}\right\} . \end{split} $

式中： $ R $为残差； $ V $为单元体积； $ A $为单元面积，下标nm表示2个单元n、m交界面的权重平均值； $ \Delta t $为时间步；上标 ${s}$和 ${s}+1$为相邻时间步.

1）运动方程.

(9) $ \nabla \cdot \boldsymbol{\sigma }+{\rho }_{{\rm{m}}}{g}={\rho }_{{\rm{m}}}\frac{\mathrm{d}{{{v}}}}{\mathrm{d}\mathrm{t}} . $

式中： $ \boldsymbol{\sigma } $为应力张量， ${\;\rho }_{{\rm{m}}}$为岩体的均值密度， ${v}$为节点速度.

2）本构方程.

(10) $ \Delta {\sigma }'={H}({\sigma }',\dot{{\boldsymbol{\varepsilon }}}\Delta {t} ) \text{，} $

(11) $ \dot{{\boldsymbol{\varepsilon}} }=\frac{1}{2}(\nabla {{\boldsymbol{v}}}+{(\nabla {{\boldsymbol{v}}})}^{\mathrm{T}} ) \text{，} $

(12) $ {\boldsymbol{\varepsilon}} =\frac{1}{2}(\nabla {{\boldsymbol{u}}}+{(\nabla {{\boldsymbol{u}}})}^{\mathrm{T}} ) \text{，} $

(13) $ \Delta {\boldsymbol{\varepsilon }}=\Delta {{\boldsymbol{\varepsilon}} }^{{\rm{e}}}+\Delta {{\boldsymbol{\varepsilon}} }^{\mathrm{p}}+\Delta {{\boldsymbol{\varepsilon}} }^{\mathrm{T}} \text{，} $

(14) $ \Delta {\boldsymbol{\varepsilon }}^{\mathrm{T}}={{\boldsymbol{I}}}\mathrm{\alpha }\Delta T_0 . $

式中： ${H}$为给定的控制方程， $ \dot{{\boldsymbol{\varepsilon}} } $为应变率张量， $ {\boldsymbol{\varepsilon }}$为应变，α为热膨胀系数.

3）有效应力方程.

(15) $ \Delta {\boldsymbol{\sigma }}'=\boldsymbol{\sigma }+{{\boldsymbol{I}}}\mathrm{\beta }p . $

式中： $\; \mathrm{\beta } $为Biot’s有效应力系数， $p$为孔隙压力.

4）求解迭代计算方程.

(16) $ {{{\boldsymbol{v}}}}_{i}^{(t+\Delta t/2)}={{{\boldsymbol{v}}}}_{i}^{(t-\Delta t/2)}+\frac{\Delta t}{m}{\sum} _{i}{{{\boldsymbol{F}}}}_{i}^{t} . $

式中： $ m $为集中节点质量， $ \boldsymbol{F} $为节点力.

利用Python语言，基于弱耦合理论，开发了TOUGH2-FLAC3D气液固耦合计算程序. TOUGH2与FLAC3D耦合计算方法，最早由Rutqvist^[11]提出. 主要用于TOUGH2中ECO2N模块与FLAC3D耦合计算，应用在CO₂地质封存领域多场耦合分析中. 该程序充分利用TOUGH2和FLAC3D在各自领域的优势，取长补短，将两者耦合以达到研究目的. TOUGH2程序计算多孔介质中多相流多组分的流动，该程序已经被很多研究机构验证和使用^[25]. FLAC3D程序主要针对岩土力学领域，能够处理单相流动状态的热力学耦合和水力耦合过程^[26]. 两者在各自领域都有着广泛的应用.

本研究通过弱耦合理论基本原理，在计算处理上作如下假定：1）多孔介质由固体骨架、液体和气体三相组成. 2）除了蒸汽和气体以外，不考虑各相之间的相互溶解. 3）忽略各个相的交界面上的热力学效应.

TOUGH2-FLAC3D耦合计算的整体框架，如图1所示^[11].

在渗流力学耦合过程中，流体在固体骨架孔隙中的运移，将会引起温度场T、孔隙压力场p以及流体饱和度S的变化. 由于流固耦合的作用，这种变化会引起应力-应变场的改变，从而引起岩土体骨架变形，进而引起诸如渗透率、毛细管压力的改变，因此，可以通过修正、更新每一时步的孔隙率、渗透率和毛管压力进行迭代计算. 如此反复循环迭代计算，直至达到设定的收敛条件. 同时，渗流力学耦合作用会引起水位面的改变，饱和度的改变会使得物理力学参数改变，比如土体的有效黏聚力、内摩擦角.

本研究所提出的耦合程序的主要特点如下：1）基于多孔介质的流体力学理论，忽略化学因素的影响. 利用TOUGH2和FLAC3D计算程序，实现了水、力、热耦合计算. 2）主要使用了2种编程语言Python和FISH. 其中Python主要用于TOUGH2和FLAC3D的交互计算，FISH则用于FLAC3D的内部操作计算. 3）TOUGH2与FLAC3D的计算网格配套，只需要一套计算网格，无须额外建立计算网格. 4）对比表明，该耦合程序在计算流固耦合问题时，计算速度远远大于目前流行的商业程序. 尤其是对于大型工程或者网格较多的模型，两者的计算速度有数量级上的差别.

本耦合计算程序的基本理论，是用TOUGH2来计算流体（气-水）的流动和热传导，将计算后的孔隙压力（孔隙水压力和孔隙气压力）和温度传递给FLAC3D程序，来进行应力应变计算，之后更新孔隙度、渗透率和毛细管压力. 将更新后的结果返回TOUGH2程序，作为TOUGH2的输入参数进行相关计算. 如此循环迭代，直至收敛，如图2所示为该耦合计算程序的结构图.

在程序发布时，仅包括一个src目录和一个initial_tough_flac.bat文件. 在计算时，首先建立一个工作目录，在完成初始化之后，程序将建立一组工作目录，结构图如图2所示. common_data目录，用来存放程序运行时所需要的公共文件和部分临时文件，包括TOUGH2的气-水两相流计算程序T2_EOS3.exe和TOUGH2的单相流计算程序T2_EOS9.exe. src目录包含全部计算程序和一些说明文件. step1_get_MESH目录，主要用于建立TOUGH2和FLAC3D耦合计算的配套的单元网格模型. 本程序网格模型的转化，是从FLAC3D中建立网格模型，然后通过转换程序，转换成TOUGH2格式的网格文件，供TOUGH2直接调用进行相关的计算. step2_get_TOUGH2_Result目录，主要用于TOUGH2中的流体计算. 里面包含2个文件夹：INCON文件夹和FLOW文件夹. 其中，INCON文件主要存储TOUGH2程序计算后得到的初始条件，包括初始温度场、初始流体压力场以及初始气水饱和度分布场. FLOW主要存储流体流动计算的结果. step3_TOUGH2_to_FLAC3D文件主要存储TOUGH2和FLAC3D的耦合计算的结果. 值得注意的是，每个目录下面都包含有一个readme.txt和step_cmd.bat文件. readme.txt文件主要用来说明各个工作目录的计算功能以及操作步骤，step_cmd.bat文件是具体的DOS批处理运行程序.

由上述分析可知，TOUGH2-FLAC3D的耦合计算，实质上就是2个程序在计算过程中进行数据传递. 整体计算分为3部分：1）网格转换，这是耦合计算的基础；2）分别计算得到各自对应的初始状态；3）耦合计算.

如图3所示，TOUGH2的数据存储在单元中心，而FLAC3D的数据存储在网格节点，为了实现耦合计算，须将两者的网格信息转换一致. 本研究将FLAC3D中的网格转换成TOUGH2计算的配套网格格式. 由于TOUGH2能够处理的网格为泰森多边形网格，目前该程序能处理的模型单元为六面体网格.

这一步主要为TOUGH2-FLAC3D的耦合计算准备初始条件，即初始的流体压力、温度、流体饱和度及初始的应力应变场. 具体计算流程，如图4所示. 具体为基于设置的边界条件和初始状态参数，由TOUGH2计算渗流场部分，由FLAC3D计算应力场部分，最后调用耦合程序计算流固耦合部分.

为了克服收敛问题，岩土体变形和孔隙参数在模拟计算中单独计算. TOUGH2和FLAC3D耦合计算程序的迭代计算流程如图5所示.

FLAC3D程序从step2_get_TOUGH2_Result/FLOW文件夹中读取由TOUGH2计算得到的流体压力p、温度T、气体饱和度S_g或液体饱和度S_l，进行力学计算，并更新孔隙度、渗透率和毛管压力，之后将孔隙度、渗透率和毛管压力返回TOUGH2程序进行迭代计算.

由于要考虑骨架的变形造成的孔隙度、渗透率和毛管压力的更新，采用如下公式进行相应更新^[27-30]：

(17) $ \mathrm{d}\phi =\left(1-{\phi }_{0}\right){\varepsilon }_{\mathrm{v}}-(1-{\phi }_{0}){\varepsilon }_{\mathrm{s}} \text{，} $

(18) $ K={K}_{0}{\left(\frac{\phi }{{\phi }_{0}}\right)}^{3}{\left(\frac{{1-\phi }_{0}}{1-\phi }\right)}^{2} \text{，} $

(19) $ {p}_{{\rm{cL}}}={p}_{{\rm{c}}}\frac{\sqrt{{K}_{0}/{\phi }_{0}}}{\sqrt{K/\phi }} . $

式中 $ ：{\varepsilon }_{\mathrm{v}} $为体积应变增量， $ {\varepsilon }_{\mathrm{s}} $为剪切应变增量， $ {\phi }_{0} $为初始孔隙度， $ {K}_{0} $为初始渗透率， ${p}_{{\rm{c}}}$为初始毛管压力.

该耦合程序主要包括2大类型，单相流耦合计算程序TOUGH2(EOS9)-FLAC3D和气-水两相流耦合计算程序TOUGH2(EOS3)-FLAC3D. 两者之间的主要区别为TOUGH2的计算模块EOS3和EOS9的区别. EOS3模块在渗流计算中能够计算真实的气水两相流的运移及其相互作用的影响. EOS9模块在渗流计算中采用被动气压力假设，在计算中程序默认系统的气体与外界大气连通，渗流过程中气压始终保持标准大气压不变.

本研究开发的TOUGH2-FLAC3D耦合程序只使用了TOUGH2的EOS3模块和EOS9模块，主要针对两相渗流问题和单相渗流问题. TOUGH2系列软件包含众多的模块，很容易将TOUGH2其他模块与FLAC3D计算程序搭接起来，在搭接过程中，只要注意TOUGH2的主要计算变量，做好流体参数与力学参数的耦合传递以及相关渗透率、孔隙度和毛管压力的更新即可. 本研究所开发的程序充分发挥TOUGH2系列软件的优势，大大拓展TOUGH2-FLAC3D软件的使用范围，模拟岩土工程中涉及到的各种多场耦合作用.

采用著名的Liakopoulos排水实验^[31]来验证TOUGH2-FLAC3D气液固三相耦合分析程序的正确性. 很多研究人员采用Liakopoulos排水实验来验证气液固三相耦合模型的有效性^[32-39].

实验模型如图6所示. 该实验采用高为1m、直径为0.1m的Del Monte砂柱，沿其高度方向划分成20个网格，砂柱的底部铺设有一层透水石，侧面采用圆筒隔水. 实验在恒定的室温（35 ℃）下进行.

在试验开始前，将水流从试样顶部连续注入，试样底部的水流自由流出，直到水流稳定，即保证初始时刻试样为完全饱和状态，在此饱和状态建立好后试件顶部停止注水，然后抽去底部隔水板，使试样在重力的作用下自由排水，此后正式开始试验，记录不同时刻不同位置的试验结果.

在数值计算中，使用的毛管压力曲线和相对渗透率曲线，是由Liakopoulos^[31]拟合得到的. 表达式分别如下：

(20) $ {S}_{{\rm{l}}}=1.0-0.101\,52\times {\left(\frac{{p}_{{\rm{CAP}}}}{9.81}\right)}^{2.427\;9} \text{，} $

(21) $ {k}_{{\rm{rl}}}=1.0-2.207\times {(1.0-{S}_{{\rm{l}}})}^{1.012\;1} . $

式中： ${p}_{{\rm{CAP}}} $为毛管压力.

在Liakopoulos^[31]的工作中，并没有涉及到气体入渗试验，因此引用BROOKS and COREY模型^[40]来描述气相渗透率. 表达式如下：

(22) $ {k}_{{\rm{rg}}}={(1.0-{S}_{{\rm{e}}})}^{2}\times {(1.0-{S}_{{\rm{e}}})}^{5/3}+1.0\times {10}^{-4} , $

(23) $ {S}_{{\rm{e}}}=({S}_{{\rm{l}}}-{S}_{{\rm{rl}}})/(1.0-{S}_{{\rm{rl}}}) . $

式中： ${S}_{{\rm{e}}}$为水的有效饱和度， ${S}_{{\rm{rl}}} $为水的相对饱和度. 式（22）中的 $ 1.0\times {10}^{-4} $主要用来近似考虑饱和状态下Del Monte砂样中的微孔隙气体运移^[32,34].

在本研究试验中没有测量模型的力学参数，借鉴Schrefler 等^[33]的参数取值，具体参数如表1所示. 表中，E为弹性模量， $ \nu $为泊松比， $ \rho $为密度， $ \phi $为孔隙度， ${p}_{{\rm{atm}}}$为大气压， $ k $为渗透率.

根据Liakopoulos实验，在数值计算中，流体边界条件定义如下：1）砂样顶部为隔水、透气边界，即 ${p}_{{\rm{g}}}={p}_{{\rm{atm}}}$， ${w}_{{\rm{ga}}}=0.999$，其中 ${p}_{{\rm{g}}}$为气相压力， ${w}_{{\rm{ga}}}$为气相中空气所占的质量分数；2）砂样底部为透气、位移全约束边界，允许水流溢出，具体设置为 ${p}_{{\rm{w}}}={1.013\times 10}^{5}\;\mathrm{P}\mathrm{a}$， ${w }_{{\rm{la}}}=0.0$，其中 ${p}_{{\rm{w}}}$为液相压力， ${w }_{{\rm{la}}}$为液相中空气所占的质量分数；3）侧向为隔水、法向约束边界，边界条件定义为 ${p}_{{\rm{g}}}={p}_{{\rm{atm}}}$， ${w }_{{\rm{ga}}}=0.999$^[33]. 计算过程中砂样中的水在重力作用下下渗.

整个实验的模拟时间为120 min. 在通过TOUGH2(EOS3)-FLAC3D气液固耦合程序计算后，重点分析砂柱不同位置、不同时刻的孔隙水压力水头、沿高度方向不同时刻孔隙水压力水头、出口处的水流速度以及竖向位移. 其中，孔隙水压力水头的计算公式为 $({p}_{{\rm{w}}}-{p}_{{\rm{atm}}})/ \left({\rho }_{{\rm{w}}}\mathrm{g}\right)$；孔隙气压力水头的计算公式为 $({p}_{{{\rm{g}}}}-{p}_{{\rm{atm}}})/\left({\rho }_{{\rm{w}}}\mathrm{g}\right)$， $({p}_{{\rm{w}}}-{p}_{0})/ ({\rho }_{{\rm{w}}}g)$. 其中， ${\;\rho }_{{\rm{w}}}$为水的密度.

如图7所示为砂柱不同位置、不同时刻的孔隙水压力水头h随时间t的变化情况. 主要监测h=0.625 m和h=0.975 m这2处的孔隙水压力水头.可以看出，模拟值和Liakopoulos试验值在排水初期稍有误差，在排水过程稳定后，耦合程序计算得到的孔隙水压力水头基本能与实验测点的数据吻合.

如图8所示为砂柱沿高度H方向不同时刻孔隙水压力水头h，模拟了t=30、60、120 min这3个时间段. 可以看出，与Liakopoulos试验测点结果相比，在模拟前期t=0~30 min这一时间段内，可能由于砂柱内的水流尚未稳定，模拟结果与实测结果有所差别，但是总体变化趋势一致；随着模拟时间的延长，砂柱的排水趋于稳定，耦合计算程序计算的结果与实验测点数据基本吻合.

如图9所示为砂柱在出口处的水流速度v随时间的变化情况，模拟总时间为t=120 min. 本研究将Ehlers等^[37]和Nagel等^[38]的数值模拟数据汇总，与本程序的计算结果进行对比.可以看出，与Liakopoulos试验测点结果相比，在排水初期，Ehlers等^[37]、Nagel 等^[38]和本研究开发的耦合计算程序计算的出口流速均与实验测点有所偏差，可能是因为渗流没有达到稳定，目前的数值模型不能较好地模拟非稳定阶段. 随着模拟时间的延长，在达到稳定渗流后，本研究的计算程序所获得的结果与实验测点值更加接近.

为了研究降雨入渗过程中降雨强度和累积降雨量对边坡的影响，林鸿州^[41]通过研制的试验设备进行了相关的滑坡模拟试验，探讨不同情况下的滑坡特性. 针对其中一个试验，用本研究开发的TOUGH2-FLAC3D耦合搭接程序进行模拟，将模拟结果与试验结果进行对比，以验证程序的有效性和适用性.

选用其中一个试验进行模拟分析，如图10所示为该试验的模型侧视图和俯视图. 如图11所示为边坡试验的模型图以及张力计埋设位置图.

非饱和渗流数值计算的相对渗透率和毛管压力都采用van Genuchen模型. 具体的计算参数如表2所示. 表中， ${\gamma }_{{\rm{d}}}$为干容重， ${G}_{{\rm{s}}}$为比重，c为有效黏聚力， $\mathrm{\phi }_{\rm{e}}$为有效内摩擦角， ${k}_{{\rm{sat}}}$为饱和渗透系数，α、n为van Genuchen模型拟合参数， $ {\varphi }_{{\rm{s}}} $为水的饱和体积分数， $ {\varphi }_{{\rm{r}}} $为水的残余体积分数.

在该试验的数值模拟中，数值计算分2步进行. 首先通过计算获得边坡土体的初始状态，在此基础上再进行降雨入渗模拟. 在进行第1步计算时，流体边界条件定义如下：1）边坡顶部为大气压边界，即 ${p}_{{\rm{g}}}={p}_{{\rm{atm}}}$， ${w }_{{\rm{ga}}}=0.999$；2）边坡底部设置为 ${p}_{{\rm{w}}}={1.062\times 10}^{5}$ Pa， ${w }_{{\rm{la}}}=0.0$；3）根据试验模型认为侧边界为零流量边界. 在第2步降雨入渗模拟中，边坡的整个上边界为降雨边界，本实验的降雨强度为80 mm/h. 在模拟实验过程中力学边界条件定义为水平和竖直方向均为固定约束.

记录了本次模拟120 min内的降雨入渗过程中孔隙水压力p与时间t的关联数据，如图12所示. 在该试验中，林鸿州^[41]采用自编的非饱和渗流分析程序进行分析计算，其计算模型是基于Richards方程来进行计算的，虽然考虑了基质吸力，但是并没有考虑气体的流动. 在实际的计算过程中，使用被动气压力的假设，即认为渗流过程中土体内的气压力等于大气压，从而忽略了气压力的变化以及对孔隙水压力的影响. 因此，由图12可以看出，初始阶段计算的孔隙水压力的结果高于实际试验测得的孔隙水压力，直到22 min时，由于气体逐渐与大气压连通，这种差别才慢慢变小. 两相流计算的结果则与实际试验结果较接近，说明两相流模型能够较为真实地描述非饱和土体内部在降雨入渗过程中气液两相流体的流动状态，而且，气液固三相耦合计算也能较为准确地描述非饱和土体内部气水两相流动，因此，其分析结果相较于传统单相流的计算而言更具有实际意义.

（1）基于多孔介质的流体力学理论，忽略化学因素的影响，开发了TOUGH2-FLAC3D气液固三相耦合搭接程序. TOUGH2-FLAC3气液固三相耦合计算程序能够模拟真实的气水二相流，能够研究渗流过程中气水的相互作用及其对渗流过程的影响.

（2）从物理机制上探讨渗流与骨架变形之间的相互影响，考虑骨架变形导致的孔隙率渗透率的变化，在该耦合计算过程中对孔隙率、渗透率和毛管压力进行更新，建立了固液气三相耦合分析平台.

（3）与经典的Liakopoulos排水试验结果进行对比，结果显示本研究程序计算的孔隙水压力水头与实验测点的数据基本吻合. 与Ehlers等^[37]、Nagel 等^[38]的计算结果进行对比，结果表明本研究开发的TOUGH2-FLAC3D气液固耦合模型计算的出口流速在初始阶段与实验测点有所偏差，可能是因为渗流没有达到稳定，在非稳定阶段目前的数值模型不能很好地模拟，但是随着模拟时间的延长，在达到稳定渗流后，本研究的计算程序所获得的结果与实验测点值更为接近. 将模型试验中数值模拟计算与实测结果进行对比，可以看出，气相在渗流过程中的影响不可忽略，考虑气相影响的模拟分析与实际情况更接近.

（4）基于TOUGH2和FLAC3D的两相渗流-应力弱耦合相关计算程序目前尚未公布，但是其应用领域广泛，相关的耦合计算程序一直是众多研究者研究的重点. 本研究详细介绍了软件开发的详细工作流程，为同类型耦合软件的开发搭接提供了一定的借鉴.