企业官网建设流程全解析

不少做数值计算的朋友第一次看到“相场法模拟水力压裂”这个标题，第一反应都是：COMSOL 里真的能算裂缝扩展吗？能和流体压力耦合吗？算出来到底准不准？我最近正好把这一系列一共整理了 6 个案例，从最基础的单一裂缝扩展做到多裂缝、非均质、孔隙弹性耦合甚至三维形态，中间踩了不少坑，也把参考文献从头到尾筛了一遍。这篇就把相场法在水力压裂里的实际建模思路、COMSOL 实现路径、案例设计逻辑和收敛经验一次讲清楚，适合正在做页岩储层压裂、地热储层改造、混凝土水力损伤这类方向的同学参考。

1. 为什么用相场法算水力压裂：从“裂缝在哪”到“裂缝怎么长”

1.1 断裂力学里的两套语言：离散裂纹 vs 弥散裂纹

传统有限元做裂缝扩展，最常用的思路是离散裂纹模型。你在几何里预设一条裂缝面，裂缝尖端两侧的网格是独立的，裂缝扩展一步，就要重新连接单元或者插入新的界面单元。这个方法物理概念清楚，但在水力压裂场景里非常吃力：裂缝方向不是预先知道的，它随应力场、流体压力和岩石非均质性实时变化，还要处理裂缝分叉、交叉、多裂缝同时扩展，每步都做网格重构，前处理成本高到让人怀疑人生。

相场法走的是另一条路。它把裂缝看成一种连续变化的损伤场，用一个标量场变量 d 表示材料在某个位置“碎”的程度，d=0 是完全完好，d=1 是完全断裂。裂缝宽度不再是几何上的开口距离，而是损伤带的宽度，这个带可以跨几个单元，厚度由人为引入的特征长度参数控制。打个比方，离散裂纹像是在地图上画一条细线，相场法更像是用热力图去显示人流密度，虽然没有一条“线”，但密集区的位置和形状清清楚楚。

1.2 水力压裂问题为什么特别适合相场法

水力压裂的核心矛盾在于：流体注入改变了孔隙压力，孔隙压力改变了有效应力，有效应力又控制着裂缝是否扩展；裂缝一旦扩展，又反过来给流体提供新的流动通道。这是一个强多物理场耦合问题，而且裂缝形态完全是自由的，还经常出现多裂缝井组之间的应力阴影效应。

相场法恰好把“裂缝拓扑变化”这个最麻烦的部分消掉了。你不需要追踪裂缝边界，不需要判断尖端在哪，不需要处理分叉时的几何拓扑变动，裂缝路径完全由能量最小化原则自动决定。这在 COMSOL 里面做多物理场耦合尤其顺，固体力学模块算位移和应力，Darcy 或孔隙弹性模块算压力，相场变量描述损伤，三套变量在同一个网格上互相影响，省掉了所有界面传递操作。

当然，代价也是明显的。相场模型的计算域必须覆盖裂缝可能到达的整个区域，所有单元都要额外求解一个相场偏微分方程，计算量比局部重网格方法大。而且相场长度尺度 l 和有限元网格尺寸 h 之间存在强约束，网格太粗裂缝会走样，网格太细又算不动，这就是我后面要专门展开讲的“网格依赖”问题。

2. 相场模型的核心方程与COMSOL实现思路

2.1 Griffith能量与相场变量 d 的含义

相场断裂模型的物理根基是 Griffith 能量准则：裂缝扩展的驱动力来自弹性应变能释放，阻力来自材料的新生表面能。把裂缝表达为弥散带后，总能量可以写成体积积分形式，其中包含弹性应变能、断裂表面能以及损伤梯度项。标量场 d 满足一个偏微分方程，最常见的形式是：

d - l^2 \nabla^2 d = 0 （当弹性驱动力超过阈值时，源项激活）

在 COMSOL 里可以直接用“系数型 PDE”接口来写这个方程，也可以使用固体力学模块内置的“脆性断裂（Brittle Fracture）”多物理场特征。内置接口的好处是它已经帮你做了历史场变量的存储和应变能分解，自己写 PDE 也不是不行，但要小心符号和单位，尤其是 l 的单位是长度，方程里各项必须严格匹配。

