企业官网建设流程全解析

1. 为什么用相场法模拟压裂：传统离散裂缝方法的困境

做压裂模拟的人，不管是研究页岩气水力压裂还是实验室里的三点弯曲断裂，多少都会遇到同一个难题：裂缝是一条真正意义上的“不连续面”，模拟软件里怎么表达这条不连续面？早期主流的做法是离散裂缝模型，比如在ABAQUS里插入Cohesive单元，在Comsol里用裂纹尖端节点重划分，或者干脆用扩展有限元XFEM。这些方法各有各的适用场景，但真到复杂压裂工况下，问题就来了。

首先是Cohesive单元，它要求裂缝面沿预设路径走。真实岩石里的天然裂缝、层理面、随机节理可不管你预设没预设，裂纹一旦偏转，Cohesive单元就失效了。XFEM不需要预设路径，但需要显式追踪裂纹尖端位置，每个时间步都要判断裂纹尖端在哪里、往哪个方向扩展、扩展多远，三维情况下还要考虑裂纹面的形状变化，数值实现非常复杂，而且多裂纹交汇、分支的时候极容易崩。我用XFEM做过几组算例，感觉得出的结论很依赖“裂纹面追踪”的实现细节，换个追踪策略结果就变，心里很不踏实。

相场法（Phase Field Method，PFM）换了一个完全不同的思路：它不把裂缝当成一条“不连续面”，而是用一个标量场变量把裂纹“抹”成一个具有一定带宽的连续过渡区。这个思路最早来自脆性断裂的变分理论，后来被扩展到水力压裂、多场耦合，甚至疲劳断裂、动态断裂。在2000年前后由Francfort和Marigo提出变分框架，Bourdin等人用数值方法验证，Miehe等人又给出了便于有限元实现的形式，之后这套方法就成为断裂模拟领域的主流方案之一。

Comsol把相场法做成了一个相对完整的物理场模块，叫做“脆性断裂（Brittle Fracture）”接口，配合固体力学、流体流动、PDE接口可以搭出压裂模型。相比传统离散裂缝方法，它的核心优势有三个：第一，不需要追踪裂纹路径，裂纹扩展方向由能量最小化自动决定；第二，多条裂纹交汇、分叉、合并都是数值计算的自然结果，不需要额外处理断裂准则；第三，可以和流体流动、温度场、电场等多个物理场直接耦合，这对水力压裂、热致断裂这类多物理过程尤其友好。本文要讲的案例，正是围绕Comsol相场法模拟裂纹扩展的完整流程展开，从理论模型到参数设置再到收敛性调试，一次讲透。

2. 裂纹相场理论的数学与物理基础

2.1 裂纹怎么“铺”进连续介质里

相场法的出发点是一个叫做相场变量（Phase Field Variable）的量，通常用φ表示。φ的物理意义是材料局部损伤程度：φ=0代表材料完好无损，φ=1代表材料完全断裂。在真实的裂缝面两侧，φ从一个很小的值连续过渡到接近1的值，这个过渡区域的宽度由参数l0控制。可以把这个过渡带理解为“裂缝的模糊化表示”——它不是一条几何上的线，而是分布在空间里的一条带子。你把l0取小了，裂纹带就变窄，更接近真实裂缝；但l0太小会导致网格尺寸要求极高，计算量指数级上升，所以实际上l0是精度和计算成本之间的一个平衡点。

有了φ之后，原有的弹性应变能就要被“惩罚”掉。常用做法是在应变能密度前面乘一个退化函数g(φ)=(1−φ)2+k，其中k是一个很小的正数，用于避免完全断裂后数值刚度为0带来的求解困难。这个退化函数的直观含义是：材料越接近断裂，φ越接近1，那么它能存储的弹性应变能就越小，相当于刚度逐渐退化。这个处理非常巧妙地绕开了“裂缝面如何几何更新”的难题——刚度退化自动完成了裂纹对结构刚度的影响。

