一种用于裂缝储层中两相流体流动水力力学建模的瞬态双孔隙度方法

《Finite Elements in Analysis and Design》：A transient dual-porosity approach for hydro-mechanical modeling of two-phase fluid flow in fractured reservoirs

　　•针对可变形裂隙多孔介质中的瞬态两相流，研究了一种水力-机械双孔隙模型。•通过耦合的水力-机械控制方程，推导出了随时间变化的润湿态与非润湿态形状因子。•所提出的瞬态形状因子能够考虑压力、饱和度、毛细效应以及介质变形等随时间变化的因素。•为提升双孔隙模型的精度同时保持计算效率，研究了一种两步法算法。•实验表明，所提出的润湿态与非润湿态形状因子能够准确描述瞬态水力-机械过程。

引言
在工程领域中，有许多问题需要通过对非饱和多孔介质进行水力-机械分析来解决，比如土体固结、地热能开发以及地下氢储存与封存等。这类问题涉及由固体骨架和孔隙空间中的两种不相溶液体组成的多相系统。在许多情况下，介质由裂隙网络构成，这为模拟工作带来了更多复杂性。目前已有多种数值方法用于模拟裂隙多孔介质[1,2]，其中最常用的包括离散裂隙法[3,4]、扩展有限元法[[5], [6], [7]]、相场法[[8], [9], [10]]、基于周期动力学的方法[[11], [12], [13]]、等几何框架[14,15]、多尺度方法[[16], [17], [18]]以及等效方法[19,20]。虽然前五种方法通常被称为直接模型，由于能够明确表示裂隙结构，其精度较高，但计算成本也相当高[21,22]；而等效方法则可通过隐式表达方式降低计算成本，不过所得结果为近似值[[23], [24], [25], [26]]。双孔隙模型是一种高效且实用的等效方法，适用于含有大量裂隙的介质建模，常被用于模拟裂隙多孔介质中单相及多相流体的质量传递过程[2,24,27,28]。

双孔隙模型将裂隙介质视为具有不同特性的两个相互叠加的连续体——裂隙相与基质相，这一概念最初由Barenblatt等人提出[29]。基质相与裂隙相之间的耦合通过“基质相到裂隙相的传递函数”来描述，该函数的计算依赖于“形状因子”参数，而形状因子正是双孔隙模型的核心。双孔隙模型结果的精度在很大程度上取决于形状因子的取值，因此已有不少研究致力于该参数的计算方法。Warren和Root[30]首次在假稳态假设下建立了单相流体中形状因子的计算公式，同时还提出了将裂隙理想化为平行正交网络的“方块糖”模型，该模型依据基质块的尺寸给出恒定数值。Zimmerman等人[31]则在不同的边界条件下求解了非线性常微分方程，他们不仅提出了单相流中裂隙相与基质相之间的非线性相互作用模型，还指出了基质相与裂隙相之间传递过程的瞬态特性。Zimmerman等人的早期研究为后续在非假稳态条件下评估形状因子奠定了基础。Lim和Aziz[32]在没有假稳态假设的情况下求解了压力扩散方程，得出了与基质相-裂隙相传递形状因子相关的结果。Hassanzadeh和Darvish[33]进一步研究了单相形状因子随时间的变化规律，他们发现边界条件会影响到单相流体流动的瞬态形状因子，且在非稳态阶段，恒定的形状因子可能会引入较大误差。Wang等人[34]研究了固体相对形状因子的影响，提出了考虑基质应力敏感性影响的瞬态水力形状因子公式。Rostami等人[35]则利用有限元方法对规则及不规则几何形状的介质的瞬态水力形状因子进行了研究。Khoei和Taghvaei[36]在单相流体流动中建立了瞬态水力-机械形状因子模型，他们的研究还将讨论范围拓展到了水力-机械单相流，并指出除了早期研究中已发现的恒定形状因子的局限性外，忽略固体变形的影响也会带来较大误差。Zimmerman和Paluszny[37]则对裂隙岩体流动与双孔隙概念进行了全面综述。

