做压裂模拟的人应该都遇到过这种场景:水力裂缝扩展几十步都好好的,一到天然裂缝附近,要么裂缝停在界面上一动不动,要么沿着天然裂缝疯狂拐弯,更气人的是模型有时直接不收敛。这个问题落到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参数的一个换算例子
我用一个具体数据串一遍。假设岩石弹性模量E=30GPa,泊松比ν=0.25,I型断裂韧性K_IC=1.5MPa√m。平面应变条件下断裂能就是:
G_IC = (1-ν²) × K_IC² / E
代入:0.9375 × 2.25e12 / 30e9,约等于70.3N/m。
如果黏聚强度T0取5MPa,线性软化下的最终分离位移δf=2G_IC/T0≈28微米。这个值很小,说明CZM区域很薄。为了避免数值问题,初始刚度通常会让损伤起始位移取δf的1/10左右,也就是2.8微米,K0=T0/δ0≈1.78e12Pa/m。II型断裂能G_IIc一般比G_IC高,岩石类材料可以取180到300N/m,混合模式用BK准则,幂指数取1.5左右比较常见。
Abaqus材料卡示意如下:
*MATERIAL, NAME=HF_CZM *ELASTIC, TYPE=TRACTION 1.78e12, 1.78e12, 7.12e11 *DAMAGE INITIATION, CRITERION=MAXS 5.0e6, 5.0e6, 1.0e6 *DAMAGE EVOLUTION, TYPE=ENERGY, MIXED MODE BEHAVIOR=BK, POWER=1.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靠区域分解。很多人设置cpus=16以为一定比cpus=8快,实际上当模型单元数不够大时,MPI通信开销早就把计算收益吃掉了。
我常用的提交命令类似:
abaqus job=frac_model cpus=8 double=both output_precision=single interdouble=both会让所有计算结果用双精度,对cohesive这种小位移、大刚度问题很有帮助,但速度会变慢,输出精度反而不需要双精度,所以output_precision保持single即可。在Linux下,我习惯用nohup后台提交,避免ssh断开导致任务中断:
nohup abaqus job=frac_model cpus=8 double=both > run.log 2>&1 &关于mpi threads,我个人的体会是:单机多核不需要刻意调MPI线程数,直接给cpus就行。跨节点并行才需要仔细配置MPI环境和主机列表。如果你发现cpus从8加到16几乎没加速,甚至更慢,多半是核数超过内存带宽或物理核心数。
4.4 常见报错和对应的处理思路
下面这几个报错是cohesive水力压裂模型里出现频率最高的:
| 报错现象 | 常见原因 | 处理手段 |
|---|---|---|
| Negative eigenvalue | cohesive单元软化段进入负刚度 | 调整*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和CSDMG
Abaqus后处理里,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%都出在参数和网格之前的“关系”上,而不是求解器本身。把天然裂缝界面看成有摩擦、有黏聚强度、有断裂能的真实弱面,把水力裂缝裂尖看成有损伤区的受力对象,模型才能给出可信的相交模式。后面再做随机多裂缝扩展、离散裂缝网络压裂这类复杂问题,也都是在这个基础上叠加和扩展。