裂纹扩展的驱动力从哪来？来自能量释放。在脆性断裂理论里，裂纹伸长的条件是系统总能量下降，其核心准则是Griffith准则：只有当裂纹尖端附近的能量释放率达到或超过材料的临界能量释放率Gc时，裂纹才会扩展。相场法把这个准则变成了一种“自动演化”的形式——不需要单独判断裂纹是否达到扩展条件，因为能量最小化过程本身就驱动着φ场演化。这就是相场法看似“智能”的根本原因。

2.2 能量控制方程与两个偏微分方程的耦合

完整的相场断裂模型至少包含两个耦合的偏微分方程。

第一个是弹性力学方程，描述位移场u在含损伤材料中的平衡：

∇·σ + F = 0

但这里的应力σ不再是纯弹性应力，而是通过退化函数和相场变量修正后的退化应力，即σ = (1−φ)2+ k) σ0，其中σ0是未损伤材料的线弹性应力。这个方程和标准固体力学方程在形式上非常接近，只是在材料本构里多了一个空间变化的退化系数。

第二个是相场演化方程，是一个Allen-Cahn/Ginzburg-Landau类型的方程，一般写成如下形式：

(φ/l0) − l0∇²φ = 2(1−φ)H（等等，具体形式取决于不同文献）

Miehe给出的普及版方程为：

1/M φ点 − (1−φ)H + l0∇²φ − φ/l0 = 0

其中M是相场迁移率，H是历史应变场变量（History Field），它取整个加载过程中驱动应变能密度的最大值。引入历史场变量H是Miehe工作的一大亮点，它保证了裂纹扩展的不可逆性——一旦材料在某处损伤，即使载荷降低，损伤也不会自动愈合。这个不可逆性在物理上非常重要，因为真实裂纹不会因为卸载就消失。

为了精确建立这个退化函数驱动相位场的公式，可以让我更详细地写出已在文献中广泛用于脆性断裂的Miehe公式：

H = max(s∈[0,t]) Ψ0(+)

其中Ψ0(+)是拉伸部分应变能密度。当耦合流体时，压力产生的体积力或孔压应力都会通过这个应力张量影响H，使得裂缝路径呈现出压裂特有的“沿最大主应力方向扩展”特征。

2.3 长度尺度参数与材料韧性

不可回避的参数是l0（长度尺度参数）和Gc（临界能量释放率）。在相场模型中，l0形如裂纹带宽；Gc形如驱动裂纹传播的能耗。有两个关系式需要格外注意。

首先是强度和l0的关系。相场模型的峰值应力通常满足关系σc ∝ √(EGc/l0)（不同退化函数形式下系数略有差异）。也就是说，如果你设定了l0过大，模型的承载强度会被人为降低，模拟出的断裂峰值力会显著小于实验值。因此l0的取法一般与网格尺寸h绑定，通常取l0=(2~4)h。另一条准则是解析解验证：固定Gc和弹性模量，把l0从0.5mm、1mm、2mm、4mm拉一遍，算出的载荷-位移曲线中峰值应力变化如果超过5%，说明l0太大，要么减小l0，要么换更细的网格。这条验证法是我实测有效的。

其次是Gc的取值。Gc的本质是单位裂缝面积扩展时耗散的能量。对于岩石类材料，Gc通常在几十到几百J/m²的范围，具体数值需要根据压痕实验、三点弯曲实验或文献检索确定。实在拿不到实验数据时，也可以用断裂韧性KIC换算：Gc=KIC²/E'（E'为平面应变弹性模量E/(1−ν²)）。这个换算关系在水力压裂参数标定中非常常用，因为很多岩土的KIC数据比Gc数据好查。

3. 完整案例的Comsol实现

3.1 案例工况设定

这个案例模拟的是一个典型的水力压裂过程。尺寸采用实验室常见的平板试件：长400mm、宽200mm，厚度按平面应变处理。试件中心预置一条长度为40mm的初始裂纹（用初始相场φ=1的窄带来实现）。加载方式有两种选择：一是压裂液从裂纹中心注入，压力随时间逐渐升高；二是右侧施加固定位移载荷模拟劈裂。本文采用前者，更贴近水力压裂场景。

