双连续介质相场模型刻画水力压裂黏滞主导区与流体滞后研究综述

双连续介质相场模型刻画水力压裂黏滞主导区与流体滞后研究综述

一、引言

水力压裂是多孔介质中受流体驱动的断裂问题，涉及固体骨架变形、裂隙—孔隙多尺度流体流动与裂纹扩展的强耦合。相场断裂模型将尖锐裂纹正则化为弥散相场带，能够统一处理裂纹分支、汇合与萌生的复杂几何，已成为水力压裂数值模拟的重要工具。然而，现有相场模型通常采用单一压力场描述整个区域，裂纹内Reynolds平板流与基质Darcy流的过渡依赖相场值或水平集函数插值，导致裂缝—基质界面处的压力被人工抹平。

图1  水力压裂的三种极限状态（ 韧性K、黏滞 M、早期 O）及水力裂缝的压力剖面、相场压力抹平问题

在韧性主导的K顶点，裂尖压力降很小，抹平压力场尚可接受；但在高注入速率、高流体黏度或低断裂韧性的黏滞主导区（M顶点），裂尖附近压力剧烈下降，流体锋面甚至可能落后于裂纹，形成所谓的流体滞后（fluid lag）现象（图1a）。连续抹平的压力场无法刻画这种不连续行为，反而会在裂尖注入点周围产生过高的局部驱动力，使损伤沿横向扩展而非沿裂纹向前推进，诱发数值失稳。因此，发展能够刻画黏滞主导区与流体滞后的相场模型是该领域的关键缺口。

Tao You、 Keita Yoshioka在“Fluid Dynamics”期刊发表了文章“A dual–continuum phase-field model for hydraulic fracturing: Viscosity-dominated regime and fluid lag”，该文提出一种基于双孔介质微观孔隙力学的双连续介质相场模型。核心思想是在宏观域中引入两个独立的流体压力变量：描述宏观裂隙网络的裂纹压力pc，以及描述基质微孔隙的孔隙压力pp。通过微观孔隙力学推导损伤依赖的有效弹性张量、Biot系数张量与Biot模量，实现裂纹压力与孔隙压力在物理上的解耦；再采用固定应力分裂算法处理两场流固耦合，并通过动态域更新与变分不等式约束刻画流体滞后。数值验证覆盖了K、M、O三个经典顶点，展示了模型对真实黏滞主导水力压裂的预测能力。

二、双孔介质相场断裂理论框架

2.1 双尺度状态方程

模型将多孔介质视为两个相互嵌套的特征尺度（图2）。在介观尺度，代表性体积单元RVEc由含流体的币状裂纹与多孔基质组成；在更细的微观尺度，RVEp由固体骨架与球形微孔组成。裂纹体积分数φc与微孔体积分数φp满足φm+φc=1与φs+φp=1。两个尺度上的流体压力分别记为pc和pp，二者不一定处于局部平衡，其差值由基质渗透率和跨尺度质量交换驱动。

图2  含流体裂纹多孔介质的双尺度分离示意图：宏观相场裂纹、介观RVEc与微观RVEp

宏观本构关系由微观孔隙力学严格导出。总应力σ=Ceff:ε(u)−αc pc−αp pp，其中有效刚度Ceff和Biot系数张量αc、αp均通过平均应变集中张量与体积分数表达。当裂纹压力与孔隙压力相等时，该关系退化为常见的单连续介质形式。不同之处在于，模型允许pc≠pp，从而在宏观上刻画裂尖附近的压力跃迁，这正是黏滞主导区的重要特征。

2.2 相场断裂变分框架

模型沿用Francfort和Marigo的变分断裂框架，将总能量写为弹性能、断裂表面能与流体自由能之和，并通过标准的相场正则化引入长度尺度参数ℓ。为了避免裂纹闭合受压时不合理地产生损伤，采用无拉应变能分解，使Ceff在损伤演化中可能呈现各向异性，对应的αc和αp也为张量形式。

流体自由能ψf显式依赖于两个压力与损伤相关的Biot模量，包含pc²/(2Ncc)、pp²/(2Npp)与pcpp/Npc三项耦合。由变分可得相场演化方程，其驱动力D不仅包含正的弹性能密度，还包括两个压力场对损伤相关Biot模量的贡献。这一结构在物理上意味着：裂纹内流体压力与基质孔隙压力均可促进相场损伤演化，但二者的耦合强度由微观结构决定，而非任意标定。

2.3 双尺度质量守恒

裂纹与微孔隙中的质量守恒分别由两条Darcy型方程描述。裂纹渗透率Kc根据裂隙开度w与单元尺寸h按平板流公式w³/(12h)放大，裂隙开度由最大主应变与单元尺寸估计；微孔渗透率Kp为基质渗透率。两个孔隙网络之间通过Kazemi形状因子描述质量交换，满足总体质量守恒。与单连续介质模型相比，双连续结构使裂纹内的快速流动与基质中的缓慢扩散解耦，为后续的流体滞后描述奠定了方程基础。

三、数值解法

3.1 双压力系统的固定应力分裂

四场耦合（位移u、裂纹压力pc、孔隙压力pp、相场v）采用交错求解：在每个时间步内先固定v，耦合求解u—pc—pp，再更新v。为维持两场流固耦合的数值稳定性，将针对单孔问题设计的固定应力分裂推广到双孔系统。分别在介观裂纹尺度和微观孔隙尺度上固定体积应力，得到关于应变—压力耦合项的简化方程，进而将αc:ε与αp:ε用上一迭代步的应力与压力增量表达。

