独立变量的热力学约束：多孔弹塑性断裂相场模型流固耦合项的重构

一、引言

孔隙介质的断裂行为与其饱和状态密切相关：干燥状态下，断裂仅由固体骨架的变形与破坏控制；或部分饱和状态下，断裂由骨架变形与孔隙中流体流动共同驱动。这类流体饱和孔隙介质中的断裂问题具有重要的工程与地质意义，例如油气储层增产改造、增强型地热系统建造、水泥基材料碳化等环境辅助断裂，以及冰隙发育、岩脉侵入等地质过程。

在水力压裂建模方面，早期基础性研究通常假定线弹性以获得解析解，多数经典数值模型也依赖线弹性断裂力学。为克服离散裂纹追踪的局限，相场（梯度损伤）公式被提出并成功应用于延性断裂、疲劳断裂等各类问题，随后被扩展到流固耦合问题：流体质量演化由含Darcy定律的质量平衡方程控制，总能量泛函则通过引入孔隙压力效应进行修正。

现有文献中总能量泛函的具体形式存在差异。主流做法是直接将孔隙压力p作为Helmholtz自由能泛函的独立变量（本文称之为“混合公式"），另一些研究则以流体含量变化ξ为独立变量。从热力学角度看，Helmholtz自由能以应变ε、流体含量ξ和内部损伤变量d为规范独立变量，而Gibbs自由能以应力σ、流体压力p和d为独立变量。混合公式在保持Helmholtz自由能形式的同时用压力p替换流体含量ξ，虽然便于计算（压力直接出现在边界条件中且易于解释），但Helmholtz与Gibbs自由能之间的Legendre变换不再成立，这种不一致可能导致推导广义驱动力时的遗漏，进而损害力学响应的精度。

另一个缺口在于，现有流固耦合相场模型主要关注纯多孔弹性介质中的裂纹萌生与扩展。在高压流体与高围压的现实条件下，塑性变形可能变得显著。对于多孔弹塑性介质中的机械诱导破坏，基于细观力学并引入强度准则的损伤模型能够灵活捕捉断裂萌生与扩展，但塑性变形对水力裂缝扩展的耦合影响仍研究不足。

张丰收团队和芮易团队在“Soft Condensed Matter"期刊发表了文章“Revisiting the hydromechanical formulation of a micromechanics-based phase-field model for poro-elastoplastic media"，本文旨在重新审视流固耦合公式并给出更严格的推导：以应变ε、流体含量ξ和相场变量d为独立变量定义Helmholtz自由能，借助代表性体积单元（RVE）确定微裂纹状态并推导拉压状态下的广义应力，采用前期工作提出的内聚型退化函数描述相场演化，并采用含关联流动法则的Drucker–Prager型强度准则描述塑性流动。

二、模型公式

图1  双尺度研究对象示意图：宏观裂纹用相场标量d∈[0,1]正则化（左），饱和RVE中随机分布币形微裂纹（中、右）

考虑饱和域Ω，边界上分别施加牵引力、位移与压力。从中提取代表性体积单元（RVE）描述细观尺度上随机分布的币形微裂纹及大量充液孔隙（图1）。根据应力状态与流体压力，微裂纹可张开或闭合；宏观尺度上断裂用相场变量d∈[0,1]描述。

在自由能框架下，总应变分解为弹性应变εe与非弹性应变εp（后者由微裂纹张开或闭合微裂纹的摩擦滑移引起）。流体含量ξ参与骨架与孔隙流体的相互作用，引入参考Biot系数α0与Biot模量M0，其演化与相场变量相关：α=α0+(1−α0)[1−g(d)]。系统的总能量由体积能与裂纹表面能组成，其中体积能量密度按微裂纹状态分为两种形式：张开微裂纹（拉伸）取ψopen=½ε:Cdam(d):ε+M(αtr[ε]−ξ)²/2，Cdam(d)=g(d)C为退化弹性张量；闭合微裂纹（压剪）取ψclose=½(ε−εp):C:(ε−εp)+½εp:H(d):εp+M0(α0tr[ε−εp]+tr[εp]−ξ)²/2，H(d)为运动硬化模量。压力由p=∂ψ/∂ξ给出，通过Coleman–Noll过程可得相应的应力-应变关系。由微裂纹开闭过渡处应力与能量的连续性可导出H(d)=g(d)C/[1−g(d)]，并定义统一的广义局部应力sp，以tr[sp]=0与tr[sp]<0分别指示微裂纹张开与闭合。

