企业官网建设流程全解析

在 COMSOL 里做裂纹相场法模拟，尤其是压裂相关的案例，最容易被卡住的往往不是软件操作本身，而是“参考文献里的方程到底该怎么落到模型里”。我从离散裂纹方法转过来时，对着能量泛函愣了大半天，后来把 Griffith 断裂理论、Miehe 等人的相场演化方程和 COMSOL 的系数型 PDE 逐项对上号，案例才算真正跑通。这篇博文就把整条路径摊开讲：相场法为什么会成为压裂模拟的常用选择、控制方程的核心逻辑、COMSOL 建模实操步骤、长度尺度和网格该怎么配，以及压裂相场方向的参考文献怎么选、怎么看。

1. 为什么压裂案例要选相场法：网格依赖问题给我的教训

1.1 传统离散裂纹方法在扩展问题上的短板

几年前我做水力压裂模拟，用的还是离散裂缝模型。那时候的思路很直观：把裂缝当作一条几何线，裂缝尖端接一个应力奇异性判定，满足条件就向前扩展一小段。这个做法在预置缝、单裂缝、路径明确的情况下很有效，可一旦涉及多条裂缝交汇、裂缝偏转、或者注入压力导致裂缝向阻力小的方向拐弯，问题就来了：裂缝扩展路径极度依赖网格划分方向。网格是斜的，裂缝就往斜着走；网格是对称的，裂缝又可能同时往两个方向分叉。这种“网格诱导的虚假裂纹路径”在学术评审里几乎一眼就能看出来，做工程更是没法接受。

简单解释一下背后的原因：离散方法把裂缝实际当作体内边界，尖端附近的应力场本身有奇异性，而有限元网格又是有限分辨率的，尖端附近真实应力分布根本算不准，只能依赖人为准则判断“下一步裂不裂、往哪里裂”。判断结果必然被离散误差牵着鼻子走。

1.2 相场法为什么能绕开这个坑

相场法的思路完全不同。它不把裂缝当成一条线，而是用一个连续变化的标量场 d(x,t) 去刻画材料损伤程度：d=0 表示完整材料，d=1 表示完全断裂。裂缝被“涂抹”成一个宽度有限的损伤带，虽然带内存在一个特征宽度参数 l，但裂缝的路径不再由网格方向决定，而由力学驱动项和断裂阻力共同决定，属于偏微分方程的固有解。网格只需要在这个损伤带内足够细就行，路径本身是自由演化的。

更关键的是，相场法给出了统一的本构框架：一个能量泛函同时包含弹性应变能、断裂表面能和外部力做功，裂缝的扩展等价于整个系统能量最小化。这样就不需要人为判断“该不该裂”，只要驱动力超过材料断裂韧性，损伤场就会自发演化。对压裂这种需要同时考虑注入压力、孔隙流体压力、地应力相互作用的强耦合问题，相场法在理论逻辑上非常干净。

1.3 压裂场景对相场法特别友好的三个理由

• 路径自由：相场法的裂缝路径由能量最小化决定，裂缝偏转、分叉、交叉都可以自然出现，不需要预设拓扑。
• 耦合顺畅：损伤场可以直接参与应力张量的退化，也能控制裂缝区域的渗透率，流体注入和裂缝扩展可以统一在同一套有限元框架里求解。
• 工程可复现：COMSOL 的系数型 PDE 接口、固体力学接口、达西流动接口可以拼出完整的压裂模型，计算结果稳定性比传统重网格方法好得多。

我后来把预置缝的算例从离散方法换成相场法，同样的几何、同样的材料参数，裂缝扩展路径不再跟着网格走了，这算是让我彻底认可了这套方法的第一个实际案例。

2. 相场法的核心方程：把 Griffith 断裂能量写进连续介质

2.1 从 Griffith 准则到裂缝表面能

相场法的起点是 Griffith 断裂准则：裂缝是否扩展，取决于裂缝尖端附近的能量释放率 G 是否达到材料临界值 Gc。这个准则物理含义很好理解，但难点在于裂缝尖端的位置和路径是未知的，且格里菲斯准则在裂缝尖端处存在应力奇异性，数值上很难直接离散。