固定应力分裂生成的等效Biot系数与模量组合，使流动子问题在每个交错迭代中仅依赖上一迭代步的应变信息，从而将强耦合问题拆分为更易求解的流动步与力学步。该分裂在形式上无条件稳定，并能在迭代收敛时恢复原始系统的解。

3.2 裂纹流动的动态域更新

裂纹压力pc的物理意义仅在已开裂区域成立；若在整个计算域上求解pc，会在未损伤区产生虚假的裂纹流动。为此引入相场相关的示性函数χ(v)：当v≤vcr时该区域为开裂区，χ=1；否则χ=0，并强制pc=pp。这样裂纹流动方程只在开裂区内非平凡，同时保留了未开裂区单连续介质行为。该动态域更新不需要任何几何预设的裂纹路径，随相场演化自适应调整。

3.3 流体滞后的变分不等式处理

流体滞后区是裂纹已张开但尚未被压裂液填充的近尖区域，其中流体压力应满足非负约束（忽略蒸气压力时pc=pp=0）。将该约束表达为凸锥上的变分不等式：在所有非负试探压力中，寻找使残差不等式成立的压力场。其KKT条件自然将区域划分为压力大于零的流体填充区与压力等于零的滞后区。数值实现借助PETSc的约化空间积极集Newton求解器，无需显式追踪流体锋面位置即可刻画滞后边界。

四、数值验证

验证算例采用平面应变KGD问题，模型尺寸450m×600m，注入点位于左侧初始裂纹，上下边界施加围压σN。网格沿裂纹扩展路径加密，单元尺寸h=0.5m，相场长度尺度ℓ=1.0m。通过调节断裂韧性Gc与流体黏度μ，分别落入韧性主导区（Km>4.13）、黏滞主导区（Km<0.7）及流体滞后区（Tm≈0）。

图3  平面应变KGD问题的几何边界条件与初始相场设置

4.1 韧性主导区验证

取Km=142.5，流体黏度极低，压力沿裂缝近似均匀，裂纹扩展主要由断裂韧性控制。模型给出的注入点净压力、裂纹开度、裂纹半长及t=100s时的压力剖面均与Garagash等解析解吻合，说明双连续介质模型能够退化为标准单连续相场水力压裂模型。

图4  韧性主导区下注入压力、开度、裂纹长度及压力剖面的解析—数值对比

4.2 黏滞主导区验证

降低断裂韧性并提高流体黏度，使Km=0.593进入黏滞主导区。此时裂尖附近出现显著压力梯度，注入压力、裂纹开度与裂纹长度的演化均与解析解一致（图5）。由于能量耗散主要由黏性耗散主导，结果对相场裂尖正则化不如K顶点敏感。

图5  黏滞主导区下注入压力、开度、裂纹长度及压力剖面的解析—数值对比

进一步比较不同时间沿裂纹的开口剖面（图6），数值结果与解析解随时间增长的演化趋势一致，验证了模型对黏滞主导区流体流动与断裂扩展耦合行为的捕捉能力。

图6  黏滞主导区中20–200s裂纹开口剖面的解析—数值对比

4.3 含流体滞后的裂纹扩展

取消围压使Tm≈0，选取四个Gc值对应Km=0.593、0.496、0.405、0.316。以Case 1为例，随时间推进裂纹开度在注入点最大并向单调减小，裂纹半长由约15m增长至约110m；压力剖面明显分为注入段的高压区与近尖的零压滞后区，比例保持稳定，符合自相似解特征。相场分布显示弥散裂纹带被良好解析，压力场仅在开裂区内非零，验证了动态域更新与变分不等式的联合作用。

图7  Case 1裂纹开口与压力剖面随时间演化，及t=140s时的相场与压力场

将四个工况的无量纲开口Ωm与压力Πm与Garagash解析族对比（图8）。随着Km减小，裂尖更易超前于流体锋面，滞后区增大，数值解落在ξf=0.3–0.6的解析族之间，趋势正确。模型略微高估流体锋面推进距离，即低估滞后区大小，可能源于质量守恒在滞后边界处的非严格约束以及数值扩散。

图8  四个工况下无量纲开口（左列）与无量纲压力（右列）的解析—数值对比

五、结论与展望

该文提出了一种变分一致的双连续介质相场模型，用于模拟双孔介质中的水力压裂。通过微观孔隙力学导出的损伤依赖有效弹性张量、Biot系数与Biot模量，使裂纹压力与孔隙压力在物理上独立；固定应力分裂保证了双压力系统流固耦合的数值稳定性；动态域更新与变分不等式约束共同实现了对裂纹流动区域和流体滞后边界的隐式刻画。KGD算例在K、M、O三个顶点的验证表明，模型能够同时复现韧性主导区行为、黏滞主导区压力梯度及含流体滞后裂纹扩展的关键特征。

未来工作可从两个方向深化：一是引入更严格的空化检测方法或近裂尖自适应加密，以减小流体滞后区尺寸的数值偏差；二是将模型推广到基质—裂纹流体交换更活跃的实际多孔岩石，并考虑多相流体情形，如超临界CO₂压裂或岩浆驱动裂纹中的气相滞后，从而拓展该框架在能源地质与地球物理中的应用范围。