一个关键点是应变能分解。压缩区域的应变能不应该驱动裂缝扩展，否则你会看到裂缝横着穿过去，这完全不物理。常见做法是把应变能拆成拉伸和压缩两部分，只有拉伸部分参与损伤演化。COMSOL 内置的脆性断裂接口默认支持这类分解，但如果你自己写控制方程，很容易漏掉这个细节，导致结果怪异。

2.2 水力驱动的“力”是怎么加进去的

水力压裂和普通脆性断裂最大的区别在于载荷来源。普通断裂实验是给边界加位移或力，水力压裂是靠注入流体在裂缝内部形成高压，高压把裂缝撬开。这个压力一方面受岩石渗透性影响，另一方面又受裂缝宽度影响——裂缝张开，流体流入更多，压力变化又继续改变裂缝形态。

最常用的耦合方案是把多孔介质中的流体压力作为独立场变量，通过 Biot 有效应力原理进入固体动量方程。有效应力等于总应力减去 Biot 系数乘以孔隙压力，孔隙压力本身由质量守恒方程控制，流动参数取决于孔隙度和渗透率。当材料局部损伤接近 1 时，可以适当增大局部渗透率，模拟裂缝作为高导流通道的效果。这套做法在 COMSOL 里可以组合“固体力学”＋“Darcy 定律”或“空隙介质中的流动”接口，再通过“多物理场耦合”节点连接。

如果你是刚入门，还可以做一个更简化的版本：不让流体在孔隙中扩散，而是把注水压力直接作为随时间变化的边界载荷，作用在初始裂缝上。这个模型算起来快，适合先把相场本身的扩展行为跑通，再逐步加入真正的流固耦合。我在第一和第二个案例里就采用了这种简化，到第四个案例才加入完整的孔隙压力耦合。

2.3 COMSOL里需要启用的模块与物理场接口

我建议最低配置是“结构力学模块”加“PDE 模块”，如果要算孔隙弹性耦合，还需要“地下水流模块”或者“多孔介质模块”。具体物理场接口根据案例版本大致如下：

• 固体力学：设置弹性模量、泊松比、密度，以及 Biot 系数和有效应力相关内容。
• 脆性断裂接口：它需要指定相场长度尺度 l 和断裂韧性 Gc，或者直接指定裂缝能量释放率。
• Darcy 定律：在孔隙弹性耦合案例中描述压力扩散，需要给渗透率、孔隙度和流体黏度。
• PDE 接口：如果你不用内置脆性断裂，可以用系数型 PDE 手动建立相场演化方程。
• 瞬态研究：水力压裂是高瞬态过程，必须用瞬态求解器，边界载荷随时间逐步加载。

需要注意的是，不同 COMSOL 版本的内置实现细节略有差异，特别是应变能分解和退化函数的写法。我的经验是：先用低版本能跑通的模型再去升级，不要一上来就追求最新接口里的所有高级选项。

3. 案例一：单一裂缝扩展的完整建模过程

3.1 几何与材料参数：最容易栽跟头的其实是量纲

案例一的目标是复现一条初始裂缝在内部压力作用下沿最大主应力方向稳定扩展。很多人以为这个案例简单，其实初始几何和单位设置一旦出错，后面全白做。

我用的几何是二维矩形岩体模型，尺寸取 100 m × 100 m，中心布置一条长度为 10 m 的预制裂缝。这里的“预制裂缝”并不是真的要画一条开口缝隙，而是把相场变量 d 的初始值在那个位置设置为 1，其余位置设置为 0。材料参数可以按一种中等强度的砂岩。弹性模量取 20 GPa，泊松比 0.25，断裂能量释放率 Gc 取 100 N/m，相场长度尺度 l 取 0.5 m。单位必须全程统一：长度用 m，应力用 Pa，能量释放率用 J/m ²，这样压力载荷的单位才是 Pa。

很多同学习惯用 MPa 或者 mm 建模，COMSOL 内部是国际单位制，你用 mm 建模也可以，但所有导出材料参数、压力单位都必须跟着换算。我见过不止一次因为把 20 GPa 写成 20 MPa，结果裂缝完全不开，或者应力场小到损伤根本不演化的例子。

3.2 初始裂缝怎么给：常用两种方式对比

初始裂缝的设置有两种做法，效果差距很大。

