做压裂模拟的人应该都遇到过这种场景水力裂缝扩展几十步都好好的一到天然裂缝附近要么裂缝停在界面上一动不动要么沿着天然裂缝疯狂拐弯更气人的是模型有时直接不收敛。这个问题落到Abaqus里核心就是水力裂缝与天然裂缝相交的cohesive行为怎么定义。它不是简单地在地层网格里塞一排cohesive单元就能解决参数、网格拓扑、孔压耦合、收敛控制全都纠缠在一起。我做过不少页岩储层水力压裂离散裂缝模拟把自己踩过的坑和能复现的经验整理成一篇完整的实操笔记给正在做类似方向的工程师和研究生参考。1. 相交行为的物理本质水力裂缝在天然裂缝面前“犹豫”的三个原因1.1 不是简单“撞上了就停”很多刚开始用Abaqus做水力裂缝相交问题的人会把天然裂缝想成一道“墙”水力裂缝扩展过去碰到墙然后选择穿过去或者停下来。实际模拟里根本不是这个画面。水力裂缝接近天然裂缝时裂尖前方会产生局部应力扰动天然裂缝作为一个力学不连续面会把这种扰动放大。界面两侧的位移可能不再连续发生滑移或张开天然裂缝内部还会因为剪切滑移产生应力重分布。这个过程非常依赖界面的摩擦特性、断裂韧性和周围地应力状态。水力裂缝与天然裂缝相遇后会形成三种典型模式穿越crossing、捕获并沿天然裂缝张开arrest/dilatation、在天然裂缝另一侧偏移重新起裂offset。这三种模式在Abaqus里不是靠人为指定路径得到的而是由cohesive单元在相交处的损伤、滑移和张开行为自然而然涌现出来的。这是cohesive方法处理相交问题最大的优势也是最难调的地方。1.2 穿越、捕获和偏移的控制因素Abaqus里能不能从模拟结果里看到合理的相交模式前提是你得理解哪些因素在真实物理中起控制作用。综合起来最核心的有三个逼近角、水平应力差、天然裂缝界面强度。逼近角对相交模式的影响很直观。水力裂缝扩展方向基本垂直于最小水平主应力当它遇到一条天然裂缝时如果天然裂缝产状与最大主应力方向的夹角近90°也就是逼近角很大裂缝倾向于直接穿越如果逼近角很小水力裂缝等于斜着撞到界面上界面上剪切分量变大很容易发生沿天然裂缝的滑移和张开。水平应力差是另一个关键变量。当最大与最小水平主应力差值很小时水力裂缝对天然裂缝的“穿透意愿”会明显下降裂缝更愿意沿天然裂缝转向。反过来应力差足够大裂尖的断裂过程区会被“压”向原来的扩展方向天然裂缝界面即使发生剪切滑移也很难让主裂缝转移能量最终实现穿越。最后是天然裂缝界面自身的力学参数。界面摩擦系数越大传递剪切应力的能力越强对穿越有利界面抗拉强度和断裂能越低越容易被流体压力和剪切应力激活水力裂缝就会被“吃掉”变成沿天然裂缝延伸。这就是cohesive行为在相交模拟中必须重点标定天然裂缝参数的原因。1.3 这跟cohesive行为到底怎么对应在Abaqus的cohesive框架里有两类“cohesive行为”容易混淆。第一类是岩石基体新断裂面的cohesive行为用来描述水力裂缝裂尖前方断裂过程区中的损伤演化第二类是天然裂缝自身界面的cohesive行为用来描述既有弱面的张开与剪切滑移。相交点附近的数值过程本质上是这两类cohesive行为在竞争。水力裂缝裂尖的应力集中达到一定水平后岩石基体会起裂但如果天然裂缝界面先发生损伤或者界面在剪切作用下发生显著滑移裂尖应力集中就会被释放掉水力裂缝没法继续往前走。所以在建模时我通常会建立两套cohesive层一套沿着水力裂缝预设扩展路径另一套沿着天然裂缝走向。两套cohesive的牵引-分离参数完全不同尤其是断裂能和黏聚强度。这个做法的直接好处是相交模式不是人工指定的而是通过两条Cohesive路径上的SDEG和CSDMG演化自动判断出来的。2. Cohesive材料卡片参数怎么定才对不能直接套默认值2.1 选什么单元类型关键看有没有孔压耦合普通力学分析里用COH2D4、COH3D8这类cohesive单元就行但水力压裂必须考虑流体在裂缝里的流动和滤失。没有孔压自由度的cohesive单元只能算“力学上断开”流体压力没法在裂缝里传递算出来的缝宽和裂缝路径可能看着对但本质上不是水力裂缝而是力学裂缝。我在二维水力压裂模型里常用COH2D4P三维模型里用COH3D8P。这类带孔压自由度的cohesive单元节点除了位移自由度之外还有孔隙压力自由度Abaqus才能把流体流动和骨架变形耦合起来。如果你的问题不关心缝内流体压力和滤失只要一个远场应力作用下的裂缝扩展路径那用普通cohesive单元做力学分析也能应付但只要涉及注液、憋压、裂缝内流体前锋推进就必须上带孔压的cohesive单元。2.2 牵引-分离曲线里的三个刚度参数Abaqus的cohesive行为默认是牵引-分离关系也就是界面上的牵引力与相对位移之间的关系。这个关系里三个最基础的值初始刚度、损伤起始应力、断裂能。初始刚度处理起来很容易出问题。它不是一个真实的岩石弹性模量而是界面在损伤之前抵抗相对位移的刚度。选得太小整个模型会变得额外“软”裂缝还没起裂远场应力加载就产生了不正常的变形选得太大又容易引起矩阵病态导致收敛困难。我的经验是先按最终位移的5到10倍估算把损伤起始位移控制在最终破坏位移的1/10到1/5量级。举个例子如果黏聚强度取5MPa断裂能取70N/m按线性软化估算最终位移是28微米那初始刚度可以取1.7e12 Pa/m量级也就是让前端弹性位移不到3微米。损伤起始应力一般由岩石的抗拉强度或者天然裂缝的抗剪强度决定。基岩的cohesive层可以取5MPa量级天然裂缝的cohesive层要根据露头、岩心或者测井解释来通常只有0.1到0.5MPa甚至更低。断裂能则由断裂韧性换算得到平面应变下可以用Gc K_IC^2(1-ν²)/E来估算。2.3 从断裂韧性到cohesive参数的一个换算例子我用一个具体数据串一遍。假设岩石弹性模量E30GPa泊松比ν0.25I型断裂韧性K_IC1.5MPa√m。平面应变条件下断裂能就是G_IC (1-ν²) × K_IC² / E代入0.9375 × 2.25e12 / 30e9约等于70.3N/m。如果黏聚强度T0取5MPa线性软化下的最终分离位移δf2G_IC/T0≈28微米。这个值很小说明CZM区域很薄。为了避免数值问题初始刚度通常会让损伤起始位移取δf的1/10左右也就是2.8微米K0T0/δ0≈1.78e12Pa/m。II型断裂能G_IIc一般比G_IC高岩石类材料可以取180到300N/m混合模式用BK准则幂指数取1.5左右比较常见。Abaqus材料卡示意如下*MATERIAL, NAMEHF_CZM *ELASTIC, TYPETRACTION 1.78e12, 1.78e12, 7.12e11 *DAMAGE INITIATION, CRITERIONMAXS 5.0e6, 5.0e6, 1.0e6 *DAMAGE EVOLUTION, TYPEENERGY, MIXED MODE BEHAVIORBK, POWER1.5 70.3, 180.0, 300.0 *DAMAGE STABILIZATION 1.0e-5这只是材料骨架不是完整inp。实际提交前还需要根据单位制检查绝不能把MPa和Pa混在一个模型里。2.4 天然裂缝和基岩cohesive层必须区别对待在相交问题里最怕的就是图省事给天然裂缝和基岩用同一套cohesive参数。这样做的结果几乎永远是错误相交模式水力裂缝在天然裂缝面前不是穿得太容易就是停得太干脆。我的习惯是把天然裂缝界面单独定义一套材料核心参数和基岩拉开两个量级。天然裂缝的黏聚强度可以取基岩的1/10甚至更低断裂能取20到40N/m摩擦系数取0.4到0.8然后通过设置摩擦属性来控制剪切滑移量。天然裂缝初始刚度和基岩cohesive层类似因为如果太低地应力加载阶段天然裂缝就会提前“凹”进去裂缝还没开泵整个模型已经产生了大量虚假位移。另外还要注意天然裂缝界面通常不是一个没有初始强度的理想面而是一个有少量胶结物或者矿物填充的弱面。所以参数取法一般介于“完全光滑”和“基岩本身”之间。具体取值要靠现场资料和实验数据但模拟规律一定要清楚黏聚强度低、断裂能低、摩擦系数低天然裂缝就容易被激活反过来天然裂缝强度接近基岩水力裂缝就倾向于穿越。3. 相交点附近的网格与拓扑处理最容易翻车也最影响收敛3.1 零厚度cohesive和有限厚度cohesive怎么选Abaqus的cohesive单元可以建造成零厚度也可以建造成有限厚度。水力压裂模拟中我几乎都用零厚度cohesive层来代表裂缝面。天然裂缝是地质弱面本身不应该有真实的“厚度”参与力学计算如果给一个有限厚度比如0.1mm虽然网格好做一点但会在界面两侧引入额外的抗弯刚度裂缝开度结果反而失真。零厚度cohesive做起来繁琐一些因为需要在裂缝路径两侧保留重合节点。Abaqus里可以在网格生成后用脚本或手工把裂缝路径上的单元重新编号让cohesive层与两侧实体单元共节点而不是靠绑定Tie。Tie虽然省事但会在交界面引入额外的约束柔度对裂缝形态非常敏感。对于带孔压的cohesive单元建议一开始就规划好单元方向保证cohesive单元厚度方向与裂缝面法向一致。弄反了裂缝开度方向会错乱后处理时看不出来但张开的几何完全不对。3.2 网格尺寸与cohesive过程区长度匹配水力裂缝的cohesive过程区很薄单元尺寸稍大数值结果就不是“裂缝扩展”而是“整排单元同时破坏”。一个快速估算方法是用E×Gc/T0²得到过程区长度量级。拿前面的例子算30GPa×70N/m÷5MPa²量级大概只有几十毫米。但实际模拟中裂缝可能长几十米甚至上百米整体网格都加密到毫米级不现实。我一般只在裂尖可能经过的路径上局部加密网格尺寸控制在过程区长度的一半左右远离裂缝路径的地方逐渐放大。天然裂缝与水力裂缝相交的区域在中点附近还要额外加密最好让交汇处两三个网格尺寸基本一致。这样能显著减少裂缝在相交点产生非自然偏转的数值伪影。3.3 交叉点拓扑两条cohesive层怎么接二维模型里水力裂缝路径和天然裂缝路径是一条线交叉交叉点本身只是一个节点。两条cohesive单元链在这个节点上共节点就可以问题不大。但三维模型里水力裂缝面和天然裂缝面在空间中交于一条线两条零厚度层在相交线上会发生严重的单元重叠和过约束导致负特征值、零主元这些典型报错。处理三维交叉线我的经验是把天然裂缝在交叉线附近做一个小半径的“截断”或者在水力裂缝路径上预先预留一个通道避免两套cohesive层在同一片区域里同时存在。这两套界面可以共享交叉线上的节点但不要让它们在交叉线附近有厚度重叠。否则不是不收敛就是算出来天然裂缝在交汇处出现不正常的刚性突起。3.4 初始闭合状态和初始应力场cohesive层建模完成后第一件事是地应力平衡。不要一上来就注液先把远场地应力和初始孔隙压力加载好让cohesive单元处于一个闭合但未损伤的状态。如果cohesive初始刚度或者材料定义有问题地应力平衡阶段SDEG就会悄悄变成0.1甚至更高后面水力裂缝还没到天然裂缝就已经损伤了。还要检查cohesive层有没有初始穿透。零厚度cohesive经常因为网格划分顺序问题在局部交叉点出现极小的几何重叠。这个问题在Abaqus里不容易从云图看出来但会不断产生严重不连续迭代。最直接的办法是在加载之前专门输出一次SDEG和接触状态如果发现初始损伤已经存在就回头检查网格拓扑。4. 为什么inp越跑越慢并行、阻尼和增量步的取舍4.1 运行慢的第一来源不是求解器而是单元刚度很多人遇到“Abaqus inp文件运行很慢”第一反应是换电脑或者加核数其实水力压裂cohesive模型卡顿的根源经常在单元本身。cohesive单元初始刚度设得过高导致整个模型的刚度矩阵条件数变差隐式求解每一步都要反复迭代。即使能收敛每步的计算量也比普通实体单元大得多。另一个常见问题是cohesive单元尺寸太小。裂缝路径附近局部加密后如果其余区域网格没做合理过渡最小网格尺寸会被实验数据拉得很低时间增量步长被最小网格卡死整个模拟几百步都推不动。所以速度优化不能只盯并行核数应先检查最小网格尺度、cohesive初始刚度和加载幅值。4.2 黏性正则化是收敛三件套的主角cohesive单元发生损伤后刚度下降进入软化段问题就变成了“负刚度”。隐式求解器在这种阶段非常容易发散。Abaqus提供了一个专门工具黏性正则化也就是在材料卡里加*DAMAGE STABILIZATION。它给损伤演化加一个黏性阻力让软化段不瞬间降到零避免严重的负特征值。黏性正则化系数不能拍脑袋。给太大裂缝扩展会被明显延迟甚至出现“伪延性”裂缝张开量失真给太小又起不到稳定作用。我一般的流程是先给1e-5算不过去就逐步加到1e-4最多到1e-3。如果系数已经到1e-3还在发散那问题大概率不是阻尼不够而是网格或拓扑出了问题。加载方式也要配合。在注液步我不建议直接给一个阶跃流速最好用Smooth Step幅值曲线把注液速率从0平滑加到目标值。这样能避免初始冲击引起的瞬态剧烈损伤尤其能让天然裂缝相交点附近避免出现一步之内SDEG突然从0变1的极端状态。4.3 并行核数、mpi threads和Linux提交细节并行设置是提速的最后一步不是第一步。Abaqus/Standard和Abaqus/Explicit的并行机制不一样。Standard主要靠稀疏矩阵求解器并行Explicit靠区域分解。很多人设置cpus16以为一定比cpus8快实际上当模型单元数不够大时MPI通信开销早就把计算收益吃掉了。我常用的提交命令类似abaqus jobfrac_model cpus8 doubleboth output_precisionsingle interdoubleboth会让所有计算结果用双精度对cohesive这种小位移、大刚度问题很有帮助但速度会变慢输出精度反而不需要双精度所以output_precision保持single即可。在Linux下我习惯用nohup后台提交避免ssh断开导致任务中断nohup abaqus jobfrac_model cpus8 doubleboth run.log 21 关于mpi threads我个人的体会是单机多核不需要刻意调MPI线程数直接给cpus就行。跨节点并行才需要仔细配置MPI环境和主机列表。如果你发现cpus从8加到16几乎没加速甚至更慢多半是核数超过内存带宽或物理核心数。4.4 常见报错和对应的处理思路下面这几个报错是cohesive水力压裂模型里出现频率最高的报错现象常见原因处理手段Negative eigenvaluecohesive单元软化段进入负刚度调整*DAMAGE STABILIZATION系数Too many attempts时间增量被压到下限检查网格最小尺寸、加载幅值、是否漏掉黏性正则化Zero pivot相交点过约束或tie重复重新检查交叉点拓扑去掉重复约束Pore pressure change too large注液速率过快压力增量步太大加密时间步采用平滑注液曲线或适当减小最小增量步The solution appears to be diverging迭代发散打开非对称求解器切回Single precision检查精度另外带孔压的cohesive模型在Standard中建议开启非对称求解器。孔压-应力耦合矩阵天然非对称用默认对称求解器会牺牲收敛性和速度。这个选项在Step设置里可以直接打开虽然单步内存占用会增加但总时长往往能降下来。5. 后处理怎么看用SDEG、孔隙压力和缝宽判断相交模式5.1 先分清SDEG和CSDMGAbaqus后处理里cohesive单元有两个经常被搞混的输出量SDEG和CSDMG。SDEG是标量刚度退化系数反映单元从完整到完全失效的整体退化程度0表示未损伤1表示完全断裂。CSDMG是损伤起始后的标量损伤变量重点看损伤什么时候、在哪个位置开始出现。判断水力裂缝是否穿过天然裂缝我习惯同时看这两个量。SDEG能告诉你裂缝主路径在哪里贯通CSDMG能告诉你天然裂缝什么时候被“激活”。只看纤维一样细的云图不够还得结合孔隙压力分布和位移场确认流体是不是真的从裂缝里流过去了。5.2 三种相交结果的典型后处理特征如果相交模式是穿越后处理画面通常是这样水力裂缝路径上的SDEG从注入点延伸到天然裂缝另一侧中间没有中断天然裂缝SDEG在交点附近可能有一小段非零区域但数值没有到1或者即使到1也没有形成连续的张开通道。孔隙压力分布上压力前锋能连续越过天然裂缝。如果相交模式是捕获并沿天然裂缝张开你会看到水力裂缝的SDEG在天然裂缝面前停下来主裂缝尖端不再向前扩展天然裂缝上的SDEG从交点向两侧扩展并且CSDMG演化区域明显比交点处宽。这时候孔隙压力云图会出现一个明显变化压力从主裂缝尖端“拐弯”进入了天然裂缝流体前锋沿天然裂缝继续推进。如果相交模式是偏移特征最微妙天然裂缝界面局部有损伤但并没有形成像捕获模式那样的大范围张开水力裂缝在天然裂缝另一侧某个偏移距离重新起裂整体路径出现一个台阶。这种模式最容易和穿越误判需要把SDEG、缝宽和孔隙压力放在同一张图里对比。5.3 缝宽和压力前沿怎么输出缝宽是判断相交模式的重要辅助量但Abaqus不会直接给一个叫“缝宽”的输出。cohesive单元完全破坏后缝宽近似等于上壁面和下壁面节点之间的法向位移差。在后处理里我一般把cohesive上下表面节点位移提取出来做一个简单的坐标差就是实际裂缝张开宽度。孔隙压力前沿的判断则要小心。裂缝内的孔隙压力应该比远场初始孔隙压力高而且沿着裂缝有压力梯度如果在cohesive单元上出现压力跳变或者局部负压先检查是不是单元方向、初始孔隙压力或者单位制出了问题。根本不是参数问题而往往是孔压初始条件没设对。5.4 结果不合理时按什么顺序排查如果模拟结果和现场微地震、井口压力曲线完全对不上我通常按固定顺序排查而不是乱调参数先看单位制这个Bug最容易藏在角落里再看初始应力平衡cohesive是否提前损伤第三看网格尺寸是否超过cohesive过程区长度第四看相交点拓扑第五看黏性正则化系数是否过大最后才调天然裂缝的黏聚强度和断裂能。这个顺序的好处是前面任何一个没检查好调后面的参数都没用反而越调越玄。我个人在实际项目里最深的体会是Abaqus里cohesive相交模拟的问题90%都出在参数和网格之前的“关系”上而不是求解器本身。把天然裂缝界面看成有摩擦、有黏聚强度、有断裂能的真实弱面把水力裂缝裂尖看成有损伤区的受力对象模型才能给出可信的相交模式。后面再做随机多裂缝扩展、离散裂缝网络压裂这类复杂问题也都是在这个基础上叠加和扩展。
Abaqus水力裂缝与天然裂缝相交的Cohesive行为模拟实战
做压裂模拟的人应该都遇到过这种场景水力裂缝扩展几十步都好好的一到天然裂缝附近要么裂缝停在界面上一动不动要么沿着天然裂缝疯狂拐弯更气人的是模型有时直接不收敛。这个问题落到Abaqus里核心就是水力裂缝与天然裂缝相交的cohesive行为怎么定义。它不是简单地在地层网格里塞一排cohesive单元就能解决参数、网格拓扑、孔压耦合、收敛控制全都纠缠在一起。我做过不少页岩储层水力压裂离散裂缝模拟把自己踩过的坑和能复现的经验整理成一篇完整的实操笔记给正在做类似方向的工程师和研究生参考。1. 相交行为的物理本质水力裂缝在天然裂缝面前“犹豫”的三个原因1.1 不是简单“撞上了就停”很多刚开始用Abaqus做水力裂缝相交问题的人会把天然裂缝想成一道“墙”水力裂缝扩展过去碰到墙然后选择穿过去或者停下来。实际模拟里根本不是这个画面。水力裂缝接近天然裂缝时裂尖前方会产生局部应力扰动天然裂缝作为一个力学不连续面会把这种扰动放大。界面两侧的位移可能不再连续发生滑移或张开天然裂缝内部还会因为剪切滑移产生应力重分布。这个过程非常依赖界面的摩擦特性、断裂韧性和周围地应力状态。水力裂缝与天然裂缝相遇后会形成三种典型模式穿越crossing、捕获并沿天然裂缝张开arrest/dilatation、在天然裂缝另一侧偏移重新起裂offset。这三种模式在Abaqus里不是靠人为指定路径得到的而是由cohesive单元在相交处的损伤、滑移和张开行为自然而然涌现出来的。这是cohesive方法处理相交问题最大的优势也是最难调的地方。1.2 穿越、捕获和偏移的控制因素Abaqus里能不能从模拟结果里看到合理的相交模式前提是你得理解哪些因素在真实物理中起控制作用。综合起来最核心的有三个逼近角、水平应力差、天然裂缝界面强度。逼近角对相交模式的影响很直观。水力裂缝扩展方向基本垂直于最小水平主应力当它遇到一条天然裂缝时如果天然裂缝产状与最大主应力方向的夹角近90°也就是逼近角很大裂缝倾向于直接穿越如果逼近角很小水力裂缝等于斜着撞到界面上界面上剪切分量变大很容易发生沿天然裂缝的滑移和张开。水平应力差是另一个关键变量。当最大与最小水平主应力差值很小时水力裂缝对天然裂缝的“穿透意愿”会明显下降裂缝更愿意沿天然裂缝转向。反过来应力差足够大裂尖的断裂过程区会被“压”向原来的扩展方向天然裂缝界面即使发生剪切滑移也很难让主裂缝转移能量最终实现穿越。最后是天然裂缝界面自身的力学参数。界面摩擦系数越大传递剪切应力的能力越强对穿越有利界面抗拉强度和断裂能越低越容易被流体压力和剪切应力激活水力裂缝就会被“吃掉”变成沿天然裂缝延伸。这就是cohesive行为在相交模拟中必须重点标定天然裂缝参数的原因。1.3 这跟cohesive行为到底怎么对应在Abaqus的cohesive框架里有两类“cohesive行为”容易混淆。第一类是岩石基体新断裂面的cohesive行为用来描述水力裂缝裂尖前方断裂过程区中的损伤演化第二类是天然裂缝自身界面的cohesive行为用来描述既有弱面的张开与剪切滑移。相交点附近的数值过程本质上是这两类cohesive行为在竞争。水力裂缝裂尖的应力集中达到一定水平后岩石基体会起裂但如果天然裂缝界面先发生损伤或者界面在剪切作用下发生显著滑移裂尖应力集中就会被释放掉水力裂缝没法继续往前走。所以在建模时我通常会建立两套cohesive层一套沿着水力裂缝预设扩展路径另一套沿着天然裂缝走向。两套cohesive的牵引-分离参数完全不同尤其是断裂能和黏聚强度。这个做法的直接好处是相交模式不是人工指定的而是通过两条Cohesive路径上的SDEG和CSDMG演化自动判断出来的。2. Cohesive材料卡片参数怎么定才对不能直接套默认值2.1 选什么单元类型关键看有没有孔压耦合普通力学分析里用COH2D4、COH3D8这类cohesive单元就行但水力压裂必须考虑流体在裂缝里的流动和滤失。没有孔压自由度的cohesive单元只能算“力学上断开”流体压力没法在裂缝里传递算出来的缝宽和裂缝路径可能看着对但本质上不是水力裂缝而是力学裂缝。我在二维水力压裂模型里常用COH2D4P三维模型里用COH3D8P。这类带孔压自由度的cohesive单元节点除了位移自由度之外还有孔隙压力自由度Abaqus才能把流体流动和骨架变形耦合起来。如果你的问题不关心缝内流体压力和滤失只要一个远场应力作用下的裂缝扩展路径那用普通cohesive单元做力学分析也能应付但只要涉及注液、憋压、裂缝内流体前锋推进就必须上带孔压的cohesive单元。2.2 牵引-分离曲线里的三个刚度参数Abaqus的cohesive行为默认是牵引-分离关系也就是界面上的牵引力与相对位移之间的关系。这个关系里三个最基础的值初始刚度、损伤起始应力、断裂能。初始刚度处理起来很容易出问题。它不是一个真实的岩石弹性模量而是界面在损伤之前抵抗相对位移的刚度。选得太小整个模型会变得额外“软”裂缝还没起裂远场应力加载就产生了不正常的变形选得太大又容易引起矩阵病态导致收敛困难。我的经验是先按最终位移的5到10倍估算把损伤起始位移控制在最终破坏位移的1/10到1/5量级。举个例子如果黏聚强度取5MPa断裂能取70N/m按线性软化估算最终位移是28微米那初始刚度可以取1.7e12 Pa/m量级也就是让前端弹性位移不到3微米。损伤起始应力一般由岩石的抗拉强度或者天然裂缝的抗剪强度决定。基岩的cohesive层可以取5MPa量级天然裂缝的cohesive层要根据露头、岩心或者测井解释来通常只有0.1到0.5MPa甚至更低。断裂能则由断裂韧性换算得到平面应变下可以用Gc K_IC^2(1-ν²)/E来估算。2.3 从断裂韧性到cohesive参数的一个换算例子我用一个具体数据串一遍。假设岩石弹性模量E30GPa泊松比ν0.25I型断裂韧性K_IC1.5MPa√m。平面应变条件下断裂能就是G_IC (1-ν²) × K_IC² / E代入0.9375 × 2.25e12 / 30e9约等于70.3N/m。如果黏聚强度T0取5MPa线性软化下的最终分离位移δf2G_IC/T0≈28微米。这个值很小说明CZM区域很薄。为了避免数值问题初始刚度通常会让损伤起始位移取δf的1/10左右也就是2.8微米K0T0/δ0≈1.78e12Pa/m。II型断裂能G_IIc一般比G_IC高岩石类材料可以取180到300N/m混合模式用BK准则幂指数取1.5左右比较常见。Abaqus材料卡示意如下*MATERIAL, NAMEHF_CZM *ELASTIC, TYPETRACTION 1.78e12, 1.78e12, 7.12e11 *DAMAGE INITIATION, CRITERIONMAXS 5.0e6, 5.0e6, 1.0e6 *DAMAGE EVOLUTION, TYPEENERGY, MIXED MODE BEHAVIORBK, POWER1.5 70.3, 180.0, 300.0 *DAMAGE STABILIZATION 1.0e-5这只是材料骨架不是完整inp。实际提交前还需要根据单位制检查绝不能把MPa和Pa混在一个模型里。2.4 天然裂缝和基岩cohesive层必须区别对待在相交问题里最怕的就是图省事给天然裂缝和基岩用同一套cohesive参数。这样做的结果几乎永远是错误相交模式水力裂缝在天然裂缝面前不是穿得太容易就是停得太干脆。我的习惯是把天然裂缝界面单独定义一套材料核心参数和基岩拉开两个量级。天然裂缝的黏聚强度可以取基岩的1/10甚至更低断裂能取20到40N/m摩擦系数取0.4到0.8然后通过设置摩擦属性来控制剪切滑移量。天然裂缝初始刚度和基岩cohesive层类似因为如果太低地应力加载阶段天然裂缝就会提前“凹”进去裂缝还没开泵整个模型已经产生了大量虚假位移。另外还要注意天然裂缝界面通常不是一个没有初始强度的理想面而是一个有少量胶结物或者矿物填充的弱面。所以参数取法一般介于“完全光滑”和“基岩本身”之间。具体取值要靠现场资料和实验数据但模拟规律一定要清楚黏聚强度低、断裂能低、摩擦系数低天然裂缝就容易被激活反过来天然裂缝强度接近基岩水力裂缝就倾向于穿越。3. 相交点附近的网格与拓扑处理最容易翻车也最影响收敛3.1 零厚度cohesive和有限厚度cohesive怎么选Abaqus的cohesive单元可以建造成零厚度也可以建造成有限厚度。水力压裂模拟中我几乎都用零厚度cohesive层来代表裂缝面。天然裂缝是地质弱面本身不应该有真实的“厚度”参与力学计算如果给一个有限厚度比如0.1mm虽然网格好做一点但会在界面两侧引入额外的抗弯刚度裂缝开度结果反而失真。零厚度cohesive做起来繁琐一些因为需要在裂缝路径两侧保留重合节点。Abaqus里可以在网格生成后用脚本或手工把裂缝路径上的单元重新编号让cohesive层与两侧实体单元共节点而不是靠绑定Tie。Tie虽然省事但会在交界面引入额外的约束柔度对裂缝形态非常敏感。对于带孔压的cohesive单元建议一开始就规划好单元方向保证cohesive单元厚度方向与裂缝面法向一致。弄反了裂缝开度方向会错乱后处理时看不出来但张开的几何完全不对。3.2 网格尺寸与cohesive过程区长度匹配水力裂缝的cohesive过程区很薄单元尺寸稍大数值结果就不是“裂缝扩展”而是“整排单元同时破坏”。一个快速估算方法是用E×Gc/T0²得到过程区长度量级。拿前面的例子算30GPa×70N/m÷5MPa²量级大概只有几十毫米。但实际模拟中裂缝可能长几十米甚至上百米整体网格都加密到毫米级不现实。我一般只在裂尖可能经过的路径上局部加密网格尺寸控制在过程区长度的一半左右远离裂缝路径的地方逐渐放大。天然裂缝与水力裂缝相交的区域在中点附近还要额外加密最好让交汇处两三个网格尺寸基本一致。这样能显著减少裂缝在相交点产生非自然偏转的数值伪影。3.3 交叉点拓扑两条cohesive层怎么接二维模型里水力裂缝路径和天然裂缝路径是一条线交叉交叉点本身只是一个节点。两条cohesive单元链在这个节点上共节点就可以问题不大。但三维模型里水力裂缝面和天然裂缝面在空间中交于一条线两条零厚度层在相交线上会发生严重的单元重叠和过约束导致负特征值、零主元这些典型报错。处理三维交叉线我的经验是把天然裂缝在交叉线附近做一个小半径的“截断”或者在水力裂缝路径上预先预留一个通道避免两套cohesive层在同一片区域里同时存在。这两套界面可以共享交叉线上的节点但不要让它们在交叉线附近有厚度重叠。否则不是不收敛就是算出来天然裂缝在交汇处出现不正常的刚性突起。3.4 初始闭合状态和初始应力场cohesive层建模完成后第一件事是地应力平衡。不要一上来就注液先把远场地应力和初始孔隙压力加载好让cohesive单元处于一个闭合但未损伤的状态。如果cohesive初始刚度或者材料定义有问题地应力平衡阶段SDEG就会悄悄变成0.1甚至更高后面水力裂缝还没到天然裂缝就已经损伤了。还要检查cohesive层有没有初始穿透。零厚度cohesive经常因为网格划分顺序问题在局部交叉点出现极小的几何重叠。这个问题在Abaqus里不容易从云图看出来但会不断产生严重不连续迭代。最直接的办法是在加载之前专门输出一次SDEG和接触状态如果发现初始损伤已经存在就回头检查网格拓扑。4. 为什么inp越跑越慢并行、阻尼和增量步的取舍4.1 运行慢的第一来源不是求解器而是单元刚度很多人遇到“Abaqus inp文件运行很慢”第一反应是换电脑或者加核数其实水力压裂cohesive模型卡顿的根源经常在单元本身。cohesive单元初始刚度设得过高导致整个模型的刚度矩阵条件数变差隐式求解每一步都要反复迭代。即使能收敛每步的计算量也比普通实体单元大得多。另一个常见问题是cohesive单元尺寸太小。裂缝路径附近局部加密后如果其余区域网格没做合理过渡最小网格尺寸会被实验数据拉得很低时间增量步长被最小网格卡死整个模拟几百步都推不动。所以速度优化不能只盯并行核数应先检查最小网格尺度、cohesive初始刚度和加载幅值。4.2 黏性正则化是收敛三件套的主角cohesive单元发生损伤后刚度下降进入软化段问题就变成了“负刚度”。隐式求解器在这种阶段非常容易发散。Abaqus提供了一个专门工具黏性正则化也就是在材料卡里加*DAMAGE STABILIZATION。它给损伤演化加一个黏性阻力让软化段不瞬间降到零避免严重的负特征值。黏性正则化系数不能拍脑袋。给太大裂缝扩展会被明显延迟甚至出现“伪延性”裂缝张开量失真给太小又起不到稳定作用。我一般的流程是先给1e-5算不过去就逐步加到1e-4最多到1e-3。如果系数已经到1e-3还在发散那问题大概率不是阻尼不够而是网格或拓扑出了问题。加载方式也要配合。在注液步我不建议直接给一个阶跃流速最好用Smooth Step幅值曲线把注液速率从0平滑加到目标值。这样能避免初始冲击引起的瞬态剧烈损伤尤其能让天然裂缝相交点附近避免出现一步之内SDEG突然从0变1的极端状态。4.3 并行核数、mpi threads和Linux提交细节并行设置是提速的最后一步不是第一步。Abaqus/Standard和Abaqus/Explicit的并行机制不一样。Standard主要靠稀疏矩阵求解器并行Explicit靠区域分解。很多人设置cpus16以为一定比cpus8快实际上当模型单元数不够大时MPI通信开销早就把计算收益吃掉了。我常用的提交命令类似abaqus jobfrac_model cpus8 doubleboth output_precisionsingle interdoubleboth会让所有计算结果用双精度对cohesive这种小位移、大刚度问题很有帮助但速度会变慢输出精度反而不需要双精度所以output_precision保持single即可。在Linux下我习惯用nohup后台提交避免ssh断开导致任务中断nohup abaqus jobfrac_model cpus8 doubleboth run.log 21 关于mpi threads我个人的体会是单机多核不需要刻意调MPI线程数直接给cpus就行。跨节点并行才需要仔细配置MPI环境和主机列表。如果你发现cpus从8加到16几乎没加速甚至更慢多半是核数超过内存带宽或物理核心数。4.4 常见报错和对应的处理思路下面这几个报错是cohesive水力压裂模型里出现频率最高的报错现象常见原因处理手段Negative eigenvaluecohesive单元软化段进入负刚度调整*DAMAGE STABILIZATION系数Too many attempts时间增量被压到下限检查网格最小尺寸、加载幅值、是否漏掉黏性正则化Zero pivot相交点过约束或tie重复重新检查交叉点拓扑去掉重复约束Pore pressure change too large注液速率过快压力增量步太大加密时间步采用平滑注液曲线或适当减小最小增量步The solution appears to be diverging迭代发散打开非对称求解器切回Single precision检查精度另外带孔压的cohesive模型在Standard中建议开启非对称求解器。孔压-应力耦合矩阵天然非对称用默认对称求解器会牺牲收敛性和速度。这个选项在Step设置里可以直接打开虽然单步内存占用会增加但总时长往往能降下来。5. 后处理怎么看用SDEG、孔隙压力和缝宽判断相交模式5.1 先分清SDEG和CSDMGAbaqus后处理里cohesive单元有两个经常被搞混的输出量SDEG和CSDMG。SDEG是标量刚度退化系数反映单元从完整到完全失效的整体退化程度0表示未损伤1表示完全断裂。CSDMG是损伤起始后的标量损伤变量重点看损伤什么时候、在哪个位置开始出现。判断水力裂缝是否穿过天然裂缝我习惯同时看这两个量。SDEG能告诉你裂缝主路径在哪里贯通CSDMG能告诉你天然裂缝什么时候被“激活”。只看纤维一样细的云图不够还得结合孔隙压力分布和位移场确认流体是不是真的从裂缝里流过去了。5.2 三种相交结果的典型后处理特征如果相交模式是穿越后处理画面通常是这样水力裂缝路径上的SDEG从注入点延伸到天然裂缝另一侧中间没有中断天然裂缝SDEG在交点附近可能有一小段非零区域但数值没有到1或者即使到1也没有形成连续的张开通道。孔隙压力分布上压力前锋能连续越过天然裂缝。如果相交模式是捕获并沿天然裂缝张开你会看到水力裂缝的SDEG在天然裂缝面前停下来主裂缝尖端不再向前扩展天然裂缝上的SDEG从交点向两侧扩展并且CSDMG演化区域明显比交点处宽。这时候孔隙压力云图会出现一个明显变化压力从主裂缝尖端“拐弯”进入了天然裂缝流体前锋沿天然裂缝继续推进。如果相交模式是偏移特征最微妙天然裂缝界面局部有损伤但并没有形成像捕获模式那样的大范围张开水力裂缝在天然裂缝另一侧某个偏移距离重新起裂整体路径出现一个台阶。这种模式最容易和穿越误判需要把SDEG、缝宽和孔隙压力放在同一张图里对比。5.3 缝宽和压力前沿怎么输出缝宽是判断相交模式的重要辅助量但Abaqus不会直接给一个叫“缝宽”的输出。cohesive单元完全破坏后缝宽近似等于上壁面和下壁面节点之间的法向位移差。在后处理里我一般把cohesive上下表面节点位移提取出来做一个简单的坐标差就是实际裂缝张开宽度。孔隙压力前沿的判断则要小心。裂缝内的孔隙压力应该比远场初始孔隙压力高而且沿着裂缝有压力梯度如果在cohesive单元上出现压力跳变或者局部负压先检查是不是单元方向、初始孔隙压力或者单位制出了问题。根本不是参数问题而往往是孔压初始条件没设对。5.4 结果不合理时按什么顺序排查如果模拟结果和现场微地震、井口压力曲线完全对不上我通常按固定顺序排查而不是乱调参数先看单位制这个Bug最容易藏在角落里再看初始应力平衡cohesive是否提前损伤第三看网格尺寸是否超过cohesive过程区长度第四看相交点拓扑第五看黏性正则化系数是否过大最后才调天然裂缝的黏聚强度和断裂能。这个顺序的好处是前面任何一个没检查好调后面的参数都没用反而越调越玄。我个人在实际项目里最深的体会是Abaqus里cohesive相交模拟的问题90%都出在参数和网格之前的“关系”上而不是求解器本身。把天然裂缝界面看成有摩擦、有黏聚强度、有断裂能的真实弱面把水力裂缝裂尖看成有损伤区的受力对象模型才能给出可信的相交模式。后面再做随机多裂缝扩展、离散裂缝网络压裂这类复杂问题也都是在这个基础上叠加和扩展。