引言

深入探索多孔介质中多相流动和固体变形特征对精准认识地下资源开发利用, 例如石油勘探开发、煤层气和页岩气开采、二氧化碳地质封存等领域, 具有重要的工程意义^[1-4]. 在流体输送过程中局部压力和饱和度梯度下的渗流作用改变了孔隙结构特征, 打破了储层内部的水力平衡; 固体空间结构的变化影响了储层内部流体的运移路径和流通能力, 具体表现为多孔介质孔隙度和渗透率等物性参数的动态响应. 发展和建立两相渗流与固体变形耦合模型及其相应的数值方法是准确描述和阐明储层内水力行为的有效手段^[5-7].

针对两相渗流和固体变形耦合模型的数值方法主要包括有限元法^[8-10]、有限体积法^[11-12]以及混合计算方法^[13-14]等. 对于两相流问题的离散和数值求解, 传统有限元能够确保流体域内整体的质量平衡, 但是无法实现流体在跨越单元间的局部质量平衡, 因此不适用于求解对流占主导的非线性不混溶问题; 有限体积法在确保局部流量平衡的同时, 具有较高的计算效率, 因此在工业界和工程尺度模拟中得到了广泛应用. 在通常情况下, 有限体积方法的计算结果仅能满足低阶有效精度; 混合计算方法是指对流体和力学控制方程分别利用不同的数值方法进行求解, 该方法集成了各种数值方法的优势. Rutqivst^[15]在耦联流体软件TOUGH2和固体变形软件FLAC3D的基础上, 开发了在油气、地热、可燃冰开采以及二氧化碳地质封存等领域广泛使用的TOUGH-FLAC耦合模拟器. 在此模拟器中, 负责流体计算的TOUGH2软件采用的是有限体积法, FLAC3D采用有限差分法计算固体变形, 此方法在网格划分和数据传递效率等方面具有一定的局限性.

间断伽辽金有限元(discontinuous Galerkin finite element method, DG-FEM)结合了有限体积与有限元两种方法的优点, 具有高阶精度和局部单元物理守恒性的特点, 容易实现局部网格加密和各单元多项式的单独选取^[16-19]. 近些年来, DG-FEM在流体和弹性力学领域得到了广泛的关注和持续的发展^[20-22]. Klieber和Riviere^[23]在解耦两相渗流微分方程组的基础上, 提出了不可压缩两相流自适应间断有限元方法. Epshteyn和Rivière^[24]提出了多孔介质不可压缩两相流完全隐式方程的不连续伽辽金方法, 该方法在结构化和非结构化网格计算中都具有较好的鲁棒性. 李子然和吴长春^[25]、陈韵骋等^[26]建立了基于线弹性力学问题的局部间断有限元方法. 在流固耦合方面, Liu等^[27-28]提出了连续和间断伽辽金耦合框架, 利用自适应加罚方法提高间断伽辽金有限元方法在模拟大型多孔弹性问题时的适用性和计算效率, 并将此方法推广至多孔介质水力压裂模拟中. 目前关于间断伽辽金有限元方法在计算可压缩两相流固模型方面的研究鲜有报道, 还有待进一步深入研究.

基于上述认识, 本文首先建立可压缩两相流固耦合控制方程, 其中包括考虑毛细压力和重力作用的可压缩两相渗流方程以及线性多孔弹性方程. 在此基础上, 推导出流固耦合方程的间断伽辽金弱形式; 然后, 通过一维Terzaghi流固问题验证模型的准确性; 最后, 分别开展二维平面应变条件下和考虑重力效应的三维流固数值算例, 分析流体压力、饱和度以及固体应力、位移的分布特征, 探讨加罚因子对模拟结果的影响.

1.   两相流固耦合控制方程

1.1   两相渗流方程

假设模型中湿润相和非湿润相为不混溶流体, 两者在多孔介质中的流速满足达西定律. 在考虑毛细压力、重力效应和固体变形的基础上, 等温可压缩两相渗流控制方程可表示为^{[2, 29]}