材料参数取一组典型砂岩数据：弹性模量E=28GPa，泊松比ν=0.22，抗拉强度3.5MPa，断裂韧性KIC=1.2MPa·m^(1/2)，由此换算出临界能量释放率Gc大约60J/m²。这个数据组和真实砂岩比较接近，算出的裂缝形态容易被后续实验验证。

这里要提醒一下：初始裂纹不要用一个纳米级的几何线切割来表示，那样网格剖分很麻烦。简便做法是把初始裂纹区域的材料相场初值直接设为1，其余地方设为0。实现方式是在“初始值”设置里给φ赋一个与坐标相关的条件，比如如果x<0.02且|y|<0.001，则φ=1。这个“初始损伤带”完全可以充当初始裂纹，而且带来的收敛麻烦比几何切割少得多。

3.2 几何、材料与参数表

在Comsol中新建二维模型，几何就是一个400×200的矩形，初始裂纹区域在模型中心靠左位置。材料节点添加弹性模量与泊松比，注意相场断裂接口需要的是“损伤材料的弹性材料”，在“脆性断裂”接口内部自带的材料节点里填写即可。

我习惯把所有控制参数集中到一个全局参数表里，后面做参数化扫描时非常省事：

3.3 物理场接口配置

Comsol 6.4里做这个案例，推荐直接在“结构力学”模块下添加“脆性断裂（brittle）”接口。这个接口会自动创建一组耦合的方程组：基础的固体力学方程、相场演化方程，以及必要的多物理场耦合节点。你不需要手动去写PDE，对多数工程应用来说内置接口已经足够。

不过有两点要手动处理。

第一是把“裂缝扩展”相关的求解变量打开。进入“脆性断裂”接口的“相场设置”子节点，确保“历史应变场”选项是启用的。历史应变场是保证裂纹不可逆扩展的钥匙，很多用户算着算着发现裂纹扩展后又缩回去了，十有八九是这里没设对。

第二是加载方式。水力压裂不是单纯在边界施加力，而是注入流体推动裂纹。严格的三维理论需要土力学中裂纹的流固耦合，但作为二维近似，你可以直接在裂纹区域边界上施加逐渐增大的“边界载荷”模拟注入压力，或者更精细一些，用Darcy定律接口算出压力场再耦合到固体力学。对于入门案例，用“边界载荷+斜坡升压”就够了。具体做法是：在裂纹带的内侧边界上施加压力p(t)=p0×t/t_pump，也就是从0在5秒内线性升到2MPa。这个斜坡函数让系统有一个准静态加载过程，对收敛友好很多。

3.4 边界条件与初始裂纹

边界条件设置需要注意夹持方式。矩形试件的左下角和右下角分别约束x方向和y方向位移，模拟实验系统的支撑。上边和右边自由。如果你要做单轴压缩或围压加载，再在上面加相应的分布载荷。

初始裂纹用“初始损伤带”的方式实现。在相场的“初始值”里写一个表达式，比如：0.5*(1-tanh((sqrt((x-0.15)^2+y^2)-0.02)/0.0005))，这个表达式的意思是：在(0.15,0)点为中心、半径20mm的圆形区域内，φ初值接近1（裂纹），区域外φ接近0（完整材料）。用tanh做光滑过渡的好处是避免了尖锐跳变引起的初始收敛困难。注意初始裂纹不要直接开在模型正中心，略微偏左一点，给裂缝扩展留出更大的自由空间，这样模拟出来的偏转路径更真实。

3.5 网格剖分策略

网格是相场法最考验耐心的环节。l0=1mm的前提下，裂纹带内至少要有3~4层单元，那么裂纹带附近网格尺寸应该在0.3mm左右。但整个400×200mm的模型不可能都用0.3mm网格，否则单元数量几十万起，算得昏天黑地。