相场法的聪明之处在于：把裂缝的几何不连续问题换成扩散边界的连续问题。假设裂缝表面能可以写成卷积形式，用损伤场 d 的梯度来近似裂缝面积：

Γ(d) = ∫_Ω [ d²/(2l) + (l/2) |∇d|² ] dΩ

这个积分里，第一项让损伤状态停留在 1 时付出能量代价，第二项让损伤区域局限在宽度约 l 的范围内。当 l 趋近于 0，这个近似就恢复到尖锐裂缝的表面能。实际计算中 l 取有限值，但没有必要无限小，因为它同时也承担了正则化参数的角色。

2.2 总能量泛函和耦合方程

把弹性能和断裂能写在一起，系统的总能量泛函是：

Ψ = ∫_Ω [ g(d) ψ₀⁺(ε) + ψ₀⁻(ε) ] dΩ + ∫_Ω Gc [ d²/(2l) + (l/2) |∇d|² ] dΩ - ∫_Ω p α ∇·u dΩ + ∫_Ω (1/2 M) ∇p · ∇p dΩ ...

这里 g(d) 是刚度退化函数，最常用的是二次型：

g(d) = (1-d)² + k

k 是残余刚度的极小值，目的只是保证数值稳定性，避免完全零刚度导致求解器奇异。应变能还要做拉伸/压缩分裂：ψ₀⁺ 对应拉伸驱动的损伤，ψ₀⁻ 对应压缩，因为材料在受压时不应该产生断裂。如果不做分裂，两个主应力都受损伤影响，会产生虚假的裂纹压闭合问题，这在压裂模拟里尤其危险——裂缝面可能被压碎，流体的流动路径也就不对了。

损伤变量的演化方程，常用形式是：

η (∂d/∂t) = Gc (l ∇²d - d/l) + 2(1-d) H

H 是历史驱动力，来自过去所有时刻拉伸应变能的最大值：

H(x,t) = max_{τ≤t} ψ₀⁺(ε(x,τ))

历史最大值的引入是为了保证裂缝扩展的不可逆性：损伤一旦产生就不能愈合。这一点在压裂注入压力卸载阶段尤其重要，否则数值上很容易出现“裂缝自动闭合、损伤自动消失”的虚假现象。

2.3 压裂案例里的流固耦合补充

如果只做纯力学断裂，模型还比较基础。真正的水力压裂案例还需要把液体压力和裂缝扩展耦合起来。常见简化做法是在能量泛函中加入孔隙压力项 p，有效应力写成：

σ_eff = g(d) σ₀ - α p I

其中 σ₀ 是排水条件下的有效应力，α 是 Biot 系数。液体流动在未损伤区按达西定律，在损伤区则要额外加上一条高渗透的裂缝通道，通常把渗透率写成损伤场的函数：

k(d) = k₀ [ (1-d)² + κ ]

κ 是裂缝残余渗透率，用来避免裂缝完全闭合后流动矩阵奇异。这个式子物理上对应一个朴素事实：裂缝越宽、损伤越严重，流体越容易通过。实际建模里，也可以不显式引入流体方程，而是把注入压力作为边界载荷直接加在预置裂缝面上，对整个模型做所谓“半耦合”简化。这种半耦合方案对理解相场法本构和网格行为非常友好，也是很多 COMSOL 教程案例的默认路径。

3. COMSOL 建模实操：从几何到系数型 PDE 的完整搭法

3.1 模型向导和物理接口选择

打开 COMSOL，建议先选二维模型做压裂案例。二维模型计算成本低，可以快速做参数扫描和网格敏感性分析，把原理跑通了再升级到三维也不迟。

物理接口选择上，我的推荐组合是：

在模型向导里逐个添加这些接口，再选择瞬态研究。如果 COMSOL 版本支持“系数型 PDE（coefficient form PDE）”，就用它；找不到的话也可以用“一般型 PDE”，但系数型 PDE 的接口排布更直观。

3.2 几何建模和裂缝初始条件的处理

几何方面，我习惯建一块矩形域代表岩体，比如长 40 m、宽 30 m，中间预留一段“初始裂缝”位置。预留方式有两种：