式中, 下标$\;\beta\; 为\;{\rm{w}}$和$ {\rm{n}}{\rm{w}} $分别表示湿润相和非湿润相, $ {S}_{\beta } $为饱和度, $\; {\rho }_{\beta } $和$ \;{\mu }_{\beta } $分别为流体密度和黏度, $ {C}_{\beta } $为流体压缩系数, $ {p}_{\beta } $为流体压力, $ {\boldsymbol{k}} $为多孔介质绝对渗透率张量, ${k}_{{\rm{r}}\beta }$为相对渗透率, $ {Q}_{\beta } $为源汇项, $ {\varepsilon }_{\rm{v}} $为体应变, $ {p}_{\rm{c}} $和$ \left|{p}_{\rm{c}}\right| $分别为毛细压力和毛细压力对湿相饱和度的导数.

为了求解方程, 需要补充如下辅助公式

式中, $ {p}_{\rm{e}} $为毛细管吸入压力, $ {S}_{\rm{e}} $为有效饱和度, $ \lambda $为孔径分布指数, $ {S}_{{\rm{r}}{\rm{w}}} $和$ {S}_{{\rm{r}}{\rm{n}}{\rm{w}}} $分别为湿润相和非湿润相的残余饱和度.

此外, 选取Genuchten−Mualem相对渗透率模型描述湿润相和非湿润相的流通能力, 表达式如下所示^[30]

模型的初始条件为

模型中Dirichlet和Neumann边界条件如下所示

式中, ${\boldsymbol{\lambda }}_{{\rm{n}}{\rm{w}}}=-{\rho }_{{\rm{n}}{\rm{w}}}\dfrac{{\boldsymbol{k}}{k}_{{\rm{r}}{\rm{n}}{\rm{w}}}}{{\mu }_{{\rm{n}}{\rm{w}}}}$, ${\boldsymbol{\lambda }}_{{\rm{w}}}=-{\rho }_{{\rm{w}}}\dfrac{{\boldsymbol{k}}{k}_{{\rm{r}}{\rm{w}}}}{{\mu }_{{\rm{w}}}}$, $ {\boldsymbol{n}} $为单位法向量, $ {p}_{{\rm{n}}{\rm{w}}0} $和$ {p}_{{\rm{n}}{\rm{w}}}^{\rm{b}} $分别为初始和边界上非湿润相孔隙压力, $ {S}_{{\rm{w}}0} $和$ {S}_{{\rm{w}}}^{\rm{b}} $分别为初始和边界上湿润相饱和度, $ {g}_{{\rm{w}}} $和$ {g}_{{\rm{n}}{\rm{w}}} $分别为边界上湿润相和非湿润相的流量.

1.2   两相渗流方程的间断伽辽金弱形式

图1为模拟几何区域、边界条件和网格示意图. 在计算域$ \varOmega $中, 边界$ \partial \varOmega $由不与内部单元接触的外边界$ \partial {\varOmega }_{\rm{o}} $和内部相邻单元共享的内边界$ \partial {\varOmega }_{\rm{i}} $两部分组成. 由图1可知, 内边界$ \partial {\varOmega }_{\rm{i}} $为相邻单元$ {E}^{+} $和$ {E}^{-} $共同拥有; 外边界$ \partial {\varOmega }_{\rm{o}} $由Dirichlet边界$ \partial {\varOmega }_{{\rm{o}}{\rm{D}}} $和Neumann边界$ \partial {\varOmega }_{{\rm{o}}{\rm{N}}} $两者构成, 即$ \partial {\varOmega }_{\rm{o}}=\partial {\varOmega }_{{\rm{o}}{\rm{D}}}\cup \partial {\varOmega }_{{\rm{o}}{\rm{N}}} $. 设$ {\cal{T}}_{\rm{h}} $为$ \bigcup \left\{E\right\} $, 是区域 $ \varOmega $的正则剖分, $ \partial E=\partial {\varOmega }_{\rm{o}}\cup \partial {\varOmega }_{\rm{i}} $表示为所有边界的集合, 则间断有限元空间定义为

图  1  几何区域、边界条件和网格示意图

Figure  1.  Domain with boundary conditions and mesh skeleton

下载: 全尺寸图片 幻灯片

式中, 离散空间${p}_{{\rm{r}}}\left(E\right)$为单元E上不超过r次的多项式的集合.

为了推导两相渗流方程的间断伽辽金弱形式, 首先在式(1)等号的两侧乘以试函数$ \widetilde {{p}_{{\rm{n}}{\rm{w}}}} $, 并在单元计算域E内进行积分, 可得到

利用分部积分和散度定理, 式(15)中等号左侧第二项可表示为