相场演化采用内聚型相场模型的退化函数与裂纹面密度函数，取线性软化律，内部长度lch随最小主应力在拉伸与压缩特征长度之间平滑过渡，从而将断裂强度作为独立材料参数引入。对张开微裂纹，本文推导的相场驱动力为−½g′(d)[ε:C:ε+(1−α0)²p²/K+2(1−α0)ptr[ε]]；而混合公式（以p为独立变量）遗漏了其中的耦合项2(1−α0)ptr[ε]，该遗漏将导致强度面不连续，是本文重点修正的对象。闭合微裂纹的摩擦滑移采用Drucker–Prager型摩擦准则并配以关联流动法则：fp=‖spdev‖+A·tr[sp]/3≤0，塑性势函数与摩擦准则相同；给定应变增量后通过返回映射算法求解塑性应变。流体流动由含Darcy定律的质量平衡方程控制，断裂区的渗透率按缝宽立方关系增强。

三、强度面分析

相场演化方程在均匀损伤（∇d=0）条件下的临界状态ed=0给出材料的强度面。对张开微裂纹与闭合微裂纹分别导出：Fopen=‖σdev‖²/2µ+(σsph+p)²/K−Gc/lch=0；Fclose=‖σdev‖+A(σsph+p)−√(Gcχ/lch)=0，其中σsph与σdev分别为平均应力与偏应力，χ=A²K+2µ。两式均与相场长度尺度无关，且孔压仅通过球应力部分进入强度面，符合有效应力原理的基本思想。为理论验证，作者从Gibbs自由能出发经Legendre–Fenchel变换重新推导，得到与上述一致的强度面表达式，确认了推导的自洽性。

图2  不同孔压（p0=0、5、10、15 MPa）三轴压缩下灰砂岩强度面与实验结果的对比（实验数据取自Zhu等，2023）

将压剪强度面与灰砂岩三轴压缩试验结果对比（图2）。试验中先施加静水围压，再注水建立低于围压的初始孔压，保持围压与孔压恒定后增加轴向载荷直至破坏。结果显示：给定孔压下峰值轴向应力σ1与围压σ3呈线性关系，且该线性关系随孔压变化而平移；所提出的模型成功捕捉了σ1–σ3的线性关系及强度面随孔压的平移。

图3  强度面对比：(a)α0=0.6时不同孔压（p=0、5、10 MPa）；(b)p=5 MPa时不同Biot系数（α0=0、0.6、1）。实线为一致性公式，虚线为混合公式，黑色点线为一致性公式的拉伸断裂强度面

若按混合公式（以p为独立变量）推导张开微裂纹的强度面，则得到Fopen,mixed=‖σdev‖²/2µ+(σsph+α0p)²/K−Gc/lch−p²(1−α0)²/K。图3a表明：干态（p=0）时一致性公式与混合公式一致；饱和态（p>0）下，混合公式导出的强度面在微裂纹开闭过渡处（tr[sp]=0）出现明显不连续，而一致性公式下孔压仅引起强度面沿球应力轴的整体平移。图3b显示：α0=1时两种公式重合；α0<1时混合公式再次在过渡处产生不连续。此外，一致性公式下拉伸断裂强度面与压剪断裂强度面相切，两种断裂模式之间过渡光滑。