• 几何凹槽法：在裂缝位置画一条细缝，真实地掏空材料。优点是最直观，缺点是后续网格和压力边界要单独处理，裂缝尖端附近网格容易畸形。
• 材料初始损伤法：几何不掏洞，而是在裂缝位置把损伤场初值设为 1，或接近 1。例如 d=0.999，这样初始状态就存在一条弱化带，压力一上来就会沿着这个区域扩展。这种方法更贴合相场法的连续场思想，也不用修改几何拓扑。

我建议新手采用第二种。具体操作是在初始条件里写一个表达式，比如：

d0 = 0.999 * exp(-(x^2 / (0.5^2))) * (y < 3)

这会在 y=0 附近沿 x 方向生成一条高斯型损伤带，模拟一条沿水平方向约 1 m 左右的初始裂缝。注意不要把整个区域都设为损伤，否则模型一开始就失去承载能力。

3.3 全局参数和变量定义

进入“全局定义 → 参数”，建议把参数全部列出来，方便后面扫描。下面是一组可用的参考值：

然后在“组件 → 定义 → 变量”里写：

• 退化函数：g_damage = (1-d)^2 + k
• 总应变能密度：psi_0 = 0.5 * solid.SE * (solid.J) ...，不过这里更稳妥的做法是直接调用 COMSOL 固体力学接口提供的应变能密度变量，再自己定义拉伸部分的投影。
• 损伤历史驱动力：H_field，先初始为 0，后续用历史变量更新。

应变能分裂这一步最容易写错。如果不想用复杂的方向性分裂，可以用体积/偏量分裂近似：把应变能分解为体积应变能和偏量应变能，只有体积拉伸部分驱动损伤。这样写起来简单，压裂场景下精度足够。

3.4 系数型 PDE 的系数怎么填

这是整个建模最核心的一步。损伤演化方程：

η (∂d/∂t) = Gc (l ∇²d - d/l) + 2(1-d) H

在 COMSOL 系数型 PDE 里，因变量设为 d，时间导数项系数 e、扩散项系数 c、吸收项系数 a、源项 f 分别填：

按 COMSOL 的约定，方程形式是：

e ∂d/∂t + ∇·(c∇d) + a d = f

所以把 c 设为-Gc*l，a 设为Gc/l，方程展开正好是：

eta ∂d/∂t - Gc*l ∇²d + Gc/l d = 2(1-d) H

移一下项就得到标准相场演化方程。注意这里的负号经常有人填反，填反了扩散项方向错误，损伤场会莫名其妙地从边界“漏”进整个域，算出来的裂缝比实际宽好几倍。

H_field 的处理我放在 3.6 节详细说，这里先假设它在材料非局部变量里已经定义好。

3.5 固体力学接口的耦合设置

在固体力学接口中，应力张量要乘上退化函数 g_damage。具体有两种做法：

• 修改弹性矩阵：在“线弹性材料”中启用“损伤”选项，把损伤变量关联到 d。COMSOL 6.x 的线弹性材料节点可以直接定义损伤刚度，这个路径最省事。
• 手动改写应力张量：把应力分量替换为(1-d)^2 * 原应力 + k * 原应力，但这种方法需要逐项修改，公式冗长且容易漏掉泊松效应。

推荐第一种做法，ROI 最高。COMSOL 6.4 的“损伤”子节点可以自由定义退化函数，最省心的设置是把退化函数表达式写成：

g_fun = (1-d)^2 + k

同时把“裂缝应力分裂”选为“应变能分裂”。如果版本太老没有这个节点，就退回手动改写。

3.6 历史变量 H 的不可逆更新

不可逆性怎么实现，是一个容易忽略的隐蔽点。我之前试过直接不用历史变量，只用当前时刻的应变能驱动损伤。结果在注入压力波动时，损伤场竟然跟着回退，裂缝自动“愈合”了。这违背物理常识，而且后续计算完全不可信。

简单可行的做法是在 COMSOL 里加一个额外的因变量 H，用全局常微分方程更新：

H_new = max(psi_plus, H_old)

在系数型 PDE 的 f 项里写2*(1-d)*max(psi_plus, H_old)也可以实现类似效果。但要注意这种表达式是强非线性的，求解器可能需要更紧的容差。更稳妥的工程做法是使用 COMSOL 离散化的“事件”功能，在每一步求解后更新 H 场，然后代入下一步。