将式(16)代入式(15), 累加所有单元区域, 可得

在间断伽辽金方法中, 变量在单元内边界$ \partial {\varOmega }_{\rm{i}} $处满足不连续设定, 因此式(17)中等号左侧第三项可以改写为

将式(18)代入式(17)中, 结合稳定项、加罚函数项以及相应边界条件, 可得

式中, $ \left\langle {} \right.{\boldsymbol{X}}\left. {} \right\rangle =\dfrac{1}{2}\left({\boldsymbol{X}}\cdot {{\boldsymbol{n}}}^{+}-{\boldsymbol{X}}\cdot {{\boldsymbol{n}}}^{-}\right) $和$ [\![{\boldsymbol{X}}]\!]={\boldsymbol{X}}\cdot {{\boldsymbol{n}}}^{+}+{\boldsymbol{X}}\cdot {{\boldsymbol{n}}}^{-} $分别为$ {\boldsymbol{X}} $在相邻单元共同边界上的平均和跳跃, $ h $为单元尺寸, $ {\delta }_{\rm{f}} $为加罚因子, $ {\omega }_{\rm{f}} $的选取决定了间断有限元方法的类型. $ {\omega }_{\rm{f}}=1 $, 表示为非对称内罚伽辽金法(non-symmetric interior penalty Galerkin method, 简称NIPG); $ {\omega }_{\rm{f}}=0 $, 表示为非完全内罚伽辽金法(incomplete interior penalty Galerkin method, 简称IIPG); $ {\omega }_{\rm{f}}=-1 $, 表示为对称内罚伽辽金法(symmetric interior penalty Galerkin method, SIPG). 在本文中, 选取$ {\omega }_{\rm{f}}=1 $^[23], 则式(19)为流体方程的非对称内罚伽辽金弱形式.

利用上述类似的推导过程和方式, 可以得到式(2)的间断伽辽金弱形式, 如下所示

式(19)和式(20)为考虑固体变形作用下可压缩两相渗流方程的间断伽辽金有限元方法的弱表达形式.

1.3   固体变形方程及其弱形式

假设多孔介质固体变形满足弹性小变形条件, 那么固体变形的平衡方程、几何方程、本构方程以及边界条件如下所示

