★ FEATURED ARTICLE 1. 为什么用相场法模拟压裂传统离散裂缝方法的困境做压裂模拟的人不管是研究页岩气水力压裂还是实验室里的三点弯曲断裂多少都会遇到同一个难题裂缝是一条真正意义上的“不连续面”模拟软件里怎么表达这条不连续面早期主流的做法是离散裂缝模型比如在ABAQUS里插入Cohesive单元在Comsol里用裂纹尖端节点重划分或者干脆用扩展有限元XFEM。这些方法各有各的适用场景但真到复杂压裂工况下问题就来了。首先是Cohesive单元它要求裂缝面沿预设路径走。真实岩石里的天然裂缝、层理面、随机节理可不管你预设没预设裂纹一旦偏转Cohesive单元就失效了。XFEM不需要预设路径但需要显式追踪裂纹尖端位置每个时间步都要判断裂纹尖端在哪里、往哪个方向扩展、扩展多远三维情况下还要考虑裂纹面的形状变化数值实现非常复杂而且多裂纹交汇、分支的时候极容易崩。我用XFEM做过几组算例感觉得出的结论很依赖“裂纹面追踪”的实现细节换个追踪策略结果就变心里很不踏实。相场法Phase Field MethodPFM换了一个完全不同的思路它不把裂缝当成一条“不连续面”而是用一个标量场变量把裂纹“抹”成一个具有一定带宽的连续过渡区。这个思路最早来自脆性断裂的变分理论后来被扩展到水力压裂、多场耦合甚至疲劳断裂、动态断裂。在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−φ)2k其中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换算GcKIC²/EE为平面应变弹性模量E/(1−ν²)。这个换算关系在水力压裂参数标定中非常常用因为很多岩土的KIC数据比Gc数据好查。3. 完整案例的Comsol实现3.1 案例工况设定这个案例模拟的是一个典型的水力压裂过程。尺寸采用实验室常见的平板试件长400mm、宽200mm厚度按平面应变处理。试件中心预置一条长度为40mm的初始裂纹用初始相场φ1的窄带来实现。加载方式有两种选择一是压裂液从裂纹中心注入压力随时间逐渐升高二是右侧施加固定位移载荷模拟劈裂。本文采用前者更贴近水力压裂场景。材料参数取一组典型砂岩数据弹性模量E28GPa泊松比ν0.22抗拉强度3.5MPa断裂韧性KIC1.2MPa·m^(1/2)由此换算出临界能量释放率Gc大约60J/m²。这个数据组和真实砂岩比较接近算出的裂缝形态容易被后续实验验证。这里要提醒一下初始裂纹不要用一个纳米级的几何线切割来表示那样网格剖分很麻烦。简便做法是把初始裂纹区域的材料相场初值直接设为1其余地方设为0。实现方式是在“初始值”设置里给φ赋一个与坐标相关的条件比如如果x0.02且|y|0.001则φ1。这个“初始损伤带”完全可以充当初始裂纹而且带来的收敛麻烦比几何切割少得多。3.2 几何、材料与参数表在Comsol中新建二维模型几何就是一个400×200的矩形初始裂纹区域在模型中心靠左位置。材料节点添加弹性模量与泊松比注意相场断裂接口需要的是“损伤材料的弹性材料”在“脆性断裂”接口内部自带的材料节点里填写即可。我习惯把所有控制参数集中到一个全局参数表里后面做参数化扫描时非常省事参数数值说明E28GPa弹性模量nu0.22泊松比KIC1.2 MPa√m断裂韧性Gc60J/m²临界能量释放率由KIC换算l01.0mm相场长度尺度参数设为网格尺寸的约3倍phi00初始除裂纹带外相场初始值p02MPa注入压力峰值t_pump5s注压时长rho_fluid1000kg/m³压裂液密度mu_fluid1e-3Pa·s压裂液动力黏度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)^2y^2)-0.02)/0.0005))这个表达式的意思是在(0.15,0)点为中心、半径20mm的圆形区域内φ初值接近1裂纹区域外φ接近0完整材料。用tanh做光滑过渡的好处是避免了尖锐跳变引起的初始收敛困难。注意初始裂纹不要直接开在模型正中心略微偏左一点给裂缝扩展留出更大的自由空间这样模拟出来的偏转路径更真实。3.5 网格剖分策略网格是相场法最考验耐心的环节。l01mm的前提下裂纹带内至少要有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是否过大。Gc60J/m²在砂岩里不算离谱但如果你给的是钢材的Gc大几十万怎么加载也不会裂。用一条简单准则判断当裂纹尖端的能量释放率G达到Gc裂缝才会扩展如果模拟压力升到3倍也扩展不了十有八九是Gc数量级不对。最后检查网格。l01mm但网格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和网格尺寸之间的比例认识不够。初学者最容易犯的错误是一上来就猛加密网格算几分钟就崩然后归咎于模型不对。实际经验告诉我先把粗网格跑通再逐级加密观察结果变化趋势反而能更快地找到可靠参数。相场法模拟是一个“精度与成本的游戏”不是说网格越细越好——你要的是一条物理合理的裂纹路径而不是无穷多的网格单元。先把这个平衡掌握住再去追求复杂工况会顺很多。
阅读完成 · 觉得有帮助？
