干水力压裂数值模拟这块儿绕不开的话题就是相场法。这两年用 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·dk₀ 是基质渗透率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_f1e-10 m²裂缝渗透率Q2e-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达西定律Darcys 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.3Gc₃/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 代替硬截断。弹性矩阵是否正定。填完横观各向同性参数后在结果里加一个“弹性矩阵特征值”的探针看看有没有负特征值。如果出现负特征值说明参数组合不合理。初始损伤区是否在应力场下产生数值振荡。初始 d1 的区域内刚度已经退化到 κ 的量级如果 κ 取得太小比如 1e-6 以下就会出现大位移振荡。把 κ 设置在 1e-4 到 1e-5 之间通常比较稳。5.2 裂纹路径看起来“糊”或者“碎”网格还是太粗损伤云图如果看起来是一条宽的模糊带而不是清晰的裂纹面第一反应不应该是调求解器而是剖网格。l₀ 取 0.005 m 时必须保证裂纹扩展路径上网格边长不超过 l₀/2。达不到这个条件相场扩散带会在数值上被人为加宽裂纹路径就会失真。网格太粗的另一个表现是裂纹扩展方向出现锯齿状偏移因为相场梯度受网格诱导会沿着网格线走“之”字形。修正办法就是细化路径网格但这个成本确实高。如果计算资源有限建议把模型尺寸缩小到 0.5 m × 0.5 ml₀ 同比例缩小这样网格数量相对可接受。5.3 注入点压力异常检查你的“注入”到底注到哪儿了有些读者反馈压力一直不涨或者涨得特别快。这两种情况背后通常是同一个原因——注入点和模型的连通关系没设对。如果注液点设成了 Dirichlet 压力边界固定压力那压力当然不会涨因为它是被固定住的。正确做法是设成 Neumann 型的注入流量边界。如果用的是达西接口在“流动”边界里选“流入/流出”并给定质量流量或体积流量。反过来如果压力涨得飞快说明流体被“堵住”无法排开裂缝多半是初始裂纹区域还没有完全损伤渗透率太低。给初始损伤区的渗透率一个高地步值比如直接用 k k_f让流体能先在这个区域流动起来压力曲线就会正常。5.4 计算时间太长几个实用降本技巧相场法的计算成本是出了名的门槛尤其是三维模型。如果只是做二维参数扫面可以从这几个方面压计算量利用对称性只建 1/4 模型相场和时间步的容差适当放宽先用粗网格试出裂纹的大致走向再在裂纹路径附近局部细化固定相场求解的次数不要每个时间步都完整求解相场方程有些步里可以让相场保持前一步的值用“分离式求解器”中的“迭代次数限制”选项把每个时间步的相场迭代上限设到 5 次左右实测下来这几招能把整体计算时间压缩到原来的三分之一而裂纹路径几乎不变。6. 给新手的一些心里话我个人在这个项目上踩过的最大一个坑就是一开始迷信“全耦合才是精确的”。实际跑下来发现对于水力压裂这种渗流-应力-损伤三重非线性耦合的问题分离式求解不仅稳而且思路更清晰——先让流体压力找到平衡再让固体和损伤对压力做出响应否则你根本分不清数值振荡到底来自哪个物理过程。另外相场法的核心不是 Comsol 操作而是你对网格尺寸、正则化长度、断裂能这三个参数之间关系的理解。l₀ 取大一点裂纹带就宽断裂能就被“稀释”了取小了网格就要密到爆炸。找到平衡点的唯一办法就是多做几组敏感性测试把 l₀ 和网格尺寸的关系摸透。如果你也是刚把横观各向同性材料引入相场水力压裂模型建议先从各向同性材料开始把注液-起裂-扩展这套流程跑通再逐步加入各向异性弹性、各向异性断裂能。每一步加一个变量出了问题也好定位。磨刀不误砍柴工这个顺序值得遵守。