第一种做法是纯几何凹陷：直接在几何里画一条很窄的长方形区域，把这个区域挖空，表示初始裂缝是真实开口。这个做法的优点是裂缝内部的流体压力可以直接加载到裂缝壁上，非常直观。缺点是尖端处容易产生严重的应力集中，而且当开口宽度接近网格尺寸时，网格质量迅速恶化，相场扩展时计算结果对网格极度敏感。

第二种做法是把初始裂缝区域设置为相场变量 d=1，不挖空几何，材料属性在损伤区域退化。这个做法更舒服，因为几何网格完全规则，应力集中被相场自然“抹平”了一部分。不过裂缝内的流体压力不能直接加在几何边界上，需要额外用一个体积载荷或者等效节点载荷，作用在损伤变量大于某个阈值的单元上。我在系列案例里默认使用第二种，因为后续扩展到多裂缝和非均质时最稳定，也是最容易复现的做法。

3.3 载荷加载与求解器设置：慢比快好

流体压力不能瞬间一步加到 10 MPa 或 20 MPa，相场模型对阶跃加载极其敏感。阶跃载荷相当于在零时刻给系统一个无穷大的加载率，非线性求解器很难收敛，即使收敛，损伤也会在奇异区域瞬间爆发，得不到稳定扩展路径。

我用的加载曲线是：

0 s 到 0.1 s，压力从 0 线性爬升到目标值，比如 15 MPa；然后是保持段。求解器选择全耦合，打开非线性控制器，最大迭代步数设置到 10 以上，初始时间步长取 1e-4 s 或者更小，后续交给自适应时间步。COMSOL 的相场模型本质上是一个数值非常刚的问题，d 的演化受长度尺度和时间尺度双重影响，宁可多花机时也要减小每一步的非线性残差。

如果计算不收敛，不要盲目减小总时间步，先看提示的是材料非线性问题还是全局牛顿法发散。通常把加载时间拉长、阻尼系数调大、或者在求解器设置里启用“辅助扫描”都能缓解。辅助扫描在这里特别管用：先把相场长度尺度设为一个偏大的值，比如 1 m，算出一个初步解作为初值，再扫描到 0.5 m，收敛性会明显改善。

3.4 结果怎么看：裂缝形态与压力响应

案例一最典型的结果输出有两样：损伤云图和注入端压力随时间的变化曲线。损伤云图可以直接看到裂缝是否沿着水平方向对称扩展，尖端是否出现树枝状分叉。正常情况下，在均匀应力场里，单裂缝应该从两端沿垂直于最小水平主应力的方向扩展，路径平直，裂缝宽度在注入点附近最大。

压力曲线同样能说明问题。如果压力持续上升，说明裂缝没有“打开”；如果压力在快速上升后突然下降，那通常是裂缝开始起裂扩展，体积空间变大，井底压力回落。这个趋势和现场压裂施工曲线非常相似，也是验证相场模型是否合理的重要指标。建议做两组敏感性测试：一组固定压力改变 Gc，看裂缝长度如何变化；另一组固定 Gc 改变长度尺度 l，看裂缝路径和端部形状如何变化。这两组测试能帮你确认参数选择落在“合理窗口”内。

4. 六套案例的设计逻辑与参数对照

4.1 案例清单总览

这一系列一共 6 个案例，我按照从简到繁、从单一物理场到多物理场耦合的顺序排布，每个案例都在前一个案例的基础上增加一个关键变量。下面把每一案例的核心内容列成了一张表，方便对照。

4.2 为什么案例要按这个顺序推进

很多人拿来一套案例就急着跑最后一个三维模型，我劝你先按顺序来。因为相场法水力压裂的“成功运行”是由多个环节共同保证的：初始条件、网格离散、非线性求解、流固耦合、历史变量保存。任何一环出问题，结果都可能完全荒谬。

案例一先把“单一裂缝能不能扩展”这个最基础的问题解决掉。案例二解决的是多裂缝之间通过应力场相互作用的问题，这在压裂设计中叫应力阴影，直接影响井间距和射孔簇位置。案例三引入材料非均质性，把相场法相比离散裂纹法的一个核心优势——不需要预设传播路径——充分展现出来。案例四才是真正意义上的水力压裂，流体压力、孔隙压力、裂缝扩展三个过程同时求解。案例五和案例六是工程场景的延伸，分别对应近井转向和三维缝网形态，也最接近实际施工关心的问题。