如果这些做法都嫌复杂，还有一个折中方案：加载时采用严格的单调递增压力曲线，并在损伤演化中加入足够大的粘性 eta。粘性项的稳定作用让它对历史依赖没那么敏感，虽然理论严谨性打折扣，但作为入门案例跑通流程完全够用。

4. 裂纹宽度、网格密度和求解器收敛：实际调试记录

4.1 长度尺度 l 的选择为什么决定了网格成本

相场法有一个绕不开的配对关系：长度尺度 l 决定损伤带宽，而网格尺寸 h 必须能分辨这个带。经验法则是：

h ≤ l / 2

也就是说，损伤带内至少要划分 2 到 3 层网格。这个约束带来的代价很现实：l 取得越小，模型越接近真实尖锐裂缝，但网格数量爆炸式增长。我见过不少新手想把 l 压到毫米级，模型直接卡死在网格剖分阶段。

实操建议是：先把 l 取到等于或略大于最小特征尺寸的 1/10 到 1/20，跑通后再逐步缩小 l，同时观察结果变化。如果 l 缩小后裂纹路径和断裂载荷变化很小，说明这个尺度已经够收敛；如果差异明显，说明先前的结果是“正则化依赖”的，还不能算物理结果。

具体到压裂案例，我在一个 30 m 宽的岩体模型里取 l=0.5 m，损伤带整体宽度大约 2~3 m，网格在裂缝路径区域加密到 0.2~0.25 m，其余区域用 1~2 m 渐变网格。整体自由度控制在 5 万以内，个人电脑几分钟到十几分钟能算完一个工况，这个量级很适合做参数扫描。

4.2 求解器配置和时间步长控制

相场方程的强非线性主要体现在源项2(1-d)H和刚度退化 g(d) 上。默认的全耦合求解器有时候会非常吃力，出现“找不到解”或者“迭代发散”的提示。

我的实用配置是：

• 启用阻尼牛顿法，初始阻尼因子设为 0.5。
• 时间步选自由步长，但设定最大步长，通常取 0.1 到 1 秒，取决于注入时间尺度。
• 相对容差设 1e-3，绝对容差根据损伤变量和位移的量级调整到 1e-4。
• 如果全耦合发散，就改用分离法：先固定 d 求位移，再固定位移求 d，反复交替迭代。这种做法看似多一个循环，但每个子问题都更接近线性，收敛稳定性显著提升。

以我的压裂案例为例，注入压力从 0 线性加载到 5 MPa，加载时间 10 秒。直接全耦合求解时，在损伤剧烈扩展、裂缝尖端起裂的那一刻总是报错。改成分离法，并且把起裂时刻的最大时间步长压到 0.02 s，就能跨过去。起裂阶段的时间步要小，这个经验非常重要。

4.3 几个典型的错误现象和对应措施

我把调试中遇到的常见现象整理成一张表，方便对照排查：

4.4 一个必做的验证步骤：网格敏感性测试

无论文献里怎么推荐参数，实际项目里都必须自己验证一次网格敏感性。做法很简单：保持几何、载荷、l 不变，把网格最大尺寸从 h 变成 h/2，再看裂纹路径和注入压力曲线。如果两条曲线基本重合，说明网格已经足够；如果相差明显，说明该加密。

不要忘了同时对比损伤带的宽度。相场法计算的断裂能释放率与 l 有关，损伤带宽度本身是正则化参数，不是绝对的“裂缝开度”。真正的裂缝开度往往要从位移场的不连续趋势去提取，或者配合流体流量来计算。这个细节经常被忽略，却会在写报告和论文时被审稿人直接问住。

5. 参考文献怎么选：理论源头、相场综述与水力压裂专门文献

5.1 奠基性文献：先读懂 Griffith 和相场正则化

做相场法模拟，最有必要精读的第一梯队文献是：

• Bourdin, Francfort, Marigo.Numerical experiments in revisited brittle fracture, JMPS, 2000。 这篇文章把 Griffith 断裂能引入 Mumford-Shah 型泛函的正则化，是相场断裂的源头。读它能理解能量最小化与损伤演化之间的联系。
• Miehe, Welschinger, Hofacker.Thermodynamically consistent phase-field models of fracture, IJNME, 2010。 这篇文章给出了热力学一致的相场断裂框架，也是目前大多数 COMSOL 案例的方程模板来源。尤其是历史场 H 的引入和应变能分裂，方案就是在这一派文献里定型的。
• Miehe, Hofacker, Welschinger.A phase field model for rate-independent crack propagation, CMAME, 2010。