推荐做法：先用“用户控制网格”划一个全局较粗的底网（4mm），再用“尺寸”节点限定裂纹扩展可能经过的区域，也就是模型中间高度±30mm的一条长带，尺寸改为0.4mm。这条细网格带就像一条“裂纹跑道”，裂缝扩展主要被限制在这条带内。当然，真实裂缝可能偏转出跑道，所以跑道宽度不要留太窄，留出上下各30mm是比较均衡的选择。

单元类型方面，用拉格朗日二次单元计算精度更高但更费钱，一次单元在三节点三角形下容易过度刚硬。我建议用三角形二次单元，同时开启“细化”选项。网格数量控制在2~5万之间，个人电脑算起来不会太吃力。

还有一个网格上的坑：不要在初始裂纹尖端使用极小网格而背景网格又很粗，网格尺寸从0.3mm直接跳到4mm会造成尖端应力场失真，裂缝启动方向可能被网格不对称性“带偏”。网格尺寸过渡系数最好控制在0.2~0.3，即相邻区域网格尺寸变化不超过3倍。

4. 移动网格与压裂流体注入的耦合处理

4.1 为什么需要用移动网格

压裂模拟和普通力学断裂模拟最大的区别在于：裂缝张开后，内部空间被流体充满，流体压力又反过来作用在裂缝面上，形成流固耦合。如果按固定网格算，裂缝面张开后节点严重错动，网格质量会迅速恶化，尤其是窄长裂缝的尖端区域，三角形单元会被拉伸成非常扁平的形状，计算精度和收敛性同时崩掉。

解决思路是在Comsol里启用“移动网格（Moving Mesh）”功能，这里的核心思想是：把材料框架坐标换算为物质框架，通过重新分配网格节点位置来适应裂缝张开带来的大变形。对于压裂问题，移动网格更像一个“网格协调器”——它不让网格完全消失，但能把畸变控制在可接受范围。顺带一提，如果你的模型涉及COMSOL的压电效应、热-力耦合等其他物理过程，移动网格同样适用于这些物理场的重新映射，这是Comsol 6.4版本很实用的能力。

4.2 移动网格设置方法

在“定义”节点下添加“移动网格”接口，然后设置“变形域”为整个矩形，并让几何变形由位移场u驱动。具体实现是：在变形域设置中选择“由实体位移驱动”，然后勾选你固体力学接口计算出的位移变量。这个操作相当于告诉Comsol：网格的每个节点随身跟着固体变形走。

同时需要设定“固定（网格）边界”——通常是模型的底部边界，作为网格锚点。如果你不设置固定边界，整个网格会随着刚体位移漂移，算到最后结果根本没法看。

移动网格配合相场法需要特别留意“网格质量检查”功能。每迭代几步，查看一下最小单元质量。如果发现有负质量或接近零质量的单元，就需要回退时间步或者改变网格约束。Comsol的自动时间步进通常会在网格崩溃前减小步长，但不代表它能完全兜底，模型大了依然有可能爆掉。

4.3 求解器与时间步控制

固体力学+相场+移动网格是非线性程度很高的多物理场耦合，直接上默认求解器很容易出现“第一步就发散”。我的经验是用以下配置：

• 求解器类型用“PARDISO”直接求解器，这个对多物理场耦合的鲁棒性比迭代求解器好，虽然内存消耗高一些，但值得。
• 开启“自动”时间步选择，初始步长设2×10⁻³秒，最大步长0.1秒。相场法对时间步非常敏感，步长太大时，量场在一步内变化剧烈，牛顿迭代极易发散。
• 在“瞬态求解器设置”里，把“回落”和“一致性”选项全部打开，并设置最大迭代次数为10。如果10次牛顿迭代仍不收敛，自动减小时间步重算。这套配置的代价是总求解时间拉长，但换来的是稳定性，我至今没遇到过直接跳不出来的情况。