当存在气体和水这两种不相溶液体时，无论是单孔隙系统还是双孔隙系统，流动动力学都会变得更加复杂[38,39]。这些流体的分布与运动受多种物理过程的共同作用影响。已有大量研究致力于模拟单孔隙系统中的不相溶两相流体流动[[40], [41], [42], [43], [44]]。Kazemi等人[45]将两相流模型拓展到了双孔隙介质中，他们所开发的基于有限差分的模拟器能够考虑流体相对流动性、重力作用、毛细压力以及储层属性的变化等因素。Coats[46]和Ueda等人[47]在双孔隙储层模拟器中采用了Kazemi等人[45]提出的形状因子，尽管他们认为这些形状因子会导致模拟结果偏低，但它们仍在裂隙多孔介质流体流动的数值模拟中得到广泛使用[48,49]。Thomas等人[50]提出了基于有限差分框架的三维三相流体流动建模方法，指出在多相系统中，从裂隙储层中采油的效果取决于基质块的大小与性质，以及裂隙系统的压力和饱和度历史。此外，他们还发现，在稳态状态下，恒定形状因子会受到周围流体类型的影响。Kazemi等人[51]利用经验传递函数，为裂隙介质中的Buckley-Leverett流动提出了解析解与数值解，他们的研究表明，采用指数函数估算基质相与裂隙相之间的传递速率，数值方法的计算效率非常高。Heinemann和Mittermeir[52]验证了Kazemi等人[51]提出的形状因子在任意矩阵几何结构、多相可压缩流体以及各向异性渗透率条件下的适用性。Zimmerman等人[53]将非线性裂隙/基质相互作用模型拓展到了集总参数双孔隙框架下的非饱和水-气两相流动中，他们的研究结果表明，所提出的非线性传递关系能够在流动的早期和后期阶段都准确描述裂隙相与基质相之间的流动情况。Penuela等人[54]提出了一个流动修正因子，用以体现单相及两相流体流动中形状因子的时间依赖性，他们发现，在单相流动条件下，随时间变化的形状因子最终会趋近于Lim和Aziz[32]所给出的数值，而且两相流动的形状因子达到稳态值的速度比单相流动的要慢。Di Donato等人[55]利用解析传递函数研究了位移对恒定两相形状因子的影响，这些函数适用于气体重力排水和逆流吸水过程。Zheng等人[56]则提出了一种新的关系式，通过考虑位移效应来模拟裂隙储层中的不可压缩两相流体流动，与以往模型不同，该模型不会在不可压缩条件下预测出为零的基质相与裂隙相之间的交换率，而是通过采用随局部速度方向变化的恒定形状因子提升了模拟精度。

在大多数两相流研究中，双孔隙模型通常采用基于几何结构的恒定形状因子，该因子在稳态和瞬态条件下均保持不变[57,58]。但在流动变化剧烈的瞬态阶段，由于孔隙间流动具有很强的时间依赖性，这种假设可能会引入显著误差。即便在针对瞬态流动的详细分析中，形状因子也往往被视为与时间无关的，只有少数两相流研究考虑了随时间变化的基质相-裂隙相传递函数[[59], [60], [61], [62], [63]]。Sarma和Aziz[64]为模拟天然裂隙储层中的多相流体流动，提出了一个瞬态水力形状因子，他们使用了方形和非方形的单一基质块，并将结果与高分辨率单孔隙模拟结果进行了对比，证明该瞬态水力形状因子能够准确描述压力扩散和吸水过程。Mei等人[65]针对经历逆流吸水的裂隙状储层，提出了一个随时间变化的形状因子，他们的研究表明，基质的弯曲程度和特征长度会影响流动的持续时间。Li等人[66]为处于逆流吸水状态、具有各向异性基质块的裂隙储层提出了随时间变化的形状因子，并通过高分辨率单基质块模拟和大规模双孔隙模拟验证了该模型的有效性，他们的模型通过将各向异性问题转化为等效的各向同性问题，能够处理排列有序和混乱的渗透率张量。总体而言，在上述大多数关于两相流的研究中，要么假设形状因子始终保持恒定，要么没有考虑固体变形对瞬态形状因子的影响。此外，通常情况下，非润湿相的形状因子会取为润湿相形状因子的相同值或其倍数[[50], [51], [52], [53], [54],61]]。