这种递进逻辑还有一个实际好处：当你改到一个新案例时，如果结果不对，你能立刻判断是哪一层引入的新机制出了问题。比如案例二跑挂了，可以想是间距问题还是载荷设置问题；案例四跑挂了，第一反应就是检查 Biot 系数、渗透率和时间尺度是否匹配，而不会去怀疑相场方程写错。

4.3 案例二和案例三：相场法最舒服的舒适区

如果你想快速体会相场法比离散裂纹法好用在哪里，直接做案例二和案例三就行。案例二里两条裂缝同时扩展，间距降到一定值后，内侧应力被释放，两条裂缝会向相反方向偏转，经典教科书里把这种现象叫裂缝排斥。用离散裂纹模型做这个案例，你必须不断判断尖端位置和偏转角，每扩展一步都要重画网格；相场法里只需要把两条预制损伤区放到几何里，然后让能量最小化自己去决定偏转路径，整个计算过程中网格拓扑完全不变。

案例三更爽。你在几何里把中间某个矩形区域的 Gc 调高或调低，相场法可以自动算出裂缝在高韧性层边缘是穿过还是止裂，不需要为界面做任何特殊处理。这种能力对研究页岩储层中不同岩性互层时的压裂高度控制非常有用，也是论文里最常出现的“三人行”类对比图来源。

5. 收敛失败、网格依赖与结果验证——实战中真正的拦路虎

5.1 相场长度尺度 l 和网格尺寸 h 怎么配合

等模型开始跑了，你才会发现：真正花时间的根本不是物理场搭建，而是解不出来、解出来不对、换个网格结果又变。这三个问题都指向同一个根源——相场长度尺度 l 与网格尺寸 h 的关系。

相场法的损伤带宽度由参数 l 控制，为了保证损伤带形态光滑，网格不能太粗，业内普遍经验是 h 不超过 l/2。也就是说，l=0.5 m 时，网格最大尺寸最好在 0.25 m 以下，再稳妥一点可以压到 0.1 m 左右。如果你的计算域是 100 m × 100 m，光这个网格就是几十万到上百万单元，这就是相场法计算量大的直接原因。

网格太粗时，裂缝路径会被拉成锯齿状，甚至出现斜向单元串扰，而不是一条平滑的曲线。我在案例一里试过用 1 m 的均匀网格去算 l=0.5 m，结果是裂缝扩展方向出现了大量 45 度折线，怎么看都像棋盘格而不是裂缝。如果你只是为了做方案对比，网格粗一点倒还能接受；但如果要出论文级别的损伤云图，就必须牺牲计算量把网格压下来。

5.2 为什么压缩区也会裂？多半是应变能分解没做好

一个很奇怪但常见的错误是：压力从裂缝内部向外推，结果裂缝不仅横向扩展，还整个岩体被“压裂”了，损伤云图一片红。这种情况十有八九是应变能分解没生效。

在脆性断裂相场模型中，开裂驱动力应该只来自拉伸主导的应变能分量，压缩主导的分量即使很大也不应该驱动损伤。COMSOL 内置接口里有对应的张量分解选项，但如果你自己写 PDE 方程，很容易在源项里直接采用总应变能密度。只要这么一写，压缩区的能量同样会驱动损伤，结果当然是算出一片虚假裂缝。检查方法很简单：单独显示总应变能密度和拉伸应变能密度，看损伤带的位置和哪个量对齐。如果和总应变能密度完全一致，那就要检查退化函数是否同时作用在了压缩分量上。

5.3 COMSOL“找不到解”时的排查顺序

我遇到过很多次同样的报错：求解器显示“找不到解”或者“At time t=... 最大迭代次数超限”。遇到这种报错，不要急着把时间步长改小一个数量级，更不要随便降低求解器精度，我从自己的实战中总结了一套排查顺序。

第一步，检查载荷是否阶跃。把压力载荷改成平滑斜坡或者用余弦过渡，很多发散问题立刻消失。第二步，检查相场变量的初始条件是否与几何位置一致，如果初始损伤区和边界载荷位置对不上，模型会先发生剧烈的局部应力调整，直接击穿 Newton 迭代。第三步，检查材料参数是否人为赋值失当，这里最常见的是 Gc 太小，导致损伤演化极快，一瞬间就整片开裂，时间步长根本追不上。第四步是打开非线性求解器诊断，查看每个物理场残差，如果只是压力场的残差发散，问题多半在 Darcy 流动的渗透率或时间尺度设置；如果是位移场残差发散，多半在弹性模量或边界条件。

