干水力压裂数值模拟这块儿，绕不开的话题就是相场法。这两年用 Comsol 做相场断裂的案例越来越多，但大多停留在各向同性介质，一旦碰上页岩、层状岩体这类横观各向同性介质，很多默认的设置就直接失灵了。这个项目就是一次完整的实战记录：从材料本构的坐标设定、相场-流动耦合的控制方程，到求解器调试、裂纹偏转路径分析，全程用 Comsol 走了一遍。如果你正准备用相场法算水力压裂，又对横观各向同性材料怎么建模心里没底，这篇可以直接当操作手册用。
1. 为什么横观各向同性介质要把压裂单独拎出来讲
1.1 横观各向同性到底是种什么材料
先把这个概念掰扯清楚。横观各向同性（Transverse Isotropy）听起来吓人，其实一句话就能说明白：材料在某个平面内所有方向性质都相同，但沿着这个平面的法线方向，性质不一样。
最典型的例子就是页岩。层理面方向（平行于沉积面）和垂直于层理面的方向，弹性模量、渗透率、断裂韧性差别都很大。实际工程测量中，页岩水平方向的杨氏模量可能达到 30 GPa，垂直方向只有 15 GPa，这直接决定了水力裂缝怎么起裂、怎么偏转。
在 Comsol 里描述这种材料，需要用到正交各向异性弹性模型，但横观各向同性比一般的正交各向异性更特殊：5 个独立弹性常数就够了。按坐标定义来说，如果取 x-y 平面为各向同性面，z 轴为材料对称轴，那么独立参数就是：
- E₁ = E₂：各向同性面内的杨氏模量
- E₃：垂直于层理面方向的杨氏模量
- ν₁₂：面内泊松比
- ν₁₃：面内主应力方向与 z 方向之间的泊松比（或者说 ν₃₁）
- G₁₃：沿 z 方向的剪切模量
注意一个容易踩坑的点：面内剪切模量 G₁₂ 不是独立参数，它由 E₁ 和 ν₁₂ 推导出来，G₁₂ = E₁ / [2(1+ν₁₂)]。如果你在 Comsol 里手动填了 G₁₂，又填了 ν₁₂，数值不自洽的话，计算得到的弹性矩阵可能不是物理上允许的。
1.2 水力压裂在这个场景下的特殊之处
各向同性介质里，水力压裂的裂纹路径比较直白，最大主应力方向控制起裂和扩展。但横观各向同性介质有两个关键差异：
第一是应力耦合效应。层理面的存在使得裂纹尖端的应力场不再对称，裂纹扩展方向会同时受三方面控制：地应力差、材料弹性各向异性、断裂韧性各向异性。三者之间的竞争关系，决定了裂缝是沿直线走、偏向层理面、还是被弱面“俘获”后沿着层理面转向。
第二是渗流-应力-损伤的强耦合。水力压裂是靠流体压力撑开裂纹的，而裂纹一旦张开，局部渗透率会从基质渗透率（比如 10⁻¹⁵ m² 量级）跳到裂缝渗透率（可能到 10⁻¹⁰ m² 量级），跨越五个数量级。这种强非线性在数值实现里非常容易引起求解器崩溃，这也是为什么很多人用扩展有限元（XFEM）做水力压裂时，裂缝内压力与位移场的耦合要写一堆自定义方程，而相场法对这一类问题的处理反而更优雅。
1.3 为什么选相场法而不是扩展有限元
我不止一次被问：既然 XFEM 也可以模拟裂纹扩展，为什么还要用相场法？我的回答是：看你的问题到底关注什么。
XFEM 的优势在于裂纹几何描述清晰，裂缝面是显式的。但它处理分叉、交叉、裂纹转向时非常麻烦，因为每个新裂纹都需要重新配置富集函数。而对于水力压裂，裂纹很可能沿着层理面发生偏转甚至分叉，这个时候 XFEM 的前处理工作量会爆炸。
相场法走的是另一条路：把裂纹面用一个连续的相场变量 d 来表示，d = 0 代表完好材料，d = 1 代表完全断裂。裂纹的起裂、扩展、分叉统统由相场演化方程自动计算，不需要显式跟踪裂纹面。代价是网格尺寸要足够细，计算量偏大，但换来的是几何处理上的极大简化。
Comsol 内置的“相场-断裂”接口（Phase Field, Fracture）就是基于这个思想实现的。配合固体力学、达西定律或者裂隙流模块，可以比较自然地搭出水力压裂的耦合框架，这也是这个项目选 Comsol 的核心原因。
2. 相场法与水压耦合的核心思路
2.1 控制方程到底长什么样
我需要先把物理框架讲透，因为很多人在 Comsol 里设置不对，根本原因是不知道自己在解什么方程。
相场断裂的核心假设，是把总能量分成两部分：弹性应变能 + 断裂能。考虑材料损伤之后，弹性应变能会被一个退化函数“打折”，这个退化函数最常用的是 g(d) = (1 - d)² + κ，κ 是一个很小的数值稳定性参数，防止完全退化后刚度矩阵奇异。
于是应力更新变成：
σ = g(d) · C⁰ : ε
其中 C⁰ 是完好材料的弹性张量。对于横观各向同性介质，C⁰ 就是由那 5 个独立弹性常数构成的正交各向异性弹性矩阵。
相场演化方程是一个梯度型方程，Comsol 内置的形式大致是：
Gc/l₀ · (d - l₀² ∇²d) = 2(1-d)H
这里的 Gc 是断裂能，单位 J/m²；l₀ 是正则化长度参数，控制裂纹扩散带的宽度；H 是历史驱动变量，由应变能中受拉部分的峰值决定。为什么叫历史变量？因为裂纹不能愈合，所以必须记录历史上达到过的最大驱动强度，而不是用当前瞬时值。这一点在 Comsol 默认实现里已经处理好了，但如果你自己写弱形式就千万别漏。
水力压裂的“水力”部分，我用的是达西定律描述孔隙流体流动，但渗透率是随相场变量变化的：
k(d) = k₀(1-d) + k_f·d
k₀ 是基质渗透率，k_f 是裂纹完全打开后的裂缝渗透率。流体压力 p 驱动裂纹张开，裂纹张开反过来改变渗透率，这是一个双向耦合。
2.2 横观各向同性如何进入相场框架
很多人在这里犯迷糊：材料各向异性到底应该放在哪个环节？
答案是：放进弹性张量 C⁰ 里。断裂能不能简单做成标量，也要区分方向。我在这个项目里把断裂能设成了两个方向不同的值：Gc₁₂ 是面内断裂能，Gc₃ 是沿层理面法向或者说层理面的断裂能。层理面往往是弱面，所以 Gc₃ 可以取 Gc₁₂ 的一半甚至三分之一。这样裂纹遇到弱面时，从能量角度来看，转向沿着弱面扩展更加“便宜”，计算出来的结果自然会出现偏转、俘获这些实际工程现象，这正是横观各向同性介质压裂模拟最想看到的东西。
2.3 为什么不能直接套用各向同性案例的设置
在 Comsol 的案例库和网上教程里，相场断裂大多是用各向同性材料演示的。直接拿过来改个各向异性材料参数，通常会出现两类问题。
第一类问题是泊松比和剪切模量的工程师输入方式。Comsol 里各向异性材料可以填工程常数，也可以填写弹性矩阵的所有分量，但软件不会替你做横观各向同性的一致性检查。如果你填了 E₁=30GPa、E₃=15GPa、ν₁₃=0.25，却没有按对称性条件修正 ν₃₁，那么矩阵实际上不是横观各向同性的，甚至可能出现负的模型刚度。
第二类问题是坐标系的取向。各向同性材料无所谓坐标，但各向异性材料极端依赖坐标。层理面平行于 x-y 平面，就要确保材料坐标系的 z 轴垂直于模型平面。如果几何体不小心旋转了或者装配坐标系出问题，裂纹扩展方向会完全乱掉。
3. Comsol 模型搭建全流程实操
3.1 几何建模与初始裂缝定义
我的模型是一个二维平面应变问题，尺寸 1 m × 1 m，中心有一条水平初始裂缝，半长 0.05 m，给一个 0.1 m 长的切缝作为初始裂纹。
关于初始裂纹的处理，有两种常见做法，效果差不多：
- 直接在几何里用分割线把初始裂缝切出来，在裂缝面上施加流体压力边界
- 不切缝，而是在初始裂缝位置预置一个初始相场损伤区，比如令 d = 1
我实际用的是第二种办法，因为更符合相场法“扩散裂纹带”的思维方式，也不用处理裂缝面上网格不连续的问题。具体操作是给几何添加一个“初始值”节点，把相场变量的初始值设置为：
d = 1，如果你在初始裂缝矩形区域内 d = 0，其余区域
这个矩形区域的宽度取一个网格单元宽度就够了。
然后施加的边界条件包括：
- 模型四周：法向位移固定，应力边界远处施加远场地应力 σx、σy
- 中心注液位置：给定一个注入速率 Q
- 相场边界：默认零通量
远场地应力我设置为各向异性的：σy（最大主应力方向）取 -6 MPa，σx 取 -3 MPa，都是压缩应力。这里有个需要注意的约定：压缩应力在岩石力学中通常带正号，但 Comsol 固体力学默认以拉伸为正，所以计算时要把实际的地应力换算成对应的分量。
3.2 材料参数与坐标设置：横观各向同性怎么填
材料参数是这个项目的核心，我列一张表，方便你直接抄作业：
| 参数 | 数值 | 说明 |
|---|---|---|
| E₁ = E₂ | 30 GPa | 层理面内杨氏模量 |
| E₃ | 15 GPa | 垂直层理面模量，通常更软 |
| ν₁₂ | 0.25 | 面内泊松比 |
| ν₁₃ | 0.25 | 面内与法向之间的泊松比 |
| G₁₃ | 8 GPa | 法向剪切模量 |
| G₁₂ | 12 GPa | 由 E₁、ν₁₂ 计算得到 |
| Gc₁₂ | 150 J/m² | 面内断裂能 |
| Gc₃ | 80 J/m² | 沿层理面断裂能，弱面 |
| l₀ | 0.005 m | 相场正则化长度 |
| k₀ | 1e-15 m² | 基质渗透率 |
| k_f | 1e-10 m² | 裂缝渗透率 |
| Q | 2e-4 m²/s | 二维注入速率 |
| ρ | 1000 kg/m³ | 流体密度 |
| μ | 5e-4 Pa·s | 流体动力黏度 |
这里要特别注意坐标约定：材料坐标系下，x-y 平面是各向同性面，z 轴沿层理面法向。如果模型是二维的 x-y 平面，那么模拟的是“垂直于层理面的截面”，也就是裂缝在这个截面里可能沿着层理面方向偏转；如果模型是 x-z 平面，同样合理，只是各向同性面变成了边界方向。我在模型里是把二维平面理解为 x-z 平面，层理面沿水平方向，这样裂缝偏转看起来就像现实中看到的水平缝转垂直缝（或者反过来）。
在 Comsol 的“固体力学 > 线弹性材料”节点里，弹性模型选“正交各向异性”，然后按表格填进去。特别提醒：ν₃₁ 不要按自己的想法乱填，让 Comsol 的对称化设置自动处理，或者在填写时单独设一个变量 ν₃₁ = ν₁₃·E₃/E₁，保证弹性矩阵对称。
3.3 物理场接口与耦合关系怎么搭
这个模型一共用到三个物理场接口：
- 固体力学（Solid Mechanics）
- 相场断裂（Phase Field, Fracture）
- 达西定律（Darcy's Law）
在 Comsol 的组合里，通常是把相场接口叠加在固体力学之上，它们之间通过“边界/域”的耦合变量自动连接。达西接口则负责算流体压力场。
关键的一步在于把手动耦合写准确。我在“变量”节点里定义了以下两个核心表达式：
压力荷载：固体域内把流体的压力场 p 转化为体积力或边界力，作用区域用损伤变量 d 加权。具体表达式是：F_p = -(∇p) 对应的稳态项要写成总应力平衡里的附加项。对于达西流和水力裂缝，更标准的是把 p 加到固体应力平衡方程里，通过“孔隙压力”子节点实现。
渗透率更新：在达西定律的渗透率节点里，把渗透率定义成 p 的表达式：k_iso = k0 + (kf-k0)*d。这样损伤越严重，渗透率越大，流体就越容易流入裂缝区域。
不过要注意，如果直接在达西定律域内用这个表达式，那么压力分布会表现为：裂缝尖端高流体压力驱动裂纹扩展，但裂缝尖端之外的区域，渗透率很低，压降明显。这既模拟了裂缝内流体压力传递，也模拟了基质的滤失效应，比单纯在初始缝面加压力边界要真实得多。
3.4 网格设置：相场法成败的分水岭
网格是相场法最需要花时间的地方。l₀ = 5 mm，那么裂纹扩散带宽度大约是 2 到 3 倍 l₀，网格尺寸必须小于 l₀ 的一半，也就是至少要 2.5 mm 以下，才能在裂纹带上分辨出相场的梯度变化。
对于 1 m × 1 m 的模型，如果全区域都细化到 2.5 mm，网格数和自由度会相当可观。我的方案是：先在初始裂缝可能扩展的带状路径上，用“映射”方式切分出一个宽度约 0.2 m 的细化区，细化区里网格尺寸控制在 1 mm 到 2 mm；其他区域自由三角形网格，尺寸 5 cm。
网格质量检查别忘了关注最小单元质量，相场求解对网格畸变很敏感，如果某个单元质量低于 0.3，收敛性会明显变差。
3.5 求解器设置与时间步控制
相场断裂是高度非线性问题，而且和达西流动耦合，直接一锅端全耦合求解非常容易不收敛。我实际跑下来的经验是分两步分离式求解更稳：
第一步：固定损伤场，求解流体压力和固体位移。 第二步：固定流体压力，求解相场变量 d 和固体位移。
在 Comsol 的“求解器配置”里，把物理场按“Darcy”和“固体+相场”两个组分别勾选，开启“分离步”，指定每个步的容差。
时间步方面，不要用太大步长。我设置的是自适应时间步长，初始步长 0.001 s，最大步长 0.05 s，总模拟时间 10 s。压裂问题的前期注入阶段，压力积累比较慢，可以把步长放小一点；一旦裂纹开始起裂扩展，压力波动加剧，步长要更小才能捕捉。
还有一个重要技巧：一开始不要让注液速率直接达到最终值，用一条斜坡加载曲线，在 0.1 s 内从 0 线性增加到 Q。这样相当于给系统一个“软启动”，大幅减少初始阶段的冲击载荷导致的不收敛问题。
4. 结果解读与典型现象验证
4.1 裂纹偏转路径：各向异性最直观的体现
模拟跑完，先不看云图，直接看损伤场 d = 0.5 的等值线，这就是裂纹面的“等效位置”。
在横观各向同性设置里，如果我只把弹性模量设成各向异性（E₁/E₃ = 2)而断裂能保持各向同性，裂纹基本还是会沿着最大主应力方向直走，只是偏转角度略有变化。这说明弹性各向异性对路径的影响是渐变式的。
但一旦把断裂能也设为各向异性，Gc₃ 远小于 Gc₁₂，裂纹前进到一定距离后，会在层理面附近发生明显转向，最终被弱面“俘获”，沿层理面方向扩展。这个现象在地质力学里就叫“裂缝沿弱面转向”，压裂施工中如果地应力差不够大、弱面强度又低，压裂液就容易顺着层理面跑，形成复杂的非平面裂缝。
从这个模拟结果你能清楚看到：相场法不需要任何额外的转向判据，能量最低原理自动决定路径，这是相场法最大的价值所在。
4.2 注入压力曲线：工程判断的关键依据
典型的相场水力压裂模拟结果，注入点压力随时间曲线长得像这样：
- 初期线性上升：流体注入，裂缝未起裂，孔隙压力积累
- 峰值点：达到起裂条件，裂纹开始扩展，压力突然下降
- 后期波动：裂纹扩展过程中，由于材料各向异性和弱面转向，压力出现锯齿状波动
如果得到的曲线没有这个“峰值后回落”的特征，大概率是模型哪里出了问题——要么是相场损伤初始值设置不当导致裂纹一开始就扩展，要么是网格不够细导致起裂提前。
这个压力曲线的工程意义非常大：压裂施工设计中，地面泵压的预测全靠这种模拟结果。相场法给出的压力曲线比某些简化模型更接近现场压裂施工的真实泵压波动规律。
4.3 参数敏感性：哪些因素对结果影响最大
我做了一轮参数扫描，重点看了三个变量：
- E₃/E₁ 比值从 1 降到 0.3
- Gc₃/Gc₁₂ 比值从 1 降到 0.4
- 远场应力差 σx-σy
结论比较有意思：E₃/E₁ 的影响主要体现在起裂位置和裂纹宽度上，模量差距越大，裂缝宽度越不均匀；但真正决定裂纹转向“干脆不干脆”的，主要是断裂能比值。当 Gc₃/Gc₁₂ 降到 0.5 以下时，裂纹被弱面俘获的概率显著上升，即使应力差有利也不一定能维持裂缝直线扩展。
这个发现其实和现场施工的经验是吻合的：页岩地层里压裂施工经常出现微地震事件分布非常离散，很大程度就是因为层理面的断裂能差异在起作用。
5. 常见问题与排查技巧实录
5.1 求解不收敛：先找三个位置
相场-水力压裂模型不收敛，绝大多数情况下不是求解器设置的问题，而是模型定义的问题。按照出现概率排序，我建议依次检查：
渗透率表达式是否引入不连续。如果 k(d) 的表达式在 d 变化时引起压力场突变，达西求解就会振荡。解决办法是在渗透率表达式里加一个很小的过渡区间，比如用 (tanh((d-0.7)/0.1)+1)/2 代替硬截断。
弹性矩阵是否正定。填完横观各向同性参数后，在结果里加一个“弹性矩阵特征值”的探针，看看有没有负特征值。如果出现负特征值，说明参数组合不合理。
初始损伤区是否在应力场下产生数值振荡。初始 d=1 的区域内，刚度已经退化到 κ 的量级，如果 κ 取得太小（比如 1e-6 以下），就会出现大位移振荡。把 κ 设置在 1e-4 到 1e-5 之间通常比较稳。
5.2 裂纹路径看起来“糊”或者“碎”：网格还是太粗
损伤云图如果看起来是一条宽的模糊带，而不是清晰的裂纹面，第一反应不应该是调求解器，而是剖网格。l₀ 取 0.005 m 时，必须保证裂纹扩展路径上网格边长不超过 l₀/2。达不到这个条件，相场扩散带会在数值上被人为加宽，裂纹路径就会失真。
网格太粗的另一个表现是裂纹扩展方向出现锯齿状偏移，因为相场梯度受网格诱导，会沿着网格线走“之”字形。修正办法就是细化路径网格，但这个成本确实高。如果计算资源有限，建议把模型尺寸缩小到 0.5 m × 0.5 m，l₀ 同比例缩小，这样网格数量相对可接受。
5.3 注入点压力异常：检查你的“注入”到底注到哪儿了
有些读者反馈：压力一直不涨，或者涨得特别快。这两种情况背后通常是同一个原因——注入点和模型的连通关系没设对。
如果注液点设成了 Dirichlet 压力边界（固定压力），那压力当然不会涨，因为它是被固定住的。正确做法是设成 Neumann 型的注入流量边界。如果用的是达西接口，在“流动”边界里选“流入/流出”并给定质量流量或体积流量。
反过来，如果压力涨得飞快，说明流体被“堵住”无法排开裂缝，多半是初始裂纹区域还没有完全损伤，渗透率太低。给初始损伤区的渗透率一个高地步值，比如直接用 k = k_f，让流体能先在这个区域流动起来，压力曲线就会正常。
5.4 计算时间太长：几个实用降本技巧
相场法的计算成本是出了名的门槛，尤其是三维模型。如果只是做二维参数扫面，可以从这几个方面压计算量：
- 利用对称性，只建 1/4 模型
- 相场和时间步的容差适当放宽，先用粗网格试出裂纹的大致走向，再在裂纹路径附近局部细化
- 固定相场求解的次数，不要每个时间步都完整求解相场方程，有些步里可以让相场保持前一步的值
- 用“分离式求解器”中的“迭代次数限制”选项，把每个时间步的相场迭代上限设到 5 次左右
实测下来，这几招能把整体计算时间压缩到原来的三分之一，而裂纹路径几乎不变。
6. 给新手的一些心里话
我个人在这个项目上踩过的最大一个坑，就是一开始迷信“全耦合才是精确的”。实际跑下来发现，对于水力压裂这种渗流-应力-损伤三重非线性耦合的问题，分离式求解不仅稳，而且思路更清晰——先让流体压力找到平衡，再让固体和损伤对压力做出响应，否则你根本分不清数值振荡到底来自哪个物理过程。
另外，相场法的核心不是 Comsol 操作，而是你对网格尺寸、正则化长度、断裂能这三个参数之间关系的理解。l₀ 取大一点，裂纹带就宽，断裂能就被“稀释”了；取小了，网格就要密到爆炸。找到平衡点的唯一办法就是多做几组敏感性测试，把 l₀ 和网格尺寸的关系摸透。
如果你也是刚把横观各向同性材料引入相场水力压裂模型，建议先从各向同性材料开始，把注液-起裂-扩展这套流程跑通，再逐步加入各向异性弹性、各向异性断裂能。每一步加一个变量，出了问题也好定位。磨刀不误砍柴工，这个顺序值得遵守。