本研究旨在为裂隙多孔介质中的两相流体流动的水力-机械分析构建一种计算型双孔隙模型。为了提升双孔隙模型中基质相与裂隙相之间质量传递模拟的精度，本研究从耦合的水力-机械控制方程中得到了随时间变化的润湿态与非润湿态形状因子，与传统恒定形状因子模型不同，这些形状因子能够考虑流体压力、饱和度、毛细效应以及固体变形等随时间的变化情况。研究提出了一种两步计算算法：首先从原始模型中选取一个代表性单元格，并将其置于一致的边界条件之下，以此计算瞬态两相形状因子；在第二阶段，将计算得到的随时间变化的形状因子引入双孔隙模型中，用于模拟大尺度的裂隙介质。通过引入符合物理规律的瞬态形状因子以及水力-机械耦合效应，该算法显著提升了双孔隙模型在瞬态阶段的模拟精度，同时仍保持了传统双孔隙方法的计算效率。本研究中的数值模拟主要针对一维和二维问题，但该框架并不局限于此，原则上可以通过采用相应的三维有限元离散化方法扩展到三维问题，这样的扩展或许能为未来在更复杂的实际应用中评估该框架的性能提供有益思路。

章节节选
可变形多孔介质中两相流的控制方程
为模拟可变形裂隙介质中的两相流体流动，需要考虑三种基本成分：可变形的固体相，以及两种不相溶的润湿态与非润湿态流体。根据Biot的多孔介质理论，该系统的控制方程包括应力平衡方程以及每种流体相的质量守恒（或连续性）方程[[67], [68], [69]]。在本研究中，所选取的主要变量包括位移u、润湿态……

水力-机械分析中的瞬态两相形状因子
为建立可变形多孔介质中两相流润湿态与非润湿态传递形状因子的通用关系，本研究在问题中假设存在一个基质块，在该基质块中，润湿态与非润湿态流体的基质压力遵循以下方程：
1Qwwmp˙wm+1Qwnmp˙nm+αSwm?.u˙m+?.[kwm(??pwm+ρw(b?u¨))]=0
1Qnwmp˙wm+1Qnnmp˙nm+αSnm?.u˙m+?.[knm(??pnm+ρn(b?u¨))]=0
此外，应用多相流的广义达西速度定律后，……

有限元离散化
水力-机械两相流的控制偏微分方程可通过有限元Galerkin方法进行离散化，控制方程的弱形式可通过常规的Galerkin方法得出。其中的未知量可以用全局形函数表示为节点值的形式：
asu=NuUpwm=Nppwm,pnm=Nppnmpwf=Nppwf,pnf=Nppnf
据此，控制方程（1）、（3）、（4）可以如下离散化：
ψu=MuU¨+KU?Qwmpwm?Qnmpnm?Fum=0
ψpwm=MwmU¨+Qw

数值模拟结果
为验证所提出的两步计算算法中随时间变化的润湿态与非润湿态形状因子在模拟可变形裂隙多孔介质中瞬态水力-机械两相流时的准确性与性能，本研究通过数值方法解决了多个算例。在第一个算例中，选择了一组二维单一基质块作为基准测试，涵盖了为一维和二维流动场景，以此评估该算法的准确性……

结论
本文提出了一种水力-机械双孔隙模型，用于模拟可变形裂隙多孔介质中的两相流体流动。该模型基于耦合的水力-机械控制方程，推导出了一组随时间变化的润湿态与非润湿态形状因子，这些形状因子能够充分考虑压力、饱和度、毛细效应以及固体变形等参数的随时间变化特性。敏感性分析表明，各种物理参数都会对结果产生影响。具体而言，Mahtab Taghvaei在写作——原始稿件、验证、软件、方法论、研究设计以及形式分析方面做出了贡献；Mohammad Khorshidi则在写作——原始稿件、可视化、验证、软件、方法论、研究设计以及形式分析方面发挥作用；Amir R. Khoei则负责写作——审稿与编辑、验证、监督、项目管理、方法论、研究设计以及概念构建工作。利益冲突声明：作者们声明自己没有已知的可能影响本文所述研究的财务利益或个人关系。Mahtab Taghvaei|Mohammad Khorshidi