图4  主应力空间中的强度面对比：(a)平面应变双轴加载条件；(b)三轴拉伸条件，不同孔压（p=0、5、10 MPa）。实线为一致性公式，虚线为混合公式

为突出混合公式产生的屈服面不连续，图4在主应力空间中分别给出平面应变双轴加载与三轴拉伸条件下的强度面：混合公式（虚线）在拉压过渡处出现明显的台阶式跳跃，而一致性公式（实线）光滑连续。

图5  考虑四阶张量分解退化与塑性流动法则的强度面对比（VD表示体积与偏量部分分别退化，AFR为关联流动法则，NAFR为非关联流动法则；实线为关联流动法则，虚线为非关联流动法则）

进一步考察个体退化函数与塑性流动法则的影响：采用Mori–Tanaka均匀化方案对刚度张量的体积与偏量部分分别退化（VD），并引入非关联流动法则（膨胀系数Aθ<A）。图5显示：无偏应力（‖σdev‖=0）时各强度面重合；随偏应力增大，是否采用个体退化的强度面相互分离但保持连续直至微裂纹开闭过渡点；在过渡点处，非关联流动法则导出的强度面出现跳跃。由于膨胀系数Aθ的作用，闭合微裂纹区的强度显著降低。为保持强度面连续，后续数值实验均采用关联流动法则（A=Aθ）。

四、数值实现

控制方程的弱形式包括塑性演化、相场演化与质量平衡三个方程，对应位移u、相场d与孔压p三个主变量。塑性变量（塑性应变εp与局部应力sp）采用局部返回映射算法更新，并给出考虑压力影响的试探应力。求解采用混合单块-交错格式：位移-压力耦合系统（p–u）以单块格式直接求解，相场d与p–u之间采用交错迭代；每个时间步内先由当前应力计算特征长度lch，再依次求解相场问题与位移-压力问题，直至三场残差同时满足收敛准则。

五、数值算例

5.1 多孔弹性介质中的流体驱动断裂（KGD基准）

图6  KGD模型的(a)几何与(b)网格示意

采用经典KGD水力压裂平面应变基准检验一致性公式与混合公式的差异。无量纲黏度M=3.8×10⁻⁷，小于临界值Mc=3.4×10⁻³，压裂过程处于韧性主导区。利用对称性将无限平面简化为半无限平面，初始半裂纹位于边界（图6a）；采用40 m×120 m的足够大矩形域、7254个四边形单元、最小单元尺寸h=0.1 m（图6b），边界法向约束、周边压力为零。

图7  KGD水力压裂的(a)注入点压力、(b)注入点缝宽、(c)缝长对比。黑色实线为解析解，橙色与绿色虚线分别为一致性公式与混合公式的模拟结果

与解析解对比（图7）：相比混合公式，一致性公式在注入点给出更低的起裂压力，与解析解更接近，表明压力-变形耦合项有效防止了压裂过程中流体压力的高估；注入点缝宽与缝长的对比同样显示一致性公式优于混合公式。

5.2 多孔弹塑性介质中的流体驱动断裂

图8  水力压裂模型的(a)几何与边界条件；(b)白色虚线框区域的局部放大（用于展示多场演化）

在80 m×40 m的二维域中比较多孔弹性与多孔弹塑性模型，初始裂纹位于域中心（图8），采用31270个四边形单元、最小单元尺寸h=0.1 m。图9给出t=5、10、15、25 s时相场d、等效塑性应变εp,eq、局部应力迹tr[sp]与孔压p的演化：裂纹邻近区域局部应力保持受压（tr[sp]<0），塑性变形在裂纹两侧有限宽度内扩展，形成环绕裂纹的鞘状塑性区；裂纹附近局部应力受拉，证实断裂为流体压力驱动的拉伸主导模式。

图9  弹塑性模型在t=5、10、15、25 s的水力压裂模拟结果：(a)相场d；(b)等效塑性应变εp,eq；(c)局部应力迹tr[sp]；(d)孔压p