5.4 验证三维模型的老办法：先归零再比较

案例六三维模型最容易出现的问题是“看起来像个缝，但数值上根本不知道对不对”。我的验证思路是先把问题退化到二维解析解能做的情况。把三维模型的厚度方向应力约束成平面应变条件，或者把初始裂缝简化成规则的椭圆，把注入速率调成恒定值，然后和经典水力压裂解析解对比裂缝宽度和半径随时间的变化。如果相对误差能控制在百分之十以内，再放宽到真正的三维地应力条件，心里就有底得多。

这里要格外提醒：三维相场模型中的损伤区域体积和单一长度尺度有关，直接用损伤云图去量“裂缝体积”会明显偏大。更合理的做法是提取 d 大于 0.9 的单元区域，再做一次体积积分，而不是把整个损伤带都算成裂缝。这个细节不处理，注入量与裂缝体积的匹配关系永远对不上。

6. 参考文献怎么组织，以及这套案例后续还能怎么改

6.1 文献按三个方向收，效率最高

这套案例之所以强调“附带参考文献”，是因为相场法水力压裂本身是一个典型的交叉课题：既要懂断裂力学，又要懂渗流力学，还得会有限元数值实现。文献如果只按关键词搜，很容易搜出一堆数学论文，读起来和生产完全脱节。我建议把参考文献分成三个方向：

第一个方向是相场断裂理论基础。这一类的核心是 Griffith 能量准则、正则化裂缝模型、应变能分解和不可逆条件，重点看它们怎么推导出相场控制方程，以及长度尺度 l 对断裂能耗散的影响。第二个方向是水力压裂的经典解析和半解析模型。这一类的代表是平面应变裂缝扩展解、径向裂缝扩展解、KGD 型和中心点源型模型，主要用来做相场结果的验证基准。第三个方向是 COMSOL 官方应用库和博客中与脆性断裂、水力压裂相关的模型文档。应用库里凡是用“Brittle Fracture”接口做的模型，都可以下载下来拆着看，比自己从零搭要快得多。

具体到案例一，我建议配套的文献至少包含：一篇相场断裂原始能量公式的综述、一篇基于应变能分解的数值实现文章、一篇水力压裂经典解析模型，以及 COMSOL 帮助文档里“脆性断裂”一节。有了这四篇，案例一背后的原理和校验基础就都齐了。后续案例再分别补充多裂缝干扰、非局部损伤和孔隙弹性耦合方向的文献。

6.2 这套系列下一步能怎么扩展

六套案例只是把主线趟通，离工程实用还有很长的路。在我看来，最有价值的扩展方向有三个。

第一个是裂缝内流体流动的细化。当前的模型用 Darcy 流或压力载荷近似，没有仔细区分流体在裂缝内的湍流/层流、支撑剂输送、滤失控制，这些对裂缝宽度和导流能力有直接决定作用。第二个是材料非线性和塑性。页岩和碳酸盐岩不是纯脆性，微裂隙、塑性屈服、各向异性都会改变损伤演化路径，可以在相场模型里加入塑性修正项，或者在弹性矩阵里加入各向异性参数。第三个是随机场和概率分析。真实地层参数是空间分布的，可以把弹性模量、断裂韧性或渗透率生成随机场，配合 COMSOL 的参数化扫描做蒙特卡洛样本，输出裂缝长度的概率分布，这样的结果在工程决策里会更有说服力。

我自己现在正在做的是把这个系列从“二维均匀模型”推向“三维非均质地层”，刚开始的几个三维模型不算稳定，但我发现只要坚持先做参数敏感性分析、再跑完整模型，绝大多数问题都能定位到具体环节。相场法做水力压裂确实门槛不低，但它的优势是逻辑统一，搭好一套框架后，后面换几何、换参数、换加载条件都非常顺手。最后再分享一个小技巧：把所有案例的 COMSOL 模型文件统一放到一个文件夹里，案例编号打头，模型文件命名带上,版本和网格数量，这样回溯以前的结果特别方便，尤其是当你回头发现某个参数需要重新跑的时候。