这篇继续聊基于Cohesive单元的二维水力压裂。上一篇我们把最基础的模型从几何搭建跑到了能出结果的阶段但很多朋友反馈说模型能跑是一回事动一个参数就崩又是另一回事。这篇我打算把Cohesive单元二维水力压裂背后那些为什么这么设参数的逻辑、实际建模时容易踩的坑、以及报错之后到底怎么排查一次性讲透。内容主要面向已经在做压裂仿真、但被参数和收敛问题反复折磨的工程师和学生如果你还完全没接触过水力压裂模拟建议先跑通一个最简单的算例再来读这篇效果会好很多。1. 为什么二维水力压裂建模绕不开Cohesive单元1.1 先分清三类主流思路做水力压裂数值模拟第一件事不是急着画网格而是先把思路理清楚。目前主流的求解思路可以粗略分成三类一类是纯连续介质配合强度准则裂缝靠单元失效来显式删除或开孔表达一类是扩展有限元XFEM思路裂缝可以在单元内部穿行不需要预设路径还有一类就是本文要讲的在预设裂缝路径上布置内聚力Cohesive单元把开裂过程抽象成上下两个开裂面之间的粘聚关系。我自己的经验是在二维水力压裂这个场景下Cohesive单元是工程出图、参数标定、跟实验对比综合性价比最高的一条路。尤其当你关心的是裂缝在给定射孔位置如何起裂、如何扩展、缝宽和缝内压力随注入时间怎么变化时它比另外两条路好用得多。三种思路的适用场景和典型短板可以这样对比技术路线适合解决的问题典型短板单元失效/删除大范围破裂、压碎区模拟删除单元会丢质量孔压路径不连续XFEM扩展有限元裂缝路径未知、随机扩展多裂缝交叉和流固耦合设置很繁琐Cohesive内聚力预设压裂路径、缝宽与压力耦合需要预先布置界面存在一定网格依赖性选择Cohesive还有一个实际原因多数水力压裂问题的关注点不是裂缝往哪条路走而是给定压裂路径后裂缝开度、缝内压力、扩展速度怎么演化。这种情况下我们完全可以把断裂路径预设在Cohesive层上把计算资源留给流固耦合本身。另外Cohesive单元的处理思路和室内实验标定路径比较一致用双悬臂梁或者巴西劈裂实验拿到的断裂参数可以直接作为输入这一点是XFEM和单元删除法很难比的。1.2 Cohesive单元本质上是把断裂过程区翻译成界面本构Cohesive单元和普通实体单元最大的区别在于实体单元描述的是材料内部一点的应力应变关系Cohesive单元描述的是两个开裂面之间单位面积的拉力与张开位移关系。你可以想象用双面胶粘两块木板再从中间撕开刚拉的时候胶层还没破坏需要给力才能拉开这个阶段近似线弹性当拉力超过胶层的强度胶层进入软化力还没降到零但裂缝已经在发展最后胶层完全脱粘拉力归零。Cohesive单元就是在数值上复现这个过程。但是水力压裂里的Cohesive单元不能只用上面这种力学本构。压裂液是要在裂缝里流动、传递压力的所以单元还必须额外具备孔压自由度。拿Abaqus来说普通实体断裂用的是COH2D4而水力压裂必须用COH2D4P这类带孔隙压力自由度的单元。可以这么理解一个带孔压的Cohesive单元同时扮演了两个角色一个是粘聚界面负责传递未开裂段的拉力和软化段的残余应力另一个是可变宽度的管道负责让流体沿着缝长方向流动并把流体压力施加到裂缝面上。后面所有参数设置其实都是围绕这两个角色展开的。2. 参数背后那些力学逻辑不搞懂很难调对2.1 牵引-分离本构起裂到开裂怎么走Cohesive本构基本上就三段线弹性段、损伤起始、损伤演化。线弹性段由法向刚度Knn和切向刚度Kss决定描述的是单元还没损伤时的力与位移关系损伤起始用最大名义应力准则或者二次名义应力准则来判断损伤演化用断裂能Gc或临界位移来定义决定了裂缝从起裂到完全脱粘这一段怎么走。在水力压裂场景裂缝以I型张开为主所以法向行为是最核心的。线弹性阶段法向应力t_n和法向张开位移δ_n之间满足t_n Knn * δ_n当t_n超过法向强度t_n0或者更准确地说当准则值CSMAXSCRT达到1时单元开始损伤之后应力按t_n (1-D) * Knn * δ_n逐步释放D是损伤变量从0变到1。D等于1时单元完全失去承载能力但注意这时候单元并没有被删除它仍然能传递流体压力这正是水力压裂模型需要保留单元的原因。很多初学者会误以为损伤之后单元就该消失这是从实体单元失效删除带过来的惯性思维。在Cohesive水力压裂里完全损伤的单元恰恰是裂缝的实体它要承担缝内流体的流动通道功能如果把它删了流体就没地方流了后面的压力计算全乱套。2.2 初始刚度、强度、断裂能之间如何互相制约新手最容易犯的错误是把初始刚度Knn、强度t_n0、断裂能Gc这三个参数当成互相独立的变量来调。实际上它们共同决定了裂缝的软硬程度和脆韧程度。断裂能Gc是牵引-分离曲线下包裹的面积如果采用线性软化模型临界张开位移δ_c就可以直接算出来大约是2Gc/t_n0。换句话说强度和断裂能给定了单元从起裂到完全脱粘能张开多少已经是定死的了。这三个参数的实际表现是t_n0越大、Gc越小软化段越陡峭开裂越脆数值上就越难收敛t_n0太小则起裂太早裂缝在流体压力还不够高的时候就开始软化Knn如果取得过大整体刚度矩阵条件数会恶化容易出现病态和零主元Knn如果取得过小裂缝两侧在未损伤时就有明显相对位移模型整体偏软和实体单元的变形不协调。工程上我一般用这个经验公式来估初始刚度Knn大约等于E除以裂缝过程区特征长度。假设岩石弹性模量E30000 MPa断裂过程区特征长度取0.1 mm量级Knn大约就是3e5 MPa/mm。实际模型里给到1e4到1e6 MPa/mm这个区间都是可以接受的关键是和实体单元模量匹配。抗拉强度可以参照岩石的抗拉强度2到6 MPa是比较常见的页岩、砂岩范围断裂能Gc则参考室内实验一般取0.05到0.2 N/mm对应50到200 J/m²。切向方向的强度和断裂能可以给大一些避免在纯I型工况下提前出现剪切损伤。2.3 孔压单元的流体方程切向流动、法向滤失与裂缝宽度带孔压的Cohesive单元流体运动分成两部分切向流动和法向流动。切向流动就是压裂液沿着缝长方向的流动它对裂缝扩展的贡献最大。切向流动的流量和裂缝宽度的三次方成正比近似满足立方定律q - (k_t * h³) / (12μ) * dp/dx。也就是说缝宽增加一倍流动能力就增加八倍所以一旦裂缝张开缝内压力分布会迅速趋于平缓这个特征和现场压裂压力曲线能对上。法向流动代表流体从裂缝面泄漏进地层的量也就是滤失。很多入门模型会把法向泄漏系数设成0表示地层致密、完全不漏这样计算简单但压力曲线会明显高于实际。建模时我建议第一阶段先设成0把模型跑平稳了再逐步加入泄漏项否则一开始就兼顾滤失的话模型调起来非常痛苦。还有一个容易被忽略的参数是初始缝隙。Abaqus里的Cohesive孔压单元切向流动中参与计算的缝宽是初始缝隙加上当前张开量即hc initial gap OPEN。如果初始缝隙设成0意味着单元只有在张开后才有流通能力这符合物理直觉但数值上会导致单元一旦开始张开流体突然涌入压力震荡。我给初始缝隙一般取0.001 mm到0.01 mm的量级它比最终缝宽小一个数量级以上对结果影响不大却能极大改善收敛性。2.4 单位换算最容易翻车的隐形坑Abaqus本身没有固定单位制用户必须保证全模型单位自洽。做水力压裂仿真我最推荐用mm、N、ton、s这套单位因为应力单位正好是MPa读数据非常直观。很多不收敛和压力数量级异常的问题最后都能回溯到单位错误。在水力压裂模拟里涉及流体参数的单位最容易乱。使用mm制时水的密度是1e-9 ton/mm³水的动力黏度是1e-3 Pa·s换算到mm制变成1e-9 MPa·s断裂能方面1 J/m²等于0.001 N/mm所以100 J/m²在软件里要填0.1 N/mm注入流量在二维模型里通常写作mm³/s如果你模型实际代表每延米或每米厚度要按出平面厚度做等量换算。这里每一项单独看都不难但它们混在一起时很容易出现压力高出几个数量级的诡异问题建议调参之前先把自己模型的单位制写在模型文件最显眼的位置。3. 二维水力压裂模型搭建实操3.1 几何处理和网格布置以Abaqus为例二维水力压裂的几何可以非常简洁一块矩形岩体中间从注入点向外延伸预置一条裂缝路径。裂缝路径两侧的网格镜像布置在路径位置插入一层Cohesive单元。这里有个分支抉择做零厚度Cohesive单元还是在几何上画一层极薄的实体单元。零厚度单元在建模时需要处理节点位置和单元偏移物理意义更接近真实裂缝薄层方法则是在几何上直接画一层比如0.0001 mm的薄层分配Cohesive截面时通过constitutive thickness参数把几何厚度归一化材料刚度不受实际薄层厚度影响。入门阶段我强烈建议先用薄层法把流程跑通后面再去折腾零厚度。薄层法最大的好处是网格生成简单不容易出现单元过度畸变数值稳定性好。网格尺寸方面裂缝路径上至少布置30到50个Cohesive单元注入点附近要更密因为那里的压力梯度最陡路径附近的实体单元以目标缝长的一半为尺度做过渡加密。网格太粗时裂缝扩展是跳跃式的每跳一格就会在压力-时间曲线上留下一个锯齿峰后处理阶段非常难受。3.2 材料参数与Cohesive截面定义材料参数这里给一份可供直接起步的参考表注意这是起步值具体项目必须通过敏感性分析微调对象参数推荐值岩石实体弹性模量E30000 MPa泊松比ν0.25密度2.7e-9 ton/mm³Cohesive层法向刚度Knn1e5~1e6 MPa/mm切向刚度Kss1e5~1e6 MPa/mm法向强度t_n03 MPa切向强度t_s010 MPa避免剪切先坏I型断裂能GIC0.1 N/mmII型断裂能GIIC0.5 N/mm混合模式BK指数1.2流体密度1e-9 ton/mm³动力黏度1e-9 MPa·s水注入流量5e-4 mm³/s量级初始缝隙initial gap0.001 mm在Abaqus里设置Cohesive材料时需要同时勾选力学Traction行为和Pore Fluid Flow行为并在Pore Fluid Flow里定义切向流动和法向泄漏参数。单元类型选择COH2D4P二维四节点带孔压Cohesive单元。有一点务必注意单元控制里的单元删除选项不要打开。一旦打开损伤后的单元会被删除孔压和流体路径随之丢失这对水力压裂模型是致命的。我见过不少人的模型裂缝走到一半压力全没了就是这个选项被默认勾选或者误开了。3.3 初始应力、孔压边界与分析步设置水力压裂模拟必须考虑初始地应力。如果完全不设初始应力裂缝在无围压情况下会过度张开注入点压力没有压裂层的参考价值。对二维剖面模型通常给一个垂向应力σv和一个水平最小主应力σh注意压应力在软件里用负值。初始应力可以用初始条件定义推荐在第一分析步之前加一个geostatic或平衡步让模型在注液前就达到初始平衡状态位移几乎为零。这一步做好的标志是没有任何载荷时模型位移场基本是干净的只有地应力引起的微小预变形。分析步选择上用Static, General并勾选孔压选项或者用专门的多孔介质流固耦合分析步。时间总长对应实际注液时间比如100秒初始增量步建议给到1e-5这个量级最大增量步长不要超过总时间的1/100最大增量步数要设大默认100往往不够我一般设5000以上。Cohesive材料里加一个小的粘性正则化系数1e-6到1e-5之间这个操作对软化段的收敛性提升非常明显而且只要系数不过分大对结果的扰动可以忽略。3.4 注入点处理和加载方式选择注入点通常选在Cohesive层的一端也就是射孔位置。施加的载荷形式可以是集中孔压流量也可以是固定孔压。流量加载更接近现场低排量起裂、高排量扩展的操作也是标定实验数据时最常见的加载方式。二维模型默认出平面方向是单位厚度流量单位是mm³/s如果模型代表的实际厚度是B毫米那要把总流量除以B换算成单位厚度流量再施加。我在项目中常用一个验证技巧先用流量控制注液让裂缝平稳起裂并有稳定扩展然后中途把加载方式切换成固定孔压观察两条路径下最终缝宽和缝长是否能对上。如果能对上说明模型对流固耦合的处理是稳健的如果差得离谱大概率是边界条件或者Cohesive截面设置有隐藏问题。这和压裂现场先看压力建立、再调排量的思路一致也是把数值模型向工程语义靠拢的一种方式。4. 常见报错与不收敛问题排查实录4.1 Negative eigenvalue和零主元不是一回事Abaqus的msg文件里出现Negative eigenvalue很多人的第一反应是模型要崩了。其实负特征值和零主元是两码事。零主元多半是约束不足比如模型忘了消除刚体位移或者某个部件没绑好负特征值则更多是切线刚度矩阵在软化段不再正定造成的Cohesive单元一旦进入损伤演化局部刚度下降出现负特征值是物理上正常的信号。判断标准很简单如果负特征值只是偶尔出现而且后续迭代还在推进、残差没有爆炸那基本可以忽略如果它持续出现位移增量疯狂增加那就是真不收敛了。处理顺序是先检查边界条件和约束再检查Cohesive初始刚度是不是过大导致尖端应力集中太剧烈最后考虑加粘性正则化。粘性正则化系数我建议调试阶段一开始就给1e-5等模型稳定了再逐步降低看结果是否对系数敏感。如果不敏感说明这个数值辅助手段没有污染结果可以放心用。4.2 裂缝起裂后立刻发散先把这三个嫌疑查一遍模型能算到起裂但损伤因子D一从0增长就立刻发散这是二维水力压裂里最常见的卡点。根据我的排查经验三个嫌疑最大第一增量步太大。起裂瞬间应力释放和流体重新分布非常剧烈相当于绷了很久的弦突然断了增量步太大会导致迭代直接飞掉。把初始增量步降到1e-6甚至更低很多情况下能直接跨过这个坎。第二材料参数组合太脆。强度和断裂能的比例关系决定了软化段斜率如果强度给得特别高、断裂能又特别小软化段几乎垂直任何数值扰动都会被放大。这时候把强度适当降低或者把断裂能提高一点收敛性会立刻改善。第三流体提前进入了未损伤的单元。Cohesive单元的切向流动基于总缝宽如果初始缝隙设得偏大或者渗透率系数设置不当流体可能在单元还没损伤时就已经串过去了形成虚拟裂缝。检查gap flow设定和损伤起始准则之间的配合确保流体主要在损伤后才大量进入。4.3 压力怎么都上不去或高得离谱模型能跑但压力曲线不符合预期这是第二阶段最常见的问题。压力上不去先查三件事孔压初始条件有没有给对特别是初始孔隙压力是不是0法向泄漏系数是不是设得太大流体全漏到地层里了注入流量是不是太小压裂液连裂缝都没填满。压力高得离谱第一嫌疑绝对是单位错误黏度少乘了1e3、流量多乘了1e6都是实际遇到过的其次看缝宽立方定律下缝宽小了流动阻力是三次方级别的上涨如果初始缝隙设得太小缝内压力会异常攀升。所以看到压力数量级不对先放弃调参数回头把所有单位列出来逐一核算。还有一种情况是压力曲线从头到尾都是平的既不升也不降。如果滤失为零、裂缝又没有扩展理论上压力就应该持续上升。曲线走平说明流体没有起到推进作用多半是在注入点附近形成了一个高压小腔裂缝本身没有有效张开。这时候用后处理看Cohesive单元的孔压分布如果压力只集中在注入点一两个单元附近基本就是这个原因。4.4 网格敏感性控制与验证Cohesive单元因为带了断裂能理论上对网格的敏感性比单元删除法小得多但流体驱动型裂缝依然受网格尺寸影响。裂缝尖端的流体压力和固体变形强耦合网格太粗会让尖端附近的压力梯度和张开量被平滑掉导致扩展偏慢。我做模型一直保留三套网格粗、中、细裂缝路径上单元数量分别约20、40、80个对比三条压力-时间曲线和裂缝半长曲线。如果粗细网格之间的差异小于5%就认为网格影响可以接受。网格太细也不是免费午餐。Cohesive单元数量上去了计算时间成倍增长而且在非常细的网格下每个单元的刚度、强度参数可能还需要重新微调并不是越细越好。更合理的做法是固定一套质量良好的网格然后只做材料参数敏感性分析这样对比出的趋势才有意义。如果一套粗网格怎么调参数都出不来合理结果那问题大概率不在参数而在几何或者边界条件本身继续加密网格只会浪费机时。4.5 too many attempts 的终极方案看到too many attempts increase的程度至少报错先别急着乱改参数。按顺序排查第一步打开msg文件看最后一步是在哪个增量步失败的是加载初期就失败还是裂缝扩展中期失败这决定了你该改加载步长还是材料参数第二步把最大增量步数从默认值调到5000以上很多模型其实是被步数上限硬生生卡死的第三步降低初始增量步长Abaqus自己会在收敛后自动放大但初始给太大它可能一开始就放弃了第四步加粘性正则化第五步检查单元畸变如果Cohesive单元被压扁或者拉成畸形考虑不加单元删除选项而是把网格调整一下。如果以上招都试过还不行那就做模型回退先做一个纯力学模型不加流体耦合验证裂缝在给定张开位移下能正常扩展再单独做一个无裂缝的孔压扩散模型验证流体方程能正常计算。这两步都过了再合在一起调试。这个分层调试法比瞎改参数有效得多至少能省掉一半的无头苍蝇时间。5. 后处理怎么读压力和裂缝宽度的正确打开方式5.1 该输出哪些变量Abaqus后处理里Cohesive单元主要看SDEG刚度退化量从0到1或者CSDMG它表示单元的损伤程度CSMAXSCRT是损伤起始准则的发挥值超过1说明已经满足起裂条件。裂缝张开量可以通过变形后的位移场来提取沿Cohesive层取上下两排节点的位移差值画出来就是缝宽沿缝长的分布曲线。注入点的孔压随时序可以直接输出这是和现场压力计对比的核心数据。我在实际项目里习惯同时输出单元积分点的孔压POR和节点位移U。因为Cohesive水力压裂的缝内压力分布应当是注入点最高、缝尖最低的形状如果压力沿整条缝几乎一样平说明流体的压力传播被过度简化了要么是缝宽过大让流动阻力变得可忽略要么是流体在未损伤阶段就已经贯通了整条缝。通过观察POR云图能快速判断模型是处于正常扩展还是出现了虚拟泄漏。5.2 为什么你的压力曲线和教科书对不上理想化的恒流量压裂压力曲线通常有三个特征段早期孔弹性阶段注入时间增加、压力近线性上升到达起裂压力后出现一个明显的下降尖峰之后是裂缝扩展阶段压力随裂缝增长而缓慢变化。但数值模型要完美复现这段曲线其实不容易最常见的原因有三个。第一初始缝隙设得太大。模型相当于一开始就有一条微裂缝早期压力积累被吞掉了第一个特征段完全消失。第二地应力设得太小。裂缝很快就进入扩展压力曲线直接只留下起裂后的平台段。第三滤失设为零。无滤失假设会让压力偏高这不一定是模型错误而是边界条件不同但和现场数据对标时要有心理预期。对照实验数据时建议先看三个特征起裂压力的绝对值、起裂出现的时间、扩展段的压力走形。哪个和实验对不上再去针对对应的参数不要一上来就追求全局完美匹配。5.3 一条快速自检的验收路径二维单缝模型最直接的验收手段是拿KGD类解析解或者文献实验数据来对标。经验做法是固定注入流量统计若干时刻的裂缝半长和注入点压力检查裂缝半长是否满足时间幂律关系、压力曲线是否在合理范围内。如果没有现成解析解也可以做网格无关性验证和参数敏感性分析来证明模型的可信度。我自己的验收标准是用中等网格跑出来的压力曲线和裂缝扩展趋势在换到细网格时不能有质的改变改变断裂能20%裂缝长度应该有可重复的响应趋势注入流量加倍起裂时间要相应缩短。如果这些基本规律都不满足那说明模型结构本身有问题后处理云图画得再漂亮也只是自嗨。到这一步模型从能跑到可靠才算真正走完。6. 写在后面几个我一直保留的调试习惯最后分享几个我做了大量Cohesive水力压裂算例之后留下来的习惯。第一开始一个新模型从来不直接全参数灌进去。我会先用单单元或者小模型验证Cohesive的牵引-分离行为一条应力-张开位移曲线拉出来看它和理论值对不对得上这一步只需要几分钟却能避免后面无数次瞎调参数。第二调参时每次只动一个变量并且把对应的压力曲线、裂缝长度、损伤云图截图存档。时间一长你对这套模型的参数敏感性会形成直觉哪种参数会让起裂压力升高、哪种会让扩展变慢一眼就能判断。第三别急着把所有复杂机制都塞进模型。滤失、非牛顿流体、多裂缝交错这些机制每一项都会显著增加收敛难度。先把单缝、牛顿流体、无滤失的基准算例做到和解析解或者文献曲线对得上再一步步往上加这样即使后面出问题你也知道该去哪里找原因。水力压裂模拟是个参数魔鬼藏在细节里的领域能把基础算例做干净比会堆叠一堆高级选项有用得多。
)