图10  多孔弹性与多孔弹塑性模型对比：(a)t=25.0 s时的缝宽分布；(b)注入点压力演化

对比缝宽分布（图10a）：塑性变形限制裂缝张开，使缝宽小于弹性情形，与已有文献报道一致。注入压力演化（图10b）凸显了传播过程中的额外塑性耗散：起裂时裂纹附近的塑性变形提供额外阻力，导致破裂压力更高；传播开始后，裂缝内需要维持更高的孔压以支撑这一额外塑性耗散。可见断裂过程伴随显著的塑性能量耗散，其流固耦合行为与纯多孔弹性预测存在实质性差异。

5.3 准脆性多孔介质中的压剪断裂（双轴压缩）

图11  双轴压缩模型示意（含孔试样，平面应变）

最后模拟含孔试样的双轴压缩试验（平面应变，35460个四边形单元、最小单元尺寸h=0.2 mm，图11）：加载分两阶段，先施加5 MPa围压，再在顶部施加位移加载。图12显示断裂过程：uy=0.105 mm时试样整体处于压剪状态（tr[sp]<0），等效塑性应变在孔口萌生并沿倾斜条带局部化形成剪切带；随塑性应变与损伤演化，断裂在孔口形核；继续加载形成典型的剪切（II型）断裂，窄剪切带与损伤区重合。

图12  初始孔压p0=1 MPa时含孔试样的断裂过程：(a)断裂前的局部应力迹tr[sp]；(b)断裂前后的等效塑性应变εp,eq；(c)断裂前后的相场

图13  双轴压缩过程中的孔压演化（单位：MPa）

图13给出孔压演化：加载初期（uy=0.005 mm）压缩使试样孔压升高；随载荷增加，超压逐渐耗散；断裂萌生前（uy=0.105 mm），塑性摩擦滑移引起的体积膨胀使孔压降至初始值以下，局部降压吸引周围流体流入损伤区；断裂后（uy=0.11 mm），剪切带内渗透率急剧增大，沿断裂形成高导流通道并出现负孔压，加速基体流体排出；继续加载至uy=0.30 mm，基体孔压恢复初始水平，达到稳态。

图14  断裂后特征对比：(a)差分位移-差分力曲线；(b)断裂角（平面应变假设下取单位厚度1 m）

图14a为不同初始孔压（p0=0、1、2、3 MPa）下的差分力-位移曲线：加载早期呈线性弹性段；随损伤与塑性累积进入硬化段直至峰值；初始孔压越低，试样承受的峰值载荷越高，与Zhu等的实验结果一致；峰后残余强度由摩擦系数与围压决定，孔压削弱残余强度，且孔压增量与残余强度下降量近似呈线性关系。图14b显示断裂角随初始孔压增大而略微减小，表明孔压可促进纯剪切断裂向混合断裂转变。

六、结论

本文重新审视了多孔弹塑性介质中细观相场断裂模型的流固耦合公式，给出了相场驱动力的自洽表达式。分析表明：非关联流动法则会在强度面中引入跳跃；以流体压力p代替流体含量ξ作为独立变量的混合公式会遗漏相场驱动力中的耦合项2(1−α0)ptr[ε]，同样破坏强度面的连续性。采用关联Drucker–Prager流动法则并补全该耦合项后，强度面在拉压过渡处保持连续，断裂驱动力精确。

所提出的细观流固耦合框架与前期提出的应力依赖内聚退化函数兼容，并能与流固耦合断裂的实验结果良好吻合。多孔弹性KGD基准表明，一致性公式提高了流体驱动断裂模拟的精度。多孔弹塑性介质中的水力压裂与双轴压缩模拟证明，该模型能够同时复现机械扰动诱导的剪切主导断裂与饱和多孔介质中流体注入驱动的拉伸主导断裂，在地热能生产与地下工程领域具有广阔的应用前景。