开了这个系列的前一篇之后后台收到不少读者的私信都是围绕一个共同话题“模型确实能跑起来但结果怎么看怎么不对”——要么裂缝形状怪要么压力曲线跟现场数据对不上甚至有人在问“为什么我的二维水力压裂模型里cohesive单元一开裂整个模型就变成一坨浆糊”这篇文章我打算换个角度不重复讲“怎么建模”这种入门内容重点聊清楚几个真正影响二维水力压裂模拟质量的细节内聚区参数怎么标定、注入条件怎么设置、收敛性崩了怎么一步步排查、以及怎么从输出结果里反推模型是否靠谱。这些都是我在实际项目中踩过坑之后沉淀下来的经验希望能帮你把“能跑通”提升到“跑得对”。1. 续篇的定位从“能跑通”到“跑得对”1.1 第一篇之后第二篇该解决什么问题先交代一下背景。上一篇系列文章里我讲了一套基于Cohesive单元的二维水力压裂建模流程从部件创建、材料赋值、分析步设置到最简单的注入模拟都跑通了。那篇主要面向刚上手的人目标是让模型“动起来”。但“动起来”远远不够。搞水力压裂数值模拟的都知道这类问题涉及流固耦合、损伤力学、大变形接触调试难度远高于普通结构分析。哪怕一个看起来很简单的二维模型当cohesive单元的刚度退化与孔压单元里的流体流动相互耦合时数值上任何一个环节匹配不到位结果就会跑偏。这篇文章正是我和我团队在多个项目里反复打磨后沉淀的一套方法论。核心围绕四个问题展开材料参数怎么设才算“合理”、注入边界条件怎么加才能接近真实工况、不收敛的时候怎么排查根因、最终结果怎么验证才能确保可信。同时文中所有内容都以二维平面应变问题为背景与标题中“二维水力压裂”这个关键词严格对齐。1.2 二维模型里的基本假设你未必真正清楚了很多人在建模的时候对二维问题的基本假设不够敏感这为后续一系列异常结果埋了雷。我建议先理清这几个前提第一个假设是平面应变条件。二维水力压裂模型默认裂缝在垂直于平面方向上无限长所有力学响应只在平面内发生变化。这意味着z方向的应变为零但应力不为零。如果你在ABAQUS里选择了平面应力单元那么缝宽计算结果会比平面应变偏大不少接近现场的解释性都很差。第二个假设是对称性。大多数模型沿裂缝轴线取一半并在对称面上施加法向约束。这个做法没有问题但注意cohesive单元的孔压边界条件要分别对待。比如对称面上cohesive单元是没有流体注入的只有裂缝内壁单元才承担注入压力。很多人把pore pressure初始条件全场统一设置导致裂缝还没起裂整个模型已经处于一个水分饱和的非物理状态。第三个假设是各向同性地层。天然页岩或砂岩总有层理面和天然裂缝但二维模型通常忽略这些。忽略的目的是简化问题不代表你可以随意把弹性模量调得很高或很低——至少要与实际岩心的力学实验数据有对应关系。这三个假设如果没在建模前想清楚后面所有参数都会像无根之木。碰巧这也是我判断一个工程师水平的入门问题。2. 内聚区参数整定这是整个模型的命门2.1 要不要显式建立裂缝区域cohesive layer的布置方案做水力压裂模拟时cohesive单元的布置通常分两种一种是在预设裂缝路径上采用zero-thickness cohesive单元另一种是采用有限厚度的cohesive层模拟具有一定宽度的天然裂缝或弱面。绝大多数水力压裂二维模型采用零厚度单元更贴近真实断裂力学行为。原因其实简单水力压裂的裂缝宽度通常只有毫米级一般0.5-3mm如果赋予cohesive单元初始几何厚度那么初始刚度会显著变化从而影响不损伤阶段整体的变形响应。用零厚度单元配合“traction-separation”本构初始刚度完全由材料参数里的Knn、Kss等设定你才能精准控制未开裂时的力学响应。这里还要注意一个细节cohesive单元所占据的几何位置必须在网格划分之前预留。经验做法是在Part模块先把地层切分为两个区域中间留一道缝再在缝的位置划分cohesive单元网格。这样网格拓扑关系清晰后续设置cohesive section时也不会出现面选错的问题。有的朋友喜欢用“tie约束”把cohesive单元两端与实体单元连接到一起这在大多数场合可行。但需要留意tie约束只传递位移和力不直接传递孔压自由度。如果cohesive单元包含孔压自由度比如使用 pore pressure cohesive elementstie约束会带来自由度不匹配的警告。我建议尽量采用共享节点的方式建网格从根源上避免这个问题。2.2 损伤起始判据与压裂液侵入的耦合关系Cohesive单元本构关系里损伤起始判据常见的有最大应力准则、最大应变准则、二次应力准则。对水力压裂问题最常用的是最大应力准则或二次应力准则因为它们与强度理论的物理含义直接挂钩。但这里有一个常常被忽视的点损伤起始并不只由外力决定压裂液侵入会改变裂缝面的有效应力状态。在孔压单元里有效应力 σ σ - αp。当压裂液进入cohesive层孔隙压力p升高有效正应力减小其实加速了拉伸损伤的起始。这意味着单纯在力学边界条件下标定出的损伤参数在孔压耦合情况下会偏保守裂缝会比真实情况更晚起裂。要应对这个问题建议在设置损伤起始应力值时适当参考室内巴西劈裂或三轴抗拉强度试验数据同时考虑现场地应力场的实际分布。如果现场有微地震监测的起裂压力数据用它反推开裂时的临界应力值是更可靠的手段。另外损伤起始判据的类型不要换来换去。有些人在调试时发现裂缝不开裂就顺手把“Quads”改成“Maxe”结果裂缝开了但开得莫名其妙。判据类型的变更相当于改变了整个物理模型调参必须要有方向性每次只动一个变量及时记录。2.3 刚度退化与渗透率演化的联动机制Cohesive单元在abaqus里开裂后刚度从初始值逐渐降到接近零。但如果你用的是带孔压自由度的cohesive单元如COH2D4P那么流体在裂缝内的流动能力与单元损伤状态直接相关。ABAQUS中cohesive单元的流体流动属性包含两个方向法向流动泄漏到地层和切向流动沿裂缝延伸方向。切向流动的渗透率通常根据裂缝开度计算采用立方定律。默认公式我印象中是 k d³ / (12μ)其中d表示裂缝开度μ表示流体黏度。这背后的物理逻辑是平板流假设基于润滑方程推导而来。如果你在cohesive单元的fluid flow面板里直接给一个恒定渗透率就忽略了裂缝开度对导流能力的自增强效应——这在很多工程简化模型里反而更常见。我的经验是二维水力压裂模拟里法向滤失系数和切向立方定律必须同时激活。前者模拟压裂液向地层的滤失后者模拟裂缝内部流动阻力。只激活其中一个都会导致缝内压力分布失真。这种情况下损伤演化参数的设定就变得极其关键。损伤演化常设置为基于断裂能的线性软化或指数软化断裂能Gc的量级一般取几十到几百J/m²取决于岩石类型。断裂能过高裂缝扩展需要更多能量表现为裂缝偏短、缝宽偏大断裂能过低裂缝偏长、缝宽偏小。因此断裂能往往需要通过实际对比KGD解析解或现场净压力数据来标定而不是直接查表取一个默认值。顺便提一个重要参数粘性系数viscosity coefficient。很多人在做动力分析时会设置这个参数来抑制数值振荡默认值0.0001往往过于小在压力剧烈波动的注入阶段会看到cohesive单元的应力场剧烈跳动。我常用1e-4到1e-3量级但这需要根据单元尺寸和时间增量大小权衡——设置过大会延缓损伤演化导致裂缝扩展明显偏慢。3. 注入边界条件的几种做法和它们的坑3.1 恒定流量注入 vs 恒定压力注入二维水力压裂模拟中注入条件无非两类恒定注入流量或恒定注入压力。恒定流量注入更接近矿场施工泵注排量基本是可控的压力是响应量。在ABAQUS里实现方式是给cohesive单元的孔压自由度添加一个pore fluid volume rate边界条件。注意这个边界条件要施加在井筒对应的节点组上一般是裂缝入口的几个节点。恒定压力注入的处理则简单得多直接在入口节点组上固定一个孔压值。但代价是没有了泵注排量的物理约束裂缝扩展行为会跟真实工况出现系统性偏差。比如压力一旦超过起裂压力裂缝会加速扩展因为系统不限制流量供给。基于这些原因我的二维模型几乎全部采用恒定流量注入。但在实际调试中恒定流量注入的麻烦在于初始压力会迅速上升往往一两个增量步内就触发起裂导致起裂瞬间的收敛性极差。如果你遇到这种情况可以先把注入流量按时间做一个斜坡加载比如前0.5秒从0线性增加到目标值。这样做的好处是让孔压场缓慢建立cohesive单元在起裂前经历一个准静态加载过程稳定性明显改善。3.2 滤失系数怎么取值法向滤失leak-off系数是水力压裂模拟里最令新手头疼的参数之一。取值过大压裂液大量漏失缝内压力涨不上去裂缝根本延伸不了取值过小所有液体都用于张开裂缝缝宽偏大压力偏高。真实工况中滤失系数与地层渗透率、压裂液黏度、压差等都有关系。Carter模型给出了一个简化方程把滤失量表达为时间平方根的倒数。ABAQUS中可以直接给cohesive单元的fluid flow属性设置一个leak-off coefficient。我的经验是如果手头没有具体地层测试数据初期模型不妨设置一个偏小的滤失系数比如1e-6 m/s^0.5量级先让模型收敛、跑通裂缝扩展。等整体框架稳定后再通过敏感性分析逐步逼近现场净压力数据。这里的关键是敏感性分析一定要跑几组值把滤失系数和净压力的关系曲线画出来你很容易看到“滤失系数超过某阈值后裂缝无法扩展”的临界行为。3.3 时间步控制策略固定增量、自动增量和流体时间增量不少人在模拟里直接用默认的固定时间增量步然后发现压力曲线在起裂点附近出现锯齿状跳动甚至直接发散。这种情况几乎都是时间增量步过大造成的。水力压裂问题包含多种时间尺度之间的竞争压裂液在cohesive层中扩散需要时间裂缝尖端应力集中更新需要时间损伤演化需要时间。当三者中任意一个过程在一个增量步内发生了“过大”的变化Newton迭代就会失败。我的建议是采用自动增量步并设置一个合理的初始增量时间initial increment比如总模拟时间的1/100到1/1000。ABAQUS的自动增量算法会自动缩步以维持收敛但如果你把初始增量设得太大同样容易直接跳到失败。同时最小增量步要设得足够小比如总时间的1e-6否则起裂瞬间无论如何都过不去。还有一个被很多人忽略的细节孔压单元的稳定时间增量估计。在太小的网格里流体扩散特征时间与力学特征时间分离严重稳定性要求的时间增量可能是力学稳定增量的十分之一甚至百分之一。你可以从warning信息里看到类似“time increment required by fluid flow is smaller than...”的提示这时候就要果断缩小时间增量下限。4. 收敛性问题一次完整的排查链路4.1 一个真实遇到的“零主元”案例一定要记录一次典型的失败案例因为“负特征值”“零主元”这类警告可以说是cohesive模拟中最高频遇到的了。有次我在模拟二维水平裂缝扩展时模型跑到大约第3秒计算就终止了报错信息显示“ZERO PIVOT”。这种报错翻译成人话就是结构出现了机构或者刚度矩阵出现了奇异。而在这个问题里本质原因是cohesive单元在该区域已经发生完全损伤刚度退化到近乎零但裂缝面上没有足够的接触约束来阻止相互穿透。当时我第一时间想到的是接触设定。ABAQUS里cohesive单元完全损伤后如果不设置“cohesive zone contact”或者通用的接触对两边的实体单元就可以自由交叉或分离数值上完全不收敛。加上接触对之后模型勉强能跑但依然存在单元过度畸变的问题原因出在网格太粗裂缝尖端应力梯度无法被良好分辨导致一个增量步内损伤区跨过多个单元。最终解决方案是两路同时推进细化了cohesive层附近的网格长度控制在预期FPZ尺寸的1/3到1/5并给cohesive单元的损伤演化增加了粘性正则化参数。之后模型能够顺畅跑到预定时间。4.2 从warnings列表反推模型设置错误很多新手看到warnings就慌了其实合理的做法是从警告信息里读取“线索”针对性修正。这里列几个水力压裂模型里最常见的警告和它们对应的根因“The plastic/creep/connector friction dissipation energy is high relative to the strain energy”如果消散能占比太高说明单元变形模式异常。常见原因是时间增量步过大或网格太粗。“Strain increment has exceeded fifty times the strain to failure”cohesive单元在本构层面被拉爆了几乎可以确定是断裂能参数过小或者增量步长控制不当。“There are excessive distortion problems”实体单元畸变。配合cohesive层附近的网格细化以及接触设置通常可以缓解。“Time increment is too small”常见于孔压-力学耦合计算此时可以增大viscosity regularization参数或检查边界条件是否被约束得过于刚性。不可否认的是warnings里有一类属于“无害信息”比如接触状态的频繁开闭、材料点的局部失效等。判断的标准是计算最终是否收敛、结果是否出现明显非物理特征。如果都正常不必过度在意。4.3 稳定化的正确打开方式水力压裂分析强烈非线性的特点决定了我们经常需要一些数值稳定化手段。但稳定化不是“哪里不行加哪里”的万金油用错了会让结果失真。在ABAQUS/Standard中自动稳定化automatic stabilization的原理是给单元增加一个阻尼矩阵相当于人为增加能量耗散。对cohesive单元问题稳定化数值设置太小不起作用设置太大则相当于给裂缝扩展加了一堵隐形阻力墙——裂缝扩不动。我的操作习惯是先做“无稳定化”的试跑观察它在哪一步发散、发散前的结果是否合理。如果发散前的结果合理再在收敛失败的分析步中添加自动稳定化并设定一个很小的阻尼因子比如2e-4到5e-4量级同时在job诊断里持续监控stabilization energy占strain energy的比例。只要占比低于1%稳定化对结果的影响就可以忽略。还有一个比较容易忽略的小技巧在Cohesive单元的单元类型选择上优先使用减缩积分单元比如COH2D4而不是COH2D4R这里要注意abaqus的cohesive平面单元只有完整的4节点或6节点版本并没有减缩积分版本。所以这个“技巧”在二维cohesive这里不适用但实体单元请务必用减缩积分避免体积锁定带来的伪刚度。5. 结果解读不要只盯着“裂缝打开了”5.1 从哪里读取裂缝宽度才是对的有相当多的初学者习惯于直接在Visualization模块里看cohesive单元的节点位移取最大张开位移当作裂缝宽度。这个做法的问题很大——节点位移包含刚性位移分量并不等同于裂缝的真实开度。更可靠的方法是看cohesive单元的状态变量特别是SDEGscalar stiffness degradation。当SDEG1时表示单元完全损伤。对于完全损伤的单元我们可以提取其上下表面的相对位移然后计算法向相对位移作为该处裂缝宽度。ABAQUS中可以直接输出COPENcontact opening字段——前提是你设置了接触对。用COPEN来提取裂缝宽度是最直接的。如果没有接触对也可以用场输出中的“relative displacement”或者在cohesive单元节点上做差值运算。但在工程实践中我还是强烈建议哪怕主要分析中没有实际装配接触需求也要在cohesive层位置设立一个自接触方便提取COPEN数据。此外裂缝宽度的提取位置也有讲究。通常应取裂缝尖端后一段稳定扩展区的宽度而不是井筒附近的最大宽度。井筒附近的宽度受注入孔眼应力集中的影响很大反映的是近井效应不是远场扩展特性。5.2 压力曲线怎么看净压力与地层响应的关系恒定流量注入下井底压力或入口孔压随时间变化的曲线是整个模型最重要的“心电图”。典型曲线大致分三个阶段憋压段、起裂段和扩展段。憋压段压力线性上升cohesive单元尚未损伤起裂段压力到达峰值后迅速回落扩展段压力趋于平稳或缓慢下降。这个曲线的形状对地层参数非常敏感。比如断裂韧性偏大时峰值压力偏高滤失系数偏大时扩展段压力明显偏低地应力较高时所有段整体上移。利用这一点你可以快速校准模型参数——通过调整断裂能、滤失系数和地应力让模拟压力曲线与现场微地震或井底压力计数据匹配。但我特别提醒不要把模拟压力曲线与现场压力数据做“逐点重合”式的匹配。二维模型本身忽略了三维裂缝高度约束、层间应力差异等压力幅值存在系统偏差。更合理的做法是对比净压力斜率log-log坐标下的斜率、峰值压力位置、以及扩展段压力变化的趋势。5.3 用KGD解析解验证二维模型的正确性在没有任何现场数据的情况下KGDKristianovich-Geertsma-de Klerk模型是验证二维cohesive单元模型可靠性的最佳标尺。它的关键设定是平面应变条件下、高度固定、裂缝沿水平方向扩展。这恰好与我们构建的二维模型在同一假设框架下。KGD理论给出的裂缝半长、缝口宽度与注入时间的关系可以写成比较简洁的近似表达式。不过我不主张死记公式而是建议把模拟结果和解析解做趋势对比模拟得到的缝长-时间曲线应表现为幂函数关系指数近似落在某个合理范围缝口宽度也随时间增大但增速逐渐放缓。如果模拟结果与解析解指数差异过大比如缝宽反而随时间缩小那基本可以判断模型里某个机制出错了常见原因包括滤失系数过大、断裂能设置过小、地应力设置方向反了。在实际项目中我还经常做另一项验证改变断裂能参数观察缝长-时间曲线的响应方向是否符合预期。如果断裂能增大缝长应该变短如果结果相反说明断裂能对扩展速度已经不具备控制作用此时可能模型已经进入了“滤失主导”或“注入量主导”的状态。这种机理一致性检查比孤零零看一张云图可靠得多。6. 二维模型常见的三个“认知陷阱”6.1 二维模型高度上的单位换算很多人用二维平面应变模型模拟压裂时把模型厚度设成1m或干脆忘记乘上厚度而在给流量边界条件时却直接用了现场施工排量结果缝内压力高了几个数量级起裂时间大幅提前。二维模型的孔压自由度是定义在单位厚度上的。ABAQUS中如果赋予单元厚度1m那么注入流量单位应该对应m³/s/m即每米厚度上的流量。现场施工排量通常是m³/min你要先转换成m³/s然后除以裂缝高度才能得到二维模型里的等效流量。这个看似小学算术但在实际交流中我见到至少三个人因为漏掉高度换算得出荒谬结果后还百思不得其解。6.2 地应力方向对裂缝扩展路径的控制二维水力压裂模型里最大和最小水平主应力通常设为预设路径的法向与切向应力。这本身没有问题但要注意cohesive单元初始应力状态需要与地应力场平衡。很多人只在cohesive单元上设置一个初始化应力场却忘了实体单元也需要相应设置。这会导致初始平衡阶段就出现人为变形——模型还没开始注入就已经产生预位移。正确的操作是对整个模型包括cohesive层和实体域设置初始应力场initial conditions, typestress然后跑一个geostatic分析步让系统达到平衡。如果你用的是ABAQUS需要小心cohesive单元对应力初始化支持程度的限制cohesive单元支持初始应力吗这个问题取决于单元类型。在ABAQUS中带孔压自由度的cohesive单元是可以初始应力场设置的但一定要选择合适的几何方向和输出坐标系。为稳妥起见建议在geostatic步后检查一下初始位移是否在1e-6量级如果位移过大说明地应力初始化有误。6.3 网格尺寸对裂缝宽度的硬约束裂缝宽度是水力压裂模拟最敏感的输出指标之一而它对网格尺寸有天然的硬约束在损伤区FPZ内至少要划分3-5个单元。如果缺口附近网格太粗cohesive单元实际计算出的张开位移会显著低于理论值而且这个问题不会随着网格细化自动消失因为断裂能的分配方式本身是按单元尺寸折算的。也就是说网格越粗数值上等效的断裂能越大裂缝越难扩展。解决这个问题有两条路线。一条是老老实实细化网格把FPZ区域的cohesive单元尺寸控制在毫米级另一条是调整输入断裂能把单元尺寸的数值影响补偿掉即Gc_input Gc_real × L₀/L_elem。这种做法在工程快速评估中常用但要清楚这只是一个修补手段不是真正的收敛解。梯度和层次感上我在实际工作中通常是先跑一个较粗网格的快速模型获得初步裂缝形态和压力响应然后加密关键区域网格作为最终验证。这种“由粗到精”的网格策略比一开始就盲目上细网格效率高很多也更容易在早期排查逻辑错误。写到这里关于二维水力压裂cohesive模型的主要技术点基本上都过了一遍。最后说点我个人在实际操作中的体会很多人觉得这类模拟难在“参数太多”但我认为真正难的是“方程组里的每个参数都不能脱离物理背景单独调”。一次只改一个参数、记录它带来的变化、再决定下一步往哪个方向走——这个笨办法比任何高级优化算法都要可靠。如果你现在正被某个不收敛问题卡住我的建议是先从时间增量上下手把初始增量再缩小一个数量级如果还是不行八成是接触或单元畸变的问题而不是材料参数的问题。排查方向对了问题往往比你想象中简单。