最近一个压裂模拟项目把我们逼回了 Matlab。团队里讨论了很久最终决定不依赖商业软件而是自己搭一套“全耦合 XFEM 水力压裂与粘聚单元法”的数值框架。说实话刚开始我对这两个词能不能拼在一起是持怀疑态度的因为 XFEM 擅长描述裂缝几何粘聚单元又自带一套损伤本构两者叠加之后自由度类型、积分策略和收敛控制都会变得非常复杂。这篇文章不是理论综述而是一份折腾记录从为什么常规有限元跑不动压裂到我在 Matlab 里如何组织自由度、装配矩阵、推进耦合迭代再到踩过的一些坑和解决办法。它适合正在做计算力学、数值岩石力学或者想把论文里的算法搬进代码的同行尤其是准备用 Matlab 手写一个中小规模压裂模拟器的人。1. 为什么常规有限元跑不动水力压裂网格依赖与双场耦合1.1 每条裂缝都要重新画网格的噩梦常规有限元的裂缝模型通常把裂缝路径限制在网格边界上。你预先知道裂缝位置就把界面单元或者接触面放在那里如果裂缝要扩展就需要每一步根据应力场重新生成由裂缝切割的网格。这种局部重构的工作量远超出一般人预期每步都要做几何布尔运算把裂缝与新网格交叉旧网格的状态变量需要映射到新网格映射本身引入误差接触面或界面单元的新节点编号、单元拓扑、面绑定全都得重来一旦三维问题局部重剖分网格的鲁棒性会成为瓶颈。XFEM 的思路是用富集函数把几何不连续嵌入到插值形函数里。裂缝穿过一个单元不再需要网格跟随你可以事先就决定哪些节点带额外的自由度。正是这个特性让我下定决心用 XFEM 做裂缝几何的载体。否则每推进一个裂尖都要重新划分网格中小规模的 Matlab 代码根本撑不住一个完整的注入过程。1.2 全耦合不是顺序求解水力压裂有三个过程同时推进岩石变形、基质孔隙渗流、裂缝内流体流动。三者通过有效应力、裂缝开度和流体压力互相咬合。如果采用顺序耦合即每个时间步先更新压力再更新位移然后再更新裂缝在储层渗透率较高或注入排量较大时往往需要非常小的时间步才能稳定而且收敛后可能会出现压力场滞后于裂缝开度扩展的情况。我最初偷懒把解耦当成两个互相交替的子问题结果发现在裂尖开始软化之后交替迭代出现典型的振荡裂缝张开突然增大压力猛地下降然后下一轮又挤回去。网格稍微粗一点就直接发散。换成全耦合 Newton-Raphson 之后矩阵里同时包含位移对压力的导数和压力对位移的导数收敛行为才正常。这也是标题里“全耦合”三个字的现实含义不是宣传话术而是求解器设计的基础容不得一点含糊。2. 方案分工XFEM 管裂缝几何粘聚单元管断裂过程2.1 XFEM 富集的核心位移和孔压都不连续在裂缝穿过一个单元后位移场本身是不连续的。物理上裂缝张开两侧可以独立位移。我在 Matlab 中采用的位移近似是主流的 Heaviside 富集加裂尖渐近富集[ u^h(x) \sum_i N_i(x) u_i \sum_{i \in S_H} N_i(x) [H(x) - H(x_i)] a_i \sum_{i \in S_T} N_i(x) \sum_l F_l(r,\theta) b_i^l ]其中 (H(x)) 是符号距离函数的 Heaviside 形式在裂缝一侧为 1另一侧为 -1。写成 (H(x) - H(x_i)) 而不是直接写 (H(x))是为了保证形函数在节点上的插值性质不被打乱这是 XFEM 最基础的“shift”技巧。(F_l) 是裂尖函数族在二维各向同性材料中常取[ F \left{ \sqrt{r}\sin\frac{\theta}{2},\ \sqrt{r}\cos\frac{\theta}{2},\ \sqrt{r}\sin\frac{\theta}{2}\sin\theta,\ \sqrt{r}\cos\frac{\theta}{2}\sin\theta \right} ]这些渐近函数用来表达裂尖附近应力场的平方根奇异性。实际编码时裂尖函数并不是每个时刻都要全部启用只有当单元包含裂尖且裂尖到单元边界的距离在一定范围内才需要加入这组富集。比位移富集更麻烦的是压力场。裂缝内流体压力与基质的孔隙压力在裂缝壁面是两套量压力穿过裂缝会发生跳跃。如果只用一个连续压力场裂缝内外的流体压力就被强行绑定成连续场注入引起的压差根本表达不出来排量和滤失的计算都不准确。因此孔压近似也要做富集[ p^h(x) \sum_i N_i(x) p_i \sum_{i \in S_P} N_i(x) [H(x) - H(x_i)] c_i ]如果不做这一步从宏观上看压力响应也还算正常但裂缝壁面的滤失和裂缝内压力分布会明显失真尤其是高排量注入场景。这个细节是我一开始没想到的后来翻了几篇经典实现才发现压力场富集几乎就是水力压裂 XFEM 的标配。2.2 粘聚单元解决的是裂尖奇异性问题传统断裂力学用应力强度因子 (K_I) 判断裂纹是否扩展。但水力压裂的裂尖附近不是纯粹理想裂纹裂缝尖端前方存在断裂过程区牵引力随张开先升后降。用内聚力模型可以把起裂、扩展和软化统一进本构关系不需要单独计算 (K_I)这一点在含流体压力的裂缝中尤其舒服。我采用的双线性牵引-分离关系很简单[ t_n K_n w \quad (w \le w_0) ] [ t_n t_{max} \left(1 - \frac{w - w_0}{w_c - w_0}\right) \quad (w_0 w \le w_c) ] [ t_n 0 \quad (w w_c) ]含义非常直观(w) 小的时候界面是线弹性刚度为 (K_n)达到最大牵引 (t_{max}) 后进入损伤软化到张开度 (w_c) 时界面完全丧失承载能力。断裂能是曲线下的面积双线性模型下等于 (t_{max} w_c / 2)。为什么选双线性而不是指数型或者梯形工程应用里双线性收敛相对容易参数仅有三个(t_{max})、(w_c)、(K_n)标定简单。对于实验数据不够充分的情况指数型牵引力的初始刚度在曲线形状过陡时会带来额外的非线性求解难度我试过一次之后就没有再用。2.3 裂缝张开度是耦合的粘合剂裂缝张开度 (w) 在这里是把两套理论粘起来的节点。XFEM 的富集自由度给出不连续位移场取法向跳跃得到 (w)粘聚本构把 (w) 映射成裂尖附近的牵引力流动方程又把 (w) 映射为裂缝内流道的水力开度因为裂缝内流量遵守立方定律[ q_f -\frac{w^3}{12\mu_f} \frac{\partial p_f}{\partial s} ]所以每次迭代我都要做三件事根据当前位移解计算 (w)把 (w) 传给裂缝内流动模块更新流导再根据 (w) 更新粘聚牵引力在右端项里的贡献。任何一步的 (w) 处理不一致矩阵都会出现不对称或刚度奇异。有一个很容易忽略的问题裂缝闭合时 (w) 可能趋近于零甚至变成负值。负开度在物理上意味着裂缝两侧接触流体流导应该置零同时要激活接触约束。如果不处理裂缝内压力稍一波动就会算出负的流导导致整个流动矩阵非正定。我在代码里对 (w) 做了下限截断和光滑过渡避免 (\frac{w^3}{12\mu}) 在零附近发生数值震荡。角色XFEM粘聚单元内聚模型描述对象裂缝几何与不连续位移/压力场裂尖前端断裂过程区的牵引-分离关系核心参数富集函数、节点富集标记(t_{max})、(w_c)、(K_n)对网格要求允许裂缝穿过单元单元尺寸需解析粘聚区长度输出裂缝张开、两侧位移/压力跳跃损伤变量、牵引力、断裂能消耗耦合贡献提供裂缝张开度与压力富集提供裂尖牵引与扩展判据3. Matlab 实现的第一步自由度、富集标记与矩阵装配3.1 先做归一化否则矩阵条件数让你崩溃有限元矩阵的病态问题在水力压裂里特别明显。如果直接使用国际单位孔隙压力单位 Pa渗透率单位 m²弹性模量单位 Pa黏度单位 Pa·s这些量相差十多个数量级。我最开始跑出来压力莫名其妙地振荡查了半天根本不是算法问题而是刚度矩阵条件数已经超过 1e14直接求解时精度完全毁了。建议在组装之前做一个无量纲化或者折中选择弹性模量用 GPa压力用 MPa位移用 mm时间用 min渗透率用 m² 且黏度用 mPa·s。这样至少把跨度压缩到两三个数量级内。如果还想更稳就引入特征长度 (L_0)、特征压力 (P_0) 和特征时间 (t_0)把方程全部除以特征量。对于刚开始搭建代码的人我的建议是先不管严格的无量纲推导直接用“GPa-MPa-mm-min”这一套折中单位跑通流程等需要对比不同工况时再回来做归一化。3.2 网格生成与富集标记对一个二维矩形压裂模型我通常用 8 节点四边形单元作为力学单元压力场用相同形函数但为了避免不可压缩锁死孔压场实际上采用连续的双线性插值。网格密度需要在裂缝路径附近加密初始裂缝长度大约 10 倍单元尺寸。每个单元给一个状态标记0普通单元1裂缝完全穿过所有节点都做 Heaviside 富集2仅裂尖处于单元内部启用裂尖富集3既有裂尖富集又有 Heaviside 富集。这个标记不是固定的。裂尖每推进一个单元就要把原来的裂尖富集单元转成 Heaviside 富集并把新的裂尖单元加上富集。在 Matlab 中这个逻辑用一个数组 status 和起点终点坐标来维护推进规则是如果裂尖到局部的裂纹扩展判据满足则在当前裂尖单元的相邻前方单元增加裂尖富集。一开始我把这部分逻辑写成了“每个时间步全局扫描一遍所有单元”效率很低。后来改成只扫描裂尖附近的几个候选单元运行时间直接降了一个数量级。3.3 自由度编号的预分配策略一个常见的坑是每步裂尖推进都重新编号自由度会导致矩阵结构反复变化稀疏模式重建成本极高。我采用预分配策略每个节点一开始预先分配所有可能的自由度普通状态下未激活的自由度对应行列为空在矩阵里只参与很小的附加项。激活与否用“自由度数”数组来标记。以平面问题为例每个节点存放标准位移自由度 2 个(u_x)、(u_y)孔压自由度 1 个可能的 Heaviside 位移富集自由度 2 个可能的 Heaviside 孔压富集自由度 1 个可能的裂尖富集自由度 8 个四个裂尖函数 × 2 个位移分量。看起来自由度总量非常大但实际激活的只占一小部分用 Matlab sparse 矩阵非常合适。下面是一段示意性代码它只是为了展示预分配的思想不是可以直接跑完的完整实现% 节点数 nNode每个节点预分配 fe_dof 个自由度 fe_dof 2; % 标准 ux, uy if addPressure fe_dof fe_dof 1; end % 富集实际允许的最大数量在组装时按激活情况决定是否列入稀疏组装 nodeDof repmat(fe_dof, nNode, 1); nTotalDof sum(nodeDof); globalDof zeros(nNode, fe_dof); idx 1; for i 1:nNode for j 1:fe_dof globalDof(i, j) idx; idx idx 1; end end % 用 activeDof 掩码决定该列是否填充未被激活的自由度不进入稀疏矩阵等式的有效行真实代码里还要区分每个节点是普通节点还是富集节点并用 activeDof 掩码决定该列是否填充。这套方案牺牲了一部分自由度总量换来了“每步裂尖推进时自由度编号不变化”的巨大便利扩展性非常好。3.4 富集单元的数值积分子区域剖分被裂缝穿过的单元形函数在单元内部不再光滑高斯积分点如果直接按矩形排布会漏掉几何不连续。我的做法是用水平集函数 (\phi(x)) 描述裂缝面。对于每个单元找到 (\phi0) 的等值线把该单元划分为一组三角形子区域。每个子区域取 3 或 7 个高斯点并进行坐标映射到母单元。这个几何剖分每次裂尖更新后要重新计算我把它预计算为积分点表和权重表存储在单元信息里避免重复运算。实测下来如果不用子区域积分哪怕 Heaviside 富集做对了裂纹扩展路径也会出现明显的锯齿状偏差因为裂纹一侧的积分贡献被错误计入另一侧。用子区域剖分之后裂缝扩展路径光滑了很多。这个步骤虽然繁琐但属于绕不过去的坎建议在代码里做成独立模块方便复用和测试。4. 流体-固体耦合推进从弱形式到 Newton 迭代4.1 控制方程的记忆锚点把方程组织成具有一定逻辑的方程组比死记公式有效岩石变形(\nabla \cdot \sigma \rho \mathbf{b} 0)用有效应力 (\sigma \sigma \alpha p_p I)岩石的伸长和缩短由位移决定基质渗流(\frac{\partial (\rho \phi)}{\partial t} \nabla \cdot (\rho \mathbf{q}) 0)其中 Darcy 速度 (\mathbf{q} -\frac{k}{\mu} \nabla p_p)裂缝内流动(\frac{\partial w}{\partial t} \frac{\partial}{\partial s} q_f q_l 0)(q_f) 用立方定律(q_l) 是滤失量粘聚/断裂(t_n f(w))控制裂缝端部破坏过程。在这个框架里未知量包括位移、基质孔隙压力、裂缝内压力以及裂缝张开度。全耦合意味着每个 Newton 步同时求解这些量而不是把它们拆分。从数学上看这些问题互相通过有效应力、流导与张开度联系在一起任何一个环节缺失都会导致结果失真。4.2 弱形式与矩阵分块结构取位移试验函数 (\delta u)、压力试验函数 (\delta p_p)、裂缝压力试验函数 (\delta p_f)将强形式乘以试验函数并在域内积分得到弱形式。离散后单个 Newton 增量满足下面的分块矩阵关系[ \begin{bmatrix} K_{uu} K_{ua} K_{up} K_{uf} \ K_{au} K_{aa} K_{ap} K_{af} \ K_{pu} K_{pa} K_{pp} K_{pf} \ K_{fu} K_{fa} K_{fp} K_{ff} \end{bmatrix} \begin{bmatrix} \delta u \ \delta a \ \delta p_p \ \delta p_f \end{bmatrix}\begin{bmatrix} r_u \ r_a \ r_p \ r_{pf} \end{bmatrix} ]这是一个非对称系统因为流体对位移的贡献中包含压力对张开度的敏感项。我的经验是在早期时间步它可以近似对称一旦裂缝扩展右侧的非对称部分会占主导。如果读者需要严格的黑箱调用应当选择不依赖对称性的求解器或先对矩阵做缩放。具体到某些子矩阵的含义(K_{uu})力学刚度矩阵包含有效应力项(K_{up})孔压与位移的耦合来自有效应力原理(K_{uf})裂缝压力对裂缝两侧的作用力(K_{pp})基质渗流的刚度矩阵(K_{pf}) 与 (K_{fp})基质-裂缝流量交换项和裂缝-基质压力交换项。组装时使用统一的单元循环先判断单元类型再叠加对应的子矩阵。对仅有 Heaviside 富集但没有裂尖富集的单元需要特别小心富集自由度的列应只在本单元内起作用不能让其他单元的节点富集自由度错误地产生贡献。4.3 全耦合 Newton 迭代与自适应时间步求解流程为给定当前状态 (U_n)、时间步 (\Delta t)初始化当前迭代解为 (U_n)计算残差 (R(U^k)) 和雅可比矩阵解 (J \Delta U -R)更新当前解如果残差范数小于容差接受本步解进入下一时间步否则继续迭代如果超过最大迭代次数把 (\Delta t) 减半重试。这里有几个实际调试点。残差范数要按“压力残差”和“位移残差”分开监控。因为两者量纲不同放在一起容易被大数掩盖细节。我经常打印 (R_u / |R_u|) 和 (R_p / |R_p|) 两个量如果位移残差已经很小而压力残差仍然很大说明问题出在渗流刚度矩阵。求雅可比可以用解析推导也可以用有限差分近似。我强烈建议先用有限差分验证解析雅可比用小的三节点网格模型做这个验证一次能发现很多 bug。否则好不容易拼出一个大矩阵结果里面某一行符号写错排查起来非常痛苦。时间步长控制上我要求每个时间步内流体前沿移动距离不超过一个单元尺寸。如果超过裂缝内压力对张开度的响应会出现锯齿。自适应时间步用最朴素逻辑就行Newton 迭代步数小于 3 就增大 1.2 倍大于 8 就减半超过 12 直接判定发散。5. 粘聚本构与参数标定最容易拖垮收敛的环节5.1 双线性牵引-分离关系的三个关键参数粘聚本构的物理含义很简单但参数标定是压裂模拟里最容易翻车的部分。需要确定三个量。最大牵引 (t_{max}) 和岩石抗拉强度直接相关。对于泥页岩单轴抗拉强度通常在 1-10 MPa 量级。这个值不能拍脑袋设置最好来自巴西劈裂或直接拉伸实验。断裂能 (G_c) 可以通过三点弯曲或半圆弯曲实验获得量级通常在几十到几百 J/m²。初始刚度 (K_n) 则没有天然实验量它只是一个数值参数。选择要满足两个互相矛盾的条件足够大以避免背景单元插入虚假柔度又不能大到导致严重病态。试了很多之后我一般取 (K_n (10\sim 50) \times E/h)这里的 (h) 是单元尺寸。一个快速标定起点可以参考下表参数典型范围依据(t_{max})1-10 MPa抗拉强度实验(G_c)20-500 J/m²断裂韧性/断裂能实验(K_n)((10\sim50) E/h)数值刚度需调初始裂缝半长网格数的 10-20 倍避免裂尖初始边界影响如果你拿到的实验数据只有断裂韧性 (K_{IC})可以用 Irwin 公式粗略换算(G_c \approx K_{IC}^2 / E)其中平面应变下 (E E/(1-\nu^2))。换算出来的断裂能只能作为初始估计真正用于模拟前最好再对照经典压裂算例做一次校准否则积分路径上的能量耗散可能差出好几倍。5.2 粘聚区长度与网格分辨率的匹配粘聚区的长度约等于[ l_{cz} \approx \frac{E G_c}{t_{max}^2} ]这个概念被很多人忽视。如果你取 (E30) GPa、(G_c50) J/m²、(t_{max}2) MPa那么 (l_{cz}) 大约是 0.375 m。也就是说单元尺寸必须远小于 0.375 m才能把这个过程区解析出来。工程压裂模型的网格往往在几十米尺度所以直接用固定网格完全解析不了粘聚区。我的解决方案是在裂尖区域做局部加密或者使用自适应网格。但对 Matlab 小模型来说直接假设一个较小的区域、加密到厘米级最省事。网格不足时可以看到的现象是损伤沿网格扩散、压力曲线抖动、断裂能消耗被高估。这些都是明显信号。5.3 流体前沿和裂尖位置必须匹配压裂过程中有一个容易被忽略的细节流体可能滞后于裂缝尖端。裂缝头端有一部分是干的流体前沿不一定到达裂尖。如果数值模型强行让裂缝内充满了水压就不能体现这种滞后。对注入排量较大或高滤失的情况流体前沿追赶裂尖很快可以假设两者重合但在低排量长渗透率工况滞后区域会显著影响压力响应。实现层面追踪流体前沿用的是裂缝内另一个方程当裂缝内压力 (p_f) 和宽度 (w) 都给定时流体的质量守恒会决定前沿位置。最简单的做法是在每个单元时间步里监测裂缝张开度和压力若某单元裂缝宽度很小但压力已为 0说明流体前沿还没到达该单元。此时该单元的流导应置零而不是用很薄的膜来模拟。这个处理对收敛影响极大。我第一次模拟时没管它裂缝尖端附近出现一个无法解释的压力高尖峰后处理时还以为是 XFEM 富集的问题最后发现只是流体前沿推进逻辑写得不对。加上滞后判断之后压力曲线立刻平滑下来。6. 实际调试记录压力振荡、零能模式与时间步控制6.1 压力振荡裂尖附近的孤点尖峰我复现过程中第一次看到压裂压力曲线在裂尖单元上跳动时第一反应是检查压力富集的实现。后来定位后发现孔压富集自由度在裂尖单元中存在近似零能模式尤其当压力试验函数与位移富集函数在同一单元重叠时会出现一些小本征值反映到结果上就是压力的非物理振荡。缓解手段只对裂缝贯穿单元启用压力 Heaviside 富集裂尖单元内不启用压力富集而采用界面局部渗流模型给裂缝压力自由度引入数值稳定项类似流线迎风对裂尖单元的位移富集在粘聚刚度中叠加一个很小的正定补充阵压制零能模式。在调试时我习惯把单调下降的排量、初始缝长、恒定注入作为基线算例。一旦压力曲线出现无物理意义的振荡立即用单个激活单元做局部检查靠二分定位到具体单元。这种方法比盯着整条曲线猜要快得多。6.2 零能模式负刚度和鬼影位移XFEM 的富集自由度有时会产生零能模式表现为单元位移场中出现“鬼影”模式即使没有外力也能产生非零位移。在纯弹性 XFEM 中通过 shift 富集函数可以缓解但在裂尖富集单元中如果裂尖函数族含某些高阶项仍然可能产生局部模式。我的一个局部策略是在裂尖单元的刚度矩阵对角上加上一个与 (E/h) 同量级的小数即所谓的“稳定项”。但是这种人为增广必须控制大小过大会改变真实刚度响应。更稳妥的办法是从源头选好富集函数当单元内裂尖接近单元边时把裂尖富集退化为 Heaviside 富集等价于裂尖进入下一个单元后再启用裂尖富集。我经过几次对比这样做直接消除了大部分鬼影位移而且对裂缝扩展路径没有产生可观察的偏差。6.3 调试流程与提速建议到目前为止我把整个调试流程固定成了四步强烈建议照着走先验证一个纯力学问题一个含初始直裂缝的板施加固定拉伸看 XFEM 的裂纹张开与已知解析解是否一致再验证一个纯流体问题给定裂缝宽度分布施加恒定入口压力看裂缝内压力分布与立方定律解是否匹配固定裂缝长度、只开注水过程看压力、张开度、滤失是否随排量单调变化最后才允许裂尖扩展。运行速度方面Matlab 循环本身不是瓶颈的主要来源真正的瓶颈是每次迭代对大矩阵的重构。所以我把所有不变化的部分形函数导数、高斯权重、单元连接全部预先集成并放在 cell 里每次迭代只更新随状态变化的项。这样一轮 Newton 迭代的矩阵组装时间可以减少一半以上。另一个提速点残差下降缓慢时先检查缩放再检查雅可比而不是盲目减少时间步。当某个区域的渗透率极低时渗流矩阵条件数会非常大此时即便雅可比正确直接求解器也可能因为舍入误差而失效。我的经验是给渗流矩阵加一个小的穿孔项或人工扩散或者把时间步适当调大而不是调小反而能改善收敛。6.4 小规模算例与观测指标我用一个 1 m × 1 m 的二维模型中心 0.2 m 初始裂缝顶部与底部固定左侧与右侧加载 3 MPa 均匀应力注入流量 (1\times10^{-4}) m²/s。网格用约 60×60 的四边形单元裂尖附近局部加密到单元尺寸 2 mm。运行大概 20 个时间步每步 3 次 Newton 迭代。需要记录的主要输出有井口压力随时间曲线裂缝张开度沿裂缝路径分布裂尖位置随时间变化粘聚单元损伤因子云图基质孔压场和裂缝压力场。有些模拟器把输出集中在最终一张图上但如果压力曲线本身不能复现解析解云图再漂亮也没有意义。所以我每次在模型里放一个观测点输出时间序列用它与经典解析解的对比曲线作为诊断工具。一旦压力曲线偏离解析解第一件事不是去调网格而是去检查本步的 Newton 迭代是否真的收敛了以及时间步是不是让流体前沿跨过了太多单元。跟 XFEM 和粘聚单元纠缠了几个月我最深的体会有两点。第一算法选型固然重要但真正决定一个压裂模拟器能不能用起来的往往是最琐碎的细节自由度怎么编号、状态变量从哪个单元继承、时间步怎么控制。第二不要一上来就追求全耦合和复杂几何。先把最简单的 2D 单裂缝模型跑通再逐步增加压力富集、流体前沿追踪和复杂载荷调试成本会低很多。如果只让我分享一个技巧我会选这个在模型中间预先放一个观测点同时输出压力、张开度和流体前沿位置的时间序列。你的代码一旦能稳定再现经典的平面应变解析解曲线后续往三维、非均匀应力场、多裂缝方向扩展基本就是按部就班的事。接下来我打算把同一套框架扩展进 Matlab 的并行求解路径初期结果看起来很有希望等稳定后再单独写一篇记录。