还有一个经验：如果加载过程很短或压力很高，惯性效应不可忽略，记得在固体力学接口中开启瞬态结构分析（含惯性项），而不是默认的“准静态”模式。准静态模式下压力陡增会导致不合理的“应力突然穿透”现象，在压裂尖端尤其明显。

5. 常见问题与收敛性排查实录

5.1 裂纹不扩展，只停在初始损伤区里

这是最常见的坑。排查顺序我是这样做的：先看H历史场是否更新正确，如果H始终为0，驱动项就为0，相场根本不会演化。多数情况下历史场没传导到相场方程，疑似耦合节点漏选。

再看Gc是否过大。Gc=60J/m²在砂岩里不算离谱，但如果你给的是钢材的Gc（大几十万），怎么加载也不会裂。用一条简单准则判断：当裂纹尖端的能量释放率G达到Gc，裂缝才会扩展；如果模拟压力升到3倍也扩展不了，十有八九是Gc数量级不对。

最后检查网格。l0=1mm但网格0.8mm时，裂纹带宽内只有1~2个单元，相场插值不充分，驱动项会被严重低估。请保证l0/h≥3，这个比值是相场法数值稳定性的经验底线。

5.2 相场在初始裂纹以外突然“糊”成一片

这通常是因为“拉伸-压缩能量分解”没有正确启用。相场断裂理论里有个细节：只有拉伸应变能用于驱动裂纹扩展，压缩应变能可以允许裂纹闭合，不驱动开裂。如果在Comsol里没有区分拉伸和压缩，那么压缩应力同样会驱动裂纹“张裂”，这明显违反物理直觉，表现就是模型受压的位置也发生损伤，裂纹像流感一样传遍全场。

解决方法是：在脆性断裂接口的“相场设置”里，把能量分解方式选为“Volumetric-Deviatoric”或“Spectral”分解。光谱分解计算成本更高但精度更好，小模型直接用Spectral没问题。

5.3 压力无法维持，应力提前松弛

模拟注压时，有时裂缝还没扩展，但压力载荷对应的位移已过大，系统刚度因退化函数降到极低，导致压力无法维持。这里大概率是退化函数中最小残余刚度k设得太小了。k通常取1×10⁻⁶~1×10⁻⁸，取小了计算更精确但刚度矩阵更容易奇异。遇到这个问题，把k从1e-8提到1e-6，往往应力分布立刻恢复物理合理性。

还有一个可能性是时间步过大导致在“加载步”内刚度退化跳变过大。把最大时间步降到0.05秒，问题也能缓解。

5.4 收敛失败时怎么“救场”

碰到多次收敛失败，我有一套“救场四级”操作：

• 第一级：减小最大时间步，关闭“高精度”几何阶次。
• 第二级：把退化函数参数k从1e-10逐步提到1e-6。
• 第三级：把相场演化方程改成显式解耦——先用固体力学算出位移，再把位移冻结，单独算一步相场，然后交替更新。这个“交替求解”牺牲了一些严格性，但超不出可接受范围，在很多复杂多物理模型中是标准操作。
• 第四级：换PARDISO求解器、降低非线性容差到1e-3，并减少最大迭代次数到6，让程序更快地“缩步”。这一招本质是“让求解器早认输，缩小步伐重新打”，通常能避免灾难性的发散。

如果四级都没用，那基本是模型本身设置有问题，回到5.1~5.3逐步排查，不要硬加时间步。

6. 参考文献、验证方法与进一步扩展

6.1 核心文献清单

要严谨地做相场压裂模拟，参考文献不能省。这里列一个最小清单，够入门和中期使用：