式中, 总应力$ {\boldsymbol{\sigma }}={\boldsymbol{\sigma }}{{'}}-\alpha {\boldsymbol{I}}p $, 则式(21)可改写

式中, $ {\boldsymbol{\sigma }}{{'}} $为有效应力张量, $ {\boldsymbol{\varepsilon }} $为应变张量, $ {\boldsymbol{u}} $为位移张量, $ {\boldsymbol{f}} $为体积力张量, $ \alpha $为Biot常数, $ {\boldsymbol{D}} $为4阶弹性张量, $ {\boldsymbol{I}} $为单位张量, 平均孔隙压力$ p={S}_{{\rm{n}}{\rm{w}}}{p}_{{\rm{n}}{\rm{w}}}+{S}_{{\rm{w}}}{p}_{{\rm{w}}} $.

利用类似式(19)和式(20)的推导过程, 结合相应的边界条件, 可以得到固体变形方程的弱形式, 如下所示

式中, $ \widetilde {{\boldsymbol{u}}} $为试函数, $ G $为剪切模量, $ {\delta }_{\rm{s}} $为加罚因子. 本文选取$ {\omega }_{\rm{s}}=0 $^[26], 则式(27)为固体变形方程的非完全内罚伽辽金弱形式.

1.4   数值方法的实现

方程式(19), 式(20)和式(27)是多孔介质中可压缩两相渗流和固体变形的弱形式. 本文利用通用有限元软件COMSOL Multiphysics内置的Weak Form PDE模块求解可压缩两相流固耦合模型. 基于二次不连续拉格朗日形函数构建对解域的间断伽辽金有限元逼近, 选取向后差分公式(backward differentiation formula, BDF)隐式求解器. 其中, BDF求解器采用自适应时间步长算法和可变离散化阶数, 离散化的BDF精度从1 ~ 5不等. 基于LU分解, 采用MUMPS (multifrontal massively parallel sparse)直接求解器对线性系统进行求解, 并选取完全耦合的求解器(fully coupled solver)将线性求解器应用于牛顿单步法的非线性问题. 默认情况下, 初始时间步长设置为结束时间的0.1%, 相对容差为0.01.

2.   模型验证

图2描述了一维Terzaghi固结问题的几何示意图以及相应的边界条件. 为了确保孔隙压力沿上下两端保持为零, 设置模型上下两端为排水边界. 在模型上边界瞬间施加均布载荷$ q=-5.00 $ MPa. 本模型长度$ 2 h=10.00 $ m, 其他参数的具体取值如下^[2]: 剪切模量为30.00 GPa, 压缩模量为50.00 GPa, 多孔介质渗透率为$ {1.00\times 10}^{-12}\;{\rm{m}}^{2} $, 孔隙度为0.25, 流体压缩系数为$ {1.00\times 10}^{-12}\;{{\rm{Pa}}}^{-1} $, 流体黏度为 $ {1.00\times 10}^{-3}\;{\rm{Pa}}\cdot {\rm{s}} $, 加罚因子$ {\delta }_{\rm{f}} $和$ {\delta }_{\rm{s}} $分别设置为1.00和10.00^[31].

图  2  Terzaghi固结问题几何示意图

Figure  2.  Sketch of Terzaghi’s consolidation problem

下载: 全尺寸图片 幻灯片

图3为Terzaghi固结问题的理论解^[32]和DG有限元数值解的比较. 结果显示, DG有限元计算得到的数值解与理论解吻合较好, 进而验证了间断伽辽金方法在求解流固耦合方程时的可行性和有效性.

图  3  Terzaghi固结问题理论解和数值解的比较

Figure  3.  Comparison of the analytical and numerical results for Terzaghi’s consolidation problem

下载: 全尺寸图片 幻灯片

3.   算例分析

3.1   二维平面两相流固算例

本节利用间断伽辽金方法开展在二维平面应变条件下两相渗流和固体变形耦合数值计算. 在本算例中, 分别选取水和氢气作为湿润相和非湿润相. 模型的几何尺度、初始条件、边界条件和监测线如图4所示. 模型的长度和宽度均设置为10.00 m. 多孔介质的弹性模量和泊松比分别为10.00 GPa和0.30; 渗透率和孔隙度分别设定为$ 1.00\times {10}^{-14}\;{\rm{m}}^{2} $和$ 0.20 $. 模拟区域的初始气体压力为4.00 MPa, 初始水饱和度为0.80. 模拟左边界和下边界为位移约束, 在上边界和右边界分别施加$ -1.00\;\times\;{10}^{7} $ MPa (负号表示压应力方向)的垂直和水平应力. 模型左下角为注入边界, 设置注入的氢气压力和水饱和度分别为5.00 MPa和0.80; 右上角为出口边界, 出口气体压力和水饱和度设置为4.00 MPa和0.20.

图  4  模拟几何尺度、初始和边界条件

Figure  4.  Simulation geometrical configuration with boundary and initial conditions

下载: 全尺寸图片 幻灯片

数值模拟分两个阶段开展. 首先, 开展储层力学和流体平衡的稳态计算. 然后, 将第一步计算得到的应力和孔压等水力信息作为含水层氢气注入模拟算例的初始条件. 整个模拟时间为1500 s, 其他变量的具体数值如表1所示.

表  1  两相流固数值模拟的属性参数

Table  1.  Parameters for two-phase flow simulation

Parameter	Value	Unit
elastic modulus, $ E $	$ 30.00 $	$ {\rm{G}}{\rm{Pa}} $
poisson's ratio, $ \nu $	$ 0.25 $	$ - $
permeability, $ k $	$ 1.00\times {10}^{-14} $	$ \;{\rm{m}}^{2} $
porosity, $ \phi $	$ 0.20 $	$ - $
water density, $ \;{\rho }_{{\rm{w}}0} $	$ 1000.00 $	$ {\rm{kg}}/{\rm{m}}^{3} $
water viscosity, $ \;{\mu }_{{\rm{w}}} $	$ 1.00\times {10}^{-3} $	$ {\rm{Pa}}\cdot {\rm{s}} $
water compressibility[32], $ {c}_{{\rm{w}}} $	$ 3.84\times {10}^{-10} $	$ {{\rm{Pa}}}^{-1} $
hydrogen density, $\; {\rho }_{{\rm{n}}{\rm{w}}0} $	$ 3.18 $	$ {\rm{kg}}/{\rm{m}}^{3} $
hydrogen viscosity, $\; {\mu }_{{\rm{n}}{\rm{w}}} $	$ 9.02\times {10}^{-6} $	$ {\rm{Pa}}\cdot {\rm{s}} $
hydrogen compressibility[32], $ {c}_{{\rm{n}}{\rm{w}}} $	$ 7.71\times {10}^{-7} $	$ {{\rm{Pa}}}^{-1} $
entry pressure, $ {p}_{\rm{e}} $	$ 1.00\times {10}^{4} $	$ {\rm{Pa}} $
distribution index, $ \lambda $	$ 0.46 $	$ - $
relative permeability variable, $ \omega $	$ 2.00 $	$ - $
variable, $ {\omega }_{\rm{f}} $	1.00	$ - $
penalty factor[27], $ {\delta }_{\rm{f}} $	1.00	$ - $
variable, $ {\omega }_{\rm{s}} $	0.00	$ - $
penalty factor[27], $ {\delta }_{\rm{s}} $	10.00	$ - $

下载: 导出CSV 

| 显示表格

图5和图6分别给出了注气第1500 s气体压力$ {p}_{{\rm{n}}{\rm{w}}} $、水饱和度$ {S}_{{\rm{w}}} $、水平位移$ u $和水平有效应力$ {\sigma }_{x}^{{{'}}} $的空间分布云图以及变量沿监测线的分布图. 由图5和图6可知, 随着氢气的持续注入, 气体压力影响范围扩大至距注入边界7 m左右, 并引起固体膨胀变形且变形影响范围扩散至边界处. 究其原因是因为在多孔弹性介质中位移或者应变可传播到压力前沿之前. 受毛细压力和气水两相物理属性差异的共同影响, 氢气主要聚集于注入边界附近. 在多孔介质内孔压在远离注入边界的过程中呈现递减趋势, 而有效应力沿监测线增加.

图  5  第1500 s气体压力$ {p}_{{\rm{n}}{\rm{w}}} $、水饱和度$ {S}_{{\rm{w}}} $、水平位移$ u $和水平有效应力$ {\sigma }_{x}^{{{'}}} $的分布云图

Figure  5.  The distribution of gas pressure $ {p}_{{\rm{n}}{\rm{w}}} $, water saturation $ {S}_{{\rm{w}}} $, horizontal displacement $ u $ and horizontal effective stress $ {\sigma }_{x}^{{{'}}} $ at t = 1500 s

下载: 全尺寸图片 幻灯片

图  6  第1500 s时气体压力$ {p}_{{\rm{n}}{\rm{w}}} $、水饱和度$ {S}_{{\rm{w}}} $、水平位移$ u $和水平有效应力$ {\sigma }_{x}^{{{'}}} $沿监测线的分布

Figure  6.  The profiles of gas pressure $ {p}_{{\rm{n}}{\rm{w}}} $, water saturation $ {S}_{{\rm{w}}} $, horizontal displacement $ u $ and horizontal effective stress $ {\sigma }_{x}^{{{'}}} $ along the monitoring line at t = 1500 s

下载: 全尺寸图片 幻灯片

3.2   加罚因子的影响

图7给出了$ {\delta }_{\rm{s}}=0.10 $, 1.00和10.00时水平位移$ u $和水平有效应力$ {\sigma }_{x}^{{{'}}} $沿监测线的分布. 图8给出了$ {\delta }_{\rm{s}}=0.10 $时水平位移$ u $和水平有效应力$ {\sigma }_{x}^{{{'}}} $的分布云图. 图9给出了$ {\delta }_{\rm{f}}=0.01 $, 0.10, 1.00和10.00时气体压力$ {p}_{{\rm{n}}{\rm{w}}} $、水饱和度$ {S}_{{\rm{w}}} $沿监测线的分布. 结果表明, 不同加罚因子条件下, 变量$ u $, $ {\sigma }_{x}^{{{'}}} $, $ {p}_{{\rm{n}}{\rm{w}}} $和$ {S}_{{\rm{w}}} $沿着监测线的变化趋势大体保持一致. 在计算固体变形时, 计算结果对加罚因子的选取更为敏感. 当$ {\delta }_{\rm{s}} $ = 0.10时, 应力和位移均在局部出现明显的波动; 当$ {\delta }_{\rm{f}} $ = 0.01时, 饱和度在小范围内发生轻微浮动. 这主要由于加罚因子$ {\delta }_{\rm{s}} $和$ {\delta }_{\rm{f}} $的降低减少了自变量位移、压力以及饱和度在单元处弱连续性的约束, 造成了变量局部波动的现象. 通过适当提高加罚因子可以加强变量在跨越单元交界面处的连续性, 但是较高的加罚值会引发病态条件矩阵的可能性.

图  7  第1500 s时不同加罚因子$ {\delta }_{\rm{s}} $条件下水平位移$ u $和水平有效应力$ {\sigma }_{x}^{{{'}}} $沿监测线的分布

Figure  7.  The profiles of horizontal displacement $ u $ and horizontal effective stress $ {\sigma }_{x}^{{{'}}} $ along the monitoring line at t = 1500 s with different values of penalty parameter $ {\delta }_{\rm{s}} $

下载: 全尺寸图片 幻灯片

图  8  第1500 s时$ {\delta }_{\rm{s}}=0.01 $水平位移$ u $和水平有效应力$ {\sigma }_{x}^{{{'}}} $分布云图

Figure  8.  The distribution of horizontal displacement $ u $ and horizontal effective stress $ {\sigma }_{x}^{{{'}}} $ with $ {\delta }_{\rm{s}}=0.01 $ at t = 1500 s

下载: 全尺寸图片 幻灯片

图  9  第1500 s时不同加罚因子$ {\delta }_{\rm{f}} $条件下气体压力$ {p}_{{\rm{n}}{\rm{w}}} $、水饱和度$ {S}_{{\rm{w}}} $沿监测线(0, 0.05)−(10, 9.95)的分布

Figure  9.  The profiles of gas pressure $ {p}_{{\rm{n}}{\rm{w}}} $ and water saturation $ {S}_{{\rm{w}}} $ along monitoring line at t = 1500 s with different values of $ {\delta }_{\rm{f}} $

下载: 全尺寸图片 幻灯片

值得注意的是, 受流体边界条件、初始条件以及流体方程的非线性特征等因素的影响, 传统有限元方法在求解本节流体模型时无法收敛, 此情况在一定程度上体现了间断伽辽金有限元在求解可压缩两相流固方程上的优越性. 本节将传统有限元的计算结果作为基准方案, 通过比较间断伽辽金有限元和传统有限元两者计算得到的水平位移和水平有效应力之间的相对误差, 分析加罚因子$ {\delta }_{{\rm{s}}} $对计算结果的影响. 水平位移和水平有效应力的计算误差$ er{r}_{\rm{u}} $和$ er{r}_{\rm{s}} $的表达式为^[33]

式中, $ l $为几何模型尺度, $ {u}_{{\rm{C}}{\rm{G}}} $和$ {\sigma }_{x{\rm{C}}{\rm{G}}}^{{{'}}} $分别为有限元方法计算得到的水平位移和水平有效应力, $ {\Delta }{u}_{{\rm{C}}{\rm{G}}} $和$ {\Delta }{\sigma }_{x{\rm{C}}{\rm{G}}}^{{{'}}} $分别为$ {u}_{{\rm{C}}{\rm{G}}} $和$ {\sigma }_{x{\rm{C}}{\rm{G}}}^{{{'}}} $极差, $ {u}_{{\rm{D}}{\rm{G}}} $和$ {\sigma }_{x{\rm{D}}{\rm{G}}}^{{{'}}} $分别为间断有限元方法计算得到的水平位移和水平有效应力.

图10给出了不同加罚因子$ {\delta }_{\rm{s}} $条件下间断伽辽金有限元和传统有限元方法计算得到的水平位移和水平有效应力的比较. 结果表明, 随着$ {\delta }_{\rm{s}} $的增加, 间断有限元与传统有限元计算得到的水平方向的变形和应力的相对误差随之减少, 这说明了加罚项的增加可以有效抑制有限元函数在跨越单元时的不连续性.

图  10  不同加罚因子$ {\delta }_{\rm{s}} $条件下间断有限元和传统有限元方法计算水平位移和水平有效应力的比较

Figure  10.  Comparison of the horizonal displacement and horizonal effective stress between DG and the CG results with different values of penalty parameter $ {\delta }_{\rm{s}} $

下载: 全尺寸图片 幻灯片

3.3   考虑重力效应的两相流固三维算例

在本节中, 将开展含水层注氢气的三维数值模拟. 在此算例中, 侧重分析流体重力效应对模拟结果的的影响. 如图11所示, 三维算例的几何尺度为100.00 m×2.00 m×10.00 m. 模拟上边界和右边界为应力边界且${t}_{x}={t}_{z}=\;-1.00\;\times\;{10}^{7}\;{\rm{M}}{\rm{Pa}}$, 其余边界均设定为位移约束. 氢气从左侧边界注入, 注入率${g}_{{\rm{nw}}}= $$ 1.00\;\times\;{10}^{-4}\;{\rm{kg}}/({\rm{m}}^{2}\cdot {\rm{s}})$, 右侧为出口边界, 边界上的气体压力和水饱和度分别为4.00 MPa和0.20, 其余边界则为不流动边界. 整个注气过程持续1 d, 其他变量的数值如表1所示.

图  11  三维模型几何尺度以及模拟初始和边界条件

Figure  11.  Geometrical configuration with boundary and initial conditions of 3D case

下载: 全尺寸图片 幻灯片

图12给出了注气1 d后气体压力$ {p}_{{\rm{n}}{\rm{w}}} $、气体饱和度$ {S}_{{\rm{n}}{\rm{w}}} $、垂直位移$ w $和水平有效应力$ {\sigma }_{x}^{{{'}}} $的分布云图. 图13为在考虑和不考虑重力效应两种情况下点a, b和c处气饱和度$ {S}_{{\rm{n}}{\rm{w}}} $随时间变化图. 由图12可知, 氢气注入1 d后, 孔压沿着水平x方向传播至距左边界约70.00 m处, 储层边界处气体饱和度和孔隙压力分别增加至约0.52和6.50 MPa, 孔压的增加引起垂直位移的提高和水平有效应力的降低. 在注入过程中, 由于氢气和水两者密度的差异, 引起气体上浮且聚集于模型顶部, 造成顶部点c处气体饱和度要明显高于底部点a处气体饱和度.

图  12  注气1 d后气体压力$ {p}_{{\rm{n}}{\rm{w}}} $、气体饱和度$ {S}_{{\rm{n}}{\rm{w}}} $、垂直位移$ w $和水平有效应力$ {\sigma }_{x}^{{{'}}} $的分布图

Figure  12.  The distribution of gas pressure $ {p}_{{\rm{n}}{\rm{w}}} $, gas saturation $ {S}_{{\rm{n}}{\rm{w}}} $, vertical displacement $ w $ and horizontal effective stress $ {\sigma }_{x}^{{{'}}} $ after 1 day of injection

下载: 全尺寸图片 幻灯片

图  13  考虑和不考虑重力效应算例中点a, b和c处气饱和度$ {S}_{{\rm{n}}{\rm{w}}} $随时间变化图

Figure  13.  The temporal evolution of gas saturation $ {S}_{{\rm{n}}{\rm{w}}} $ at points a, b and c in the cases with and without gravity

下载: 全尺寸图片 幻灯片

4.   结论

本文建立了考虑毛细压力和重力效应的可压缩两相渗流与孔隙介质变形耦合作用的力学模型. 在此基础上, 分别推导出变形孔隙介质中两相渗流方程的非对称内罚伽辽金弱形式以及固体变形方程的非完全内罚伽辽金弱形式. 主要结论如下:

(1)通过对比一维Terzaghi固结问题的解析解和数值解, 验证了间断伽辽金有限元方法在求解流固耦合问题方面的可行性和准确性;

(2)开展了二维平面应变条件下和考虑重力效应作用的三维两相流固数值算例. 结果显示, 随着气体的注入, 气体饱和度和孔压随之增加, 导致多孔介质膨胀变形和有效应力降低. 由于重力的影响, 气体会上浮聚集于顶部边界;

(3)在二维平面算例中, 讨论了加罚因子对计算结果的影响. 结果表明, 在相同网格条件下, 加罚因子$ {\delta }_{\rm{f}} $和$ {\delta }_{\rm{s}} $的降低会导致流体和固体信息在局部出现不同程度的波动. 固体变形参量对加罚因子$ {\delta }_{\rm{s}} $更为敏感, 提高加罚因子$ {\delta }_{\rm{s}} $可以有效抑制位移和应力在跨越单元时的不连续性.