这些文献里的方程基本就是我前文写的演化方程的原型。如果你发现论文里的符号体系跟 COMSOL 不一样，不要慌，按能量泛函的物理项对应着迁移就行。

5.2 综述文献：快速建立全局图景

相场法近几年发展迅速，分支很多，建议看两篇综述：

• Ambati, Gerasimov, De Lorenzis.A review on phase-field models of brittle fracture and a new fast hybrid formulation, European Journal of Mechanics A/Solids, 2015。 这篇文章对各种应变能分裂方案做了对比，区分了各向同性、体积偏量、谱分解等不同做法，对理解“为什么我的裂缝会穿过压缩区”这类问题帮助很大。
• Wu, Nguyen-Thanh.From the classical theories to a phase-field model in brittle fracture, Advances in Engineering Software, 2018。 这篇对统一相场理论框架做了梳理，适合想更进一步做本构改进的人。

综述的价值在于帮你快速定位自己的问题属于哪一类，而不是把时间花在重复推导旧方程上。

5.3 水力压裂专属文献：从岩石断裂到裂缝网络

压裂方向还要单独看一批岩石水力压裂相场文献：

• Wheeler, Wick, Wollner.An augmented-Lagangian method for the phase-field approach for pressurized fractures, Computer Methods in Applied Mechanics and Engineering, 2014。 这篇是水压裂缝压力边界处理的经典，很多后续工作都建立了这个基准。
• Miehe, Mauthe, Teichtmeister.Minimization principles for the coupled problem of Darcy-Biot-type fluid transport, diffusion and fracture in porous media, 2015。 这是一篇把 Biot 孔隙弹性、达西渗流和相场断裂做统一变分原理的文献，压裂耦合的理论根基就在这里。
• Bourdin, Chukwudozie, Yoshioka.A variational approach to the numerical simulation of hydraulic fracturing, SPE Journal, 2012。

读这些文献时，重点看三件事：流体压力是怎么进入裂缝面的、渗透率退化怎么定义的、以及时间尺度上注入速率与裂缝扩展速度的关系。理解了这三个点，再看 COMSOL 案例就不会只是改参数，而是能真正调整模型假设。

5.4 我自己的文献阅读顺序建议

给刚入坑的人一个可以照抄的阅读路线：

1. 先读 Ambati 的综述，理解相场法有哪些主流分支和常见数值坑。
2. 再精读 Miehe 2010 的方程推导，确认历史变量和应变能分裂的来龙去脉。
3. 然后看 Wheeler 2014 的压力边界基准算例，尝试复现它的裂纹路径和压力曲线。
4. 最后根据项目方向选一篇最新的水力压裂相场论文，对比对方的材料参数、长度尺度选择和网格策略。

每换一个新软件、新版本、新物理过程，我都会习惯性回到这些文献里重新核对一遍方程和假设。这个习惯帮我避免了好几次“照着别人模型抄参数”的低级错误——尤其要注意材料断裂能 Gc 的单位，有的文献用 J/m²，有的用 MPa·m，差一个数量级结果就是完全不同的两种裂缝形态。

拿我自己最近跑通的一个案例来说，最深的体会是：相场法不是“换了个软件模块”，而是一整套关于断裂力学的思考方式。如果你只是想要一个能出图的压裂案例，照着我上面第三节的参数填一遍就能跑通；但如果你想让算例经得起推敲，一定要在长度尺度、网格敏感性、历史不可逆这三个地方做足功夫，并且把文献里的方程与 COMSOL 的 PDE 系数一一对应起来。最后再提醒一个容易忽略的操作习惯：每改一次长度尺度 l，就要重新做一遍网格敏感性测试，否则裂纹路径的变化到底是物理结果还是网格结果，你根本分不清楚。这套“原理-建模-调试-文献”的流程走完一遍之后，后面不管是换岩性参数、加水平地应力、还是升级到三维模型，都会顺手很多。