做COMSOL水力压裂模拟这件事我前前后后折腾了大半年从照着案例库一步步点到后来能自己搭一套固体力学加达西定理的耦合模型中间踩过的坑比想象的要多。这个题目听起来像是教科书目录但真正做出来的东西是要回答一个非常现实的问题把高压液体注入到地层里孔隙压力怎么涨、有效应力怎么降、岩石骨架往哪个方向变形以及这个应力场的变化会不会达到岩石破坏的条件。适用的人群也很明确做储层改造的工程师、研究岩石力学的学生以及所有想把多物理场耦合从理论变成可计算结果的人。1. 水力压裂模拟的整体思路与物理场耦合1.1 为什么选COMSOL做这件事模拟水力压裂市面上不是没有专用软件。专门做裂缝扩展的商业软件确实很多它们在后端的断裂力学处理上很成熟但大多是个黑盒子你输入参数它给你输出裂缝形态中间过程看不清楚。而COMSOL的优势恰恰在于它把物理过程拆开了你能看到孔压场怎么演化、应力场怎么重分布、位移场怎么响应并且固体力学和达西定律这两个接口都是内置的。这一点在工程分析中特别有用。现场压裂施工前的设计阶段往往需要回答的不是裂缝能扩展多长这种终极问题而是更基础的这个注入压力下井筒周围的有效应力是多少孔隙压力扩散有多快会不会造成压穿或者滤失异常这些问题本质上就是固体力学和渗流力学的耦合问题用COMSOL恰好能给出透明、可检查的过程结果。从我个人的经验看COMSOL的价值还在于几何适应性。现场井筒不是点储层有层理、有非均质性这些都可以在几何里体现。虽然我们一开始建议从二维平面模型做起但COMSOL的架构允许你在一个模型里逐步增加复杂度而不用换工具。这一点对于需要反复试算的项目来说效率提升非常可观。1.2 固体力学和达西定理是怎么耦合的很多人把固体力学达西定理理解成两个物理场各算各的最后把结果拼在一起。这是最大的误解。水力压裂的模拟关键不是算得准不准而是耦合机制有没有建立正确。耦合的核心是有效应力原理。在油气藏岩石这种多孔介质里外部总应力由两部分分担一部分压在岩石骨架上叫有效应力另一部分由孔隙流体承担叫孔隙压力。我在模型里设置的就是这样把孔压场算出来以后通过有效应力关系作用到固体力学的应力平衡方程中反过来固体的变形会改变孔隙体积和渗透率进而影响流动。这种相互作用的数学表达就是Biot孔隙弹性理论。实际在COMSOL中操作时耦合有两种程度。一种是单向耦合孔压作为体载荷加到固体力学里面岩石变形算完就结束了不回头影响流动。这种模式适合初步分析速度很快但忽略了应力敏感效应。另一种是双向耦合需要自己定义渗透率随应力或应变的变化关系比如经典的指数衰减模型。我在做实际项目时通常先跑单向耦合验证模型稳定性再加上双向耦合做精细分析这是个很实用的策略。需要注意的是COMSOL中的多物理场节点并不是自动把所有物理场都耦合起来。你要在多物理场部分明确添加耦合特征把两个接口关联起来。我在初学时就犯过错两个物理场都算出了结果但没有建立耦合节点结果就是各算各的白白浪费了几天时间。这算是第一个要认真提醒的点。1.3 模拟边界这一模型能回答什么问题必须诚实地划一条线。纯固体力学加达理定理的连续介质模型本身并不能模拟裂缝的扩展路径和缝宽演化。真正那种裂缝从井筒延伸几十米甚至上百米的动态过程需要断裂力学、扩展有限元或者相场法来处理。那完全是另一个复杂度级别。那么这套模型的意义在哪里在我看来是工程设计和机理分析的基础层。它可以回答这几个问题注入期间井周应力分布怎么变化哪一个位置最容易先达到破坏条件。压力扩散的速率和范围受哪些参数控制。不同注入排量、不同初始地应力条件下的孔压累积规律。施工参数突变时有效应力的响应趋势。这些正是水力压裂设计里压不开压窜砂堵等工程问题的力学本质。在这个基础上如果后续要加入裂缝扩展也可以把这里的应力场和孔压场作为初始条件导入衔接是很自然的。2. 模型构建从几何到参数的完整细节2.1 几何简化的思路我建议第一次尝试时直接用二维平面应变模型不要一上来就搞三维。一方面现场尺度的三维模型网格量极易失控动辄上百万自由度调试周期很长另一方面平面应变模型已经能保留水力压裂的核心物理机制做参数敏感性分析完全够用。在实际建模中我常用一个方形区域代表储层岩石边长取200米左右。这个范围在给定渗透率条件下能够保证模拟时间内孔压扩散不触碰到边界避免边界效应干扰结果。井筒放在区域中心可以直接用一个圆孔代表井眼半径给0.075米到0.1米这样能规避点源奇异性的问题后面我会详细说。如果你关注的是单一裂缝起裂问题可以考虑利用对称性只取四分之一模型在对称面上设定合适的边界条件能显著减少计算量。几何简化不是越简单越好也不是越复杂越真。我看过很多失败模型问题都不是物理场配置而是几何上把一个200米尺度的区域和一个0.1米尺度的井筒放在一起过渡网格不好处理。所以几何设计时你一定要想清楚我的关注区域在哪里加密区放多大范围这样才能给后面网格划分留出合理的操作空间。2.2 材料与流体参数怎么定参数用错是这类模拟里最常见的问题。数值模型的结果再好输入参数不符合实际一切都白搭。以下是我在某致密砂岩储层的模拟项目中用过的参考参数范围分享出来供大家作为起点。注意这只是基准值不同区块的岩心测试数据可能差出一个数量级。参数代表值说明弹性模量2040 GPa致密岩石偏高软岩偏低泊松比0.180.30页岩常见0.20.25渗透率1e-171e-19 m²换算约0.011 mD孔隙度0.030.10储层越低渗孔隙度越低流体粘度0.52 mPa·s压裂液通常按清水近似Biot系数0.70.9岩石骨架可压缩性决定材料参数上有一点很难确定弹性模量和渗透率在储层纵向上差异可能很大。如果手头没有测井解释数据我的建议是先用单一均质参数跑通流程再做参数的敏感性扫描。不要一上来就试图把所有非均质性都装进模型里那只会让你连收敛问题都无法定位是在哪个环节出的错。流体的压缩系数也是个经常被忽略的参数。压裂液虽然是不可压缩流体但含气或者泡沫压裂液的可压缩性不能忽略。在达西定理这个物理场里你需要设置流体的密度随压力变化的规律最简单的做法是给一个较小的压缩系数比如4e-10 /Pa代表水的弹性。忽略这一点会导致孔压传播速度被显著高估算出来的扩散云图会偏大。2.3 初始地应力场怎么处理初始地应力场处理不当是这批模拟中最隐蔽的坑。很多新手直接从默认零初始应力开始计算结果得到的位移场和应力场都是围绕一个没有初始内力的状态展开的这跟真实地层条件完全不符。正确做法是分两步。第一步先解一个稳态问题只施加远场地应力边界条件不激活注入过程让模型得到初始应力场和初始孔压场的平衡解。第二步把稳态解作为瞬态计算的初始条件再开始注入模拟。这一步在COMSOL里的操作是先运行稳态研究步然后在瞬态研究的设置中调用前一步的解作为初值。初始地应力的量级怎么定取决于构造背景。我常用的方案是设最大水平主应力60兆帕、最小水平主应力45兆帕、初始孔压25兆帕。这三个数字不是随便填的它决定了井壁处的应力集中程度。如果你用的是不同区块的测井数据务必通过一个小型的初始应力平衡检查来验证在没有注入的情况下模型应该处于完全静止状态位移场应该是零应力场应该等于你施加的初始值。如果稳态求解后发现位移不为零说明边界施加和材料参数之间有冲突需要先调整。还有一个细节如果初始孔压和地应力是分别设置的要确认有效应力在数值上是合理的。有效应力等于总应力减去孔压在储层深度条件下应该为正值否则岩石处于拉伸状态不合常理。我在一次项目里发现某区域的等效拉应力超过了岩石抗拉强度排查后发现是初始孔压给高了修正后就正常了。3. 边界条件与注入过程的设置要点3.1 位移边界与应力边界怎么选固体力学模块的边界条件有两种类型很多人选错是因为没搞清楚两者的区别。应力边界条件等价于远场地应力加载位移边界条件则是固定边界节点的位置。在水力压裂模拟这种涉及现场尺度的问题上我通常建议施加应力边界而不是位移边界。原因很简单地层在构造应力作用下是有初始内力的你要模拟的是这个有内力的系统受到局部扰动后的响应。如果四个边界都用位移固定那么初始应力只能靠初始应力变量来指定算出来的位移场是相对于初始状态的增量相当于人为把系统钉死了。而如果在边界上直接施加远场应力让岩石在边界约束下达到应力平衡那模拟的就是真实的原地应力状态。当然位移边界条件也有它的应用场景。比如你要模拟多条平行裂缝之间的应力阴影效应用周期性的对称条件约束会更方便。但常规的单井压裂分析我推荐的外边界组合是两侧边界施加最大和最小水平主应力上下边界也按地应力梯度施加角点附近做适当约束以避免刚体位移。小技巧在地应力加载时不要把应力直接法向加载到一个倾斜的边界上。COMSOL中的边界载荷是沿边界法向的而地应力的主方向是全局坐标系方向。你需要用组件耦合或者边界载荷的分量表达式来按全局方向施加应力否则算出的初始应力场会发生旋转后期应力云图会非常奇怪。3.2 注液点的处理方式注液条件的设置直接决定孔压场的形态。最粗暴的做法是在井筒中心设置一个点源项单位时间的注入量作为一个常数。但这种做法在数值上会带来奇异性点源处的孔压会趋于无穷大导致该位置的应力计算完全失真。我的建议是不要让注入量作用在一个数学点上至少要放在一个有限尺寸的孔洞边界上。用圆孔代表井筒时注入条件可以在井筒边界上设置为流量边界条件流量大小按照压裂施工的排量换算。换算关系不难完井施工中常见排量是每分钟2到3立方米折算到二维平面模型里就是除以储层的厚度。比如储层厚度按30米算那么二维模型的边界流量就是2.5除以30再除以井筒周长换算成单位面积的通量。还有一种处理方式是直接在井筒边界上给定压力边界条件模拟的是井底压力恒定的工况。这在研究不同施工压力下有效应力分布这类问题时非常好用。工程实际中压裂泵注通常是排量控制但模型调试阶段用压力控制更容易收敛。我的经验是先跑压力控制模型拿到稳定的初场再切成排量控制模式做精度更高的分析。要注意注入速率不是越接近现场真实值越好特别是在二维平面应变简化下。如果换算后的实际通量超过画网格可承载的范围第一迭代步就可能让牛顿迭代发散。此时可以适当降低注入速率观察孔压演化的相对规律再逐步恢复到真实值。3.3 多物理场耦合节点的配置顺序COMSOL里建立耦合通常是在多物理场中点击添加耦合然后选择耦合类型。这里有个非常实际的配置顺序问题。正确的逻辑是先分别独立地把两个物理场的方程和边界条件设置好检查各自单独计算都没有语法或设置错误之后再去添加耦合。不要上来就建立耦合否则一旦算不下去了你根本没法判断问题是出在固体那边、达西那边还是耦合这边。添加耦合时的内容是两个方向。第一个方向达西定律算出的孔压p需要通过体载荷作用到固体力学方程中。在COMSOL中这个作用体现在有效应力上即在固体力学中把孔压定义为来自达西定理的变量并在应力平衡中减去孔压对总应力的贡献。第二个方向如果需要考虑应力敏感效应就要把渗透率定义为孔隙应变或者体应变的函数然后让达西定律使用这个修正后的渗透率。我在第一次搭双向耦合模型时把渗透率表达式写在了一处全局定义里结果发现达西定律并没有自动去调这个更新后的渗透率因为全局定义在物理场初始化时已经被快照了。后来我改成在达西定理的物理场节点内部直接输入渗透率表达式问题就解决了。这个细节在文档里很难找到卡了我两整天的时间。4. 网格、求解与后处理实操全流程4.1 网格划分策略网格是整个模型里最需要手感的部分。水力压裂模拟的问题在于井筒尺寸和储层尺寸相差数百倍如果网格均匀布置要么近场加密不够导致应力集中丢失要么整体加密导致自由度爆炸。我常用的划分策略是以井筒为中心设置一个半径大约20米的加密区域用三角形网格单元尺寸控制在0.5米以内加密区外围过渡到中等网格单元尺寸逐渐放大到5米在200米模型边界附近用10米左右的粗网格。这样处理之后整个二维模型的单元数通常在数万量级求解时间只有几分钟调试效率非常高。网格密度的重要性可以通过一个简单例子说明井壁处的环向应力集中依赖于裂缝或者圆孔边界的离散精度粗网格会把应力峰值明显低估导致你得出起裂压力远低于实际值的错误结论。我在一次参数敏感性分析中就吃过这个亏网格从2米加密到0.5米后井壁处的最大主应力数值变化了将近15%这个量级的误差足以改变工程判断。解决方法是做一个网格无关性验证。把网格逐步加密观察目标物理量比如井壁处最大主应力或者孔压扩散半径是否收敛。如果加密到某个程度后数值变化小于2%到3%就说明网格已经足够。这个方法虽然枯燥但是避免后期被质疑模型合理性的最有效手段。4.2 求解器配置与收敛控制水流体注入过程是典型的瞬态问题。在COMSOL中瞬态求解器的时间步进设置不能随便用默认值因为注入初始阶段的孔压变化速度极快而后期扩散速度逐渐放缓用一个恒定步长要么前面步长太大导致不收敛要么后期步长太小让计算时间无谓拉长。我的做法是采用自适应步长初始步长设置在0.001小时左右最大步长限制在1小时。这样能保证注入开始那一瞬间的急剧变化被捕捉到。同时把非线性求解方法和容差设置结合起来求解器中的非线性方法用默认的牛顿法即可但要适当放宽容差到0.01或者0.05严格的话反而容易因为局部震荡造成迭代失败。一个有效的收敛辅助手段是打开辅助扫掠或者分阶段加载。也就是说把注入速率写成时间的函数前几分钟先线性爬坡到目标值而不是第一步就全量注入。这既符合现场施工中泵车启动排量逐渐增加的实际过程又能显著提高数值稳定性。我在多个模型里都用了这个方式收敛率从勉强过线提升到稳定收敛。如果算到某个时间点突然中断第一件事不是去调整求解器设置而是去查看这个时刻孔压场或者应力场的分布看看是不是某个区域已经出现了病态值。通常你会发现是应力集中区域的局部网格畸变过大或者渗透率模型输出了一个接近负数的值。物理场出问题导致的数值崩溃用纯粹的求解器设置是救不回来的。4.3 后处理读云图也要讲方法后处理很多人以为就是出一张漂亮的云图实际上中间门道不少。孔压云图要重点观察的是井筒附近的孔压梯度。在注入初期孔压的扩散前缘会呈现一个清晰的梯度带这个梯度带推进的速度就反映了储层的导压能力。如果梯度带推进太快说明渗透率设置偏大或者流体压缩性设置过小如果推进太慢则反过来检查。应力场结果的处理上有三个量要区分清楚总应力、有效应力、应力增量。很多人直接看总应力云图但水力压裂的破坏判断使用的是有效应力。在岩石力学中抗拉强度比较的也是有效应力。如果你们在分析中直接拿总应力做判断那结论很容易出错。位移场我通常以矢量图和变形标量图结合来看。要特别注意在COMSOL的默认显示中变形缩放因子可能被自动放大导致看起来非常夸张的位移。有些新手看到几毫米的位移被缩放到肉眼可见的程度就以为自己算出了米级的裂缝开度闹出笑话。正确的做法是查看具体的位移数值以数据为准。5. 常见问题与排查技巧实录5.1 不收敛的几类标准剧本不收敛是这类模拟的高频问题但多数情况都有固定的剧本。最典型的三种是一是初始网格太粗井壁处应力集中剧烈局部单元变形过大牛顿迭代无法在容差范围内找到解。二是初始条件不连续比如初始孔压场和初始应力场之间本身就不满足静力平衡导致求解器在第一时间就发现矛盾。三是时间步长设置过大特别是在注入起止的瞬态阶段前几步就跳过了关键的物理转折点。排查的顺序我建议是先缩小时间步不行再加密网格最后检查初始条件。反过来操作的话往往会耗费大量时间在错误的方向上。另外我强烈建议在两组不同网格、不同时间步设置下各跑一次对照确认结果一致这比在单个模型上调参更能说明问题的根源。还有一个不太容易想到的点材料参数的量级跨度。岩石的弹性模量可能是几十吉帕渗透率却是10的负18次方量级这个跨度会让矩阵的数值条件数变得极其糟糕。这种情况下可以尝试调整单位系统比如把长度单位从米改为毫米压力单位从帕改为兆帕让矩阵元素处于更接近的数值范围。这不是投机取巧而是有限元计算中非常实用的预处理技巧。5.2 结果异常如何排查如果你算出了结果但结果明显不合物理需要一套系统的排查思路。我遇到过最典型的情况是注入很长时间后孔压仍然只在井筒附近扩散半径只有几米这与现场认识完全不符。排查后发现是渗透率单位换算出错模型里实际使用的是平方毫米而不是平方米差了6个数量级。单位制在COMSOL中默认是国际单位但你自己写的表达式和外部数据的换算最容易出错。应力场的异常也常见。有一次我算出的最小主应力变成了负值也就是出现了拉应力状态而且在远离井筒的区域也这样这明显不符合自重应力场规律。排查后发现初始应力场输入时只设置了水平方向两个主应力忽略了垂直方向应力分量导致面内应力状态不完整模型自然无法得到合理平衡。补上垂直应力后分布恢复正常。还有一个经常被忽略的细节就是接触面或者界面处的不连续。如果你在几何里设置了多个域来代表不同岩性那么域之间的孔压和位移默认是连续的但应力可能不连续。很多人在查看结果时看到跨越岩层界面的应力跳变怀疑是数值问题其实这是材料参数差异导致的物理跳变不能强行抹平。理解这个才能正确解释模拟结果。另外我提醒大家注意一个容易误判的场景在井壁附近的拉伸应力集中区域如果拉伸应力已经接近甚至超过抗拉强度这个区域发生破坏只是一种可能性并不代表裂缝已经起裂。因为破坏之后应力要重新分布后续的裂缝扩展路径需要断裂力学才能描述。连续介质模型报出的超限区更应该理解为优先起裂区而不是裂缝形态。5.3 几个容易忽略的单位与量纲问题单位问题听起来是老生常谈但在实际模型评审中我几乎每隔一段时间就能看到一次。最常见的坑有三个渗透率单位、流体黏度单位、流量与通量的混淆。渗透率的国际标准单位是平方米但石油行业习惯用毫达西也就是1e-15平方米。做模拟时如果你从岩心报告里直接拿了一个几十毫达西的数值当作几十乘以10的负12次方写进去那结果就会大错特错。正确的做法是统一换算成国际单位制至少在模型注释中记录原始数值和换算系数。我习惯在做模型的第一天就建一个简单的换算表贴在模型文档的顶部避免两天后忘了当初填的是什么。流量与通量的混淆也很常见。你从施工设计中拿到的排量是每分钟若干立方米单位是立方米每分钟而达西定理边界条件需要的通量是立方米/秒/平方米。中间的换算必须除以井筒壁的面积和换算系数。我在3.2节给过一个例子这里要再强调一次换算错误导致的注入速率如果是真实值的几十倍井筒附近的孔压会在几秒内达到不合理的量级然后模型崩溃。如果你用的是自定义的表达式或者函数务必在每个表达式后面手动写上量纲检验注释。COMSOL本身会显示表达式的单位但如果你的表达式逻辑错误有时单位检查也发现不了比如把应变这么一个无量纲量写成了有单位的量而单位检查只会提醒单位不匹配不会帮你判断物理逻辑是否正确。我在调试经验中总结出一条任何自定义变量的物理量纲都要在模型文档中和它所代表的物理含义一一对应列表记录。另一个不常被提及但在多尺度模型里非常关键的点网格尺寸和渗透率导致的临界时间步限制。当渗透率很低、网格加密时达西流动的局部Courant条件会非常苛刻瞬态时间步长被迫降到极小值模型会陷入计算一天推进十秒的困境。如果遇到这种情况可以考虑先做稳态孔压场的某一时刻快照或者适当提高渗透率到研究敏感性所需的范围不必强行追求和现场完全一致的绝对扩散速度。写在最后的一点实战体会这套模型跑通之后我最大的感触是水力压裂模拟的前期建模工作量只占三成另外七成的精力都花在验证和排查上。尤其是耦合机制我必须反复提醒自己固体力学算出的变形和达西定律算出的孔压到底是通过哪条数学路径相互作用每一步都要有清晰的逻辑闭环。做模拟的目的不是为了生成一张漂亮的云图而是为了能回答工程现场的问题这个区块的施工压力窗口在哪里注入排量再提百分之二十会不会压穿隔层储层的应力敏感效应会让孔压扩散显著变慢吗。我建议你从最简单的二维平面模型入手把固体力学和达西定理的耦合机制摸清楚再去接更复杂的断裂扩展模块。先按这个思路把基础模型跑通你后面面对任何复杂地质条件都会有一种手上有图、心里有底的感觉。