• Francfort G A, Marigo J J. Revisiting brittle fracture as an energy minimization problem[J]. Journal of the Mechanics and Physics of Solids, 1998, 46(8): 1319-1342. 这是相场断裂的奠基文献，变分框架的起源。
• Bourdin B, Francfort G A, Marigo J J. Numerical experiments in revisited brittle fracture[J]. Journal of the Mechanics and Physics of Solids, 2000, 48(4): 797-826. 第一个数值实现，很多数值细节可以追溯到这里。
• Miehe C, Hofacker M, Welschinger F. A phase field model for rate-independent crack propagation: Robust algorithmic implementation based on operator splits[J]. Computer Methods in Applied Mechanics and Engineering, 2010, 199(45-48): 2765-2778. 工程实现的分水岭，历史场变量、应变能分解都是这篇奠定的。
• Miehe C, Mauthe S. Phase field modeling of fracture in multi-physics problems. Part II: Coupled brittle-to-ductile failure criteria and propagation of interfaces[J]. Computer Methods in Applied Mechanics and Engineering, 2015, 296: 286-325. 多物理耦合扩展，对流体-力学耦合有详细推导。
• Wilson Z A, Landis C M. Phase-field modeling of hydraulic fracture[J]. Journal of Applied Mechanics, 2016, 83(4): 041001. 这篇是把相场法用到水力压裂的比较早的论文，流固耦合处理有启发。

有精力的朋友建议再去找Miehe那篇“operator split”论文附录里的离散化公式，照着它才能对Comsol内置接口的每项设置心知肚明。

6.2 验证方法与参考解

模拟做得再漂亮，不验证就没有说服力。我推荐的验证路线有三条：

• 弹性力学标准解验证：拿一个单边缺口板，施加拉伸载荷，算出峰值载荷，与线弹性断裂力学的解析解对比。当Gc、E、ν、初始裂纹长度代入公式，得到的峰值拉力相差不超过5%时，说明模型的基础参数和边界条件队形没问题。Kanninen的《Advanced Fracture Mechanics》里有标准算例。
• 三点弯曲或紧凑拉伸实验数据验证：如果你有实验室条件，直接做一组成岩试件的三点弯曲实验，把实验载荷-位移曲线和模拟曲线叠在一起对比，看峰值和软化段是否重合。这套验证能有效检验参数标定是否合理。
• 相场不可逆性验证：卸载后检查φ场是否保持不变。如果φ场在卸载后减小了（说明裂纹愈合了），历史场设置一定有问题。

很多投稿人拿“裂纹形态与实验照片一致”来证明模型可靠，但这更像定性验证。真正的可靠性验证必须落到载荷-位移曲线、峰值力、应变场等定量指标上。Comsol的后处理里导出历史应变分布和裂纹形态都比较方便，建议把每一步的φ场云图和实验测量对照着看。

6.3 从入门到复数裂纹扩展的进阶方向

跑通单条裂纹案例后，可以往里加的东西很多。可以试离散裂缝网络与相场混合法，把天然裂缝网络预设成多个初始损伤带，模拟压裂液沿随机裂缝网络扩展的过程，观察分支和交汇。可以把流动部分从边界压力换成达西流场/布里奥流场耦合，模拟孔隙压力传播和裂纹扩展之间的相互影响。可以做参数化扫描，研究注入速率、流体黏度、压裂液温度对裂纹形态的影响，这组算例对文献参考很有价值。

如果想继续往深处凿，建议研究一下“多物理场相场压裂”的三维版本。三维比二维多了体积约束问题，网格量翻几倍，但物理上能把裂缝的平面外扩展和偏转都算出来。方法框架不变，只是对机器内存和耐心要求高了不少。Comsol 6.4在三维相场模块上的性能优化做得越来越好了，值得尝试。

我个人做下来最大的体会是：相场法看起来“高大上”，但它的参数标定和网格质量敏感度都非常高，很多“没跑通”的案例不是理论问题，而是对l0和网格尺寸之间的比例认识不够。初学者最容易犯的错误是一上来就猛加密网格，算几分钟就崩，然后归咎于模型不对。实际经验告诉我，先把粗网格跑通，再逐级加密，观察结果变化趋势，反而能更快地找到可靠参数。相场法模拟是一个“精度与成本的游戏”，不是说网格越细越好——你要的是一条物理合理的裂纹路径，而不是无穷多的网格单元。先把这个平衡掌握住，再去追求复杂工况，会顺很多。