做煤层气强化开采的数值仿真最近几年几乎绕不开COMSOL Multiphysics。尤其是注气增采这一块单纯用黑箱模型已经很难满足现场需求了大家都开始往THM耦合这个方向走也就是热-流-力三场耦合。这个项目就是基于COMSOL平台把煤层注气过程中的温度场、渗流场和应力场统一装进一个模型里去评估注气对甲烷解吸、运移和产出的综合影响。如果你正在做煤层气、页岩气或者致密砂岩气的开发方案论证或者你手头有COMSOL但不知道THM耦合怎么做这篇内容可以给你一条相对完整的参考路径。先说清楚这套思路能解决什么问题现场注气强化甲烷开采本质上不是一个单一的渗流问题。注入气体比如氮气或者二氧化碳会改变孔隙压力孔隙压力变化会引起有效应力调整应力调整又会影响裂隙开度和渗透率渗透率变化反过来改变气体流动路径。如果注入的是热气或者注气后伴有温度变化温度场还会通过热膨胀、吸附/解吸热效应进一步扰动这个系统。这就是THM耦合要描述的核心场景。COMSOL在这里的优势很明显它不需要你从底层写有限元代码物理场接口和耦合项都内置了你要做的是把地质模型、边界条件和耦合机制配置对。下面我从项目思路、机理方程、建模实操到参数标定一条条展开讲。1. 为什么在这个时间点做煤层注气THM仿真1.1 煤层气强化开采的现实难题煤层气开采和常规天然气差别很大。甲烷在煤里主要以吸附态存在于微孔表面只有不到两成是游离气。常规排水降压能解吸一部分但越往后解吸速率越慢单井产量衰减很快。很多矿区的井打到中后期日产气连500方都保不住这时候如果还想提升采收率注气是公认的方向之一。原理上不复杂注入二氧化碳或者氮气可以降低甲烷在煤表面的分压通过竞争吸附或者浓度驱替把甲烷置换出来。但问题在于注气不是简单的“往里灌就行”注入气在裂隙网络中的波及范围、压力传播速度、对煤体渗透率的动态影响这些都会决定增采效果。更麻烦的是煤层本身是强应力敏感性的多孔介质。注气抬升孔隙压力后有效应力降低裂隙会张开渗透率随之增大这是有利的一面。但压力升高到一定程度后可能诱发局部塑性变形甚至破坏裂隙又会起变化。温度场参与进来以后情况更复杂煤对温度很敏感升温能促进解吸但同时热应力会改变裂隙开度。这一连串的耦合效应靠室内实验一条条测出来不现实靠现场试错成本太高数值仿真就成了重要的评估手段。1.2 为什么COMSOL适合干这个活COMSOL在THM耦合建模上优势非常明显。它的多物理场耦合机制不是简单的“数据传递”而是在同一个求解框架里把各个物理场的方程联立在一起这对强耦合问题很友好。煤层的THM问题恰恰属于强耦合温度变化影响应力应力变化影响渗透率渗透率变化影响压力分布压力变化又反馈到有效应力和吸附平衡这种多向交互用弱耦合方式计算很容易发散或者失真。再一个就是自定义灵活。煤储层和常规岩层不太一样它需要考虑Langmuir吸附、基质收缩、裂隙压缩等特殊机制。COMSOL里虽然预置了达西流动、多孔介质传热、固体力学这些接口但真正的煤层气模型通常需要在材料属性、源项或者耦合变量里加入自定义公式这一点COMSOL做得很顺手。项目里我把渗透率随应力的变化、气体吸附引起的煤基质变形都做成了自定义表达式整体集成度很高后处理出图也方便给甲方汇报的时候可视化效果比Fortran程序跑出来的结果直观得多。另外还要说一句THM耦合并不仅限于煤层气。多年冻土路基、地热开采、核废料处置这些都是典型的THM场景COMSOL在这些领域已经有大量案例可以借鉴。模型框架搭好之后改材料参数和边界条件换一个工程场景并不难这也是我一开始就选COMSOL而不是其他专业油藏模拟器的核心原因。2. THM耦合机理三场如何互相咬合2.1 热-流-力三场的物理图像理解THM耦合先得有清晰的物理图像。你可以把煤体想象成一块干燥的海绵里面浸满了水但这里的水换成了压缩气体海绵内部的骨架还会随着温度和压力的变化自己变形。温度场负责描述这块“海绵”内部温度的分布和传递热源来自注入气体的温度差异还来自吸附/解吸产生的热量。渗流场负责描述气体和水在裂隙网络里的流动流动的驱动力是压力梯度流动的阻力取决于渗透率和流体黏度。应力场负责描述煤骨架在外力和内部压力作用下的变形与破坏受力状态改变以后裂隙的张开程度会发生变化进而改变渗透特性。这三场不是各自独立跑一遍再拼在一起而是每一步都在互相影响。注气后压力升高孔隙压力变化直接改变有效应力有效应力变化改变裂隙宽度裂隙宽度改变渗透率渗透率改变流量和压力分布。温度升高时煤基质受热膨胀但裂隙边缘的煤块膨胀方向不同会挤压或者拉开裂隙面同时温度升高会加速解吸解吸又带走一部分热量。这种循环依赖关系就是THM耦合的本质。2.2 控制方程与耦合路径在COMSOL里搭THM模型底层方程并不需要你从零推但你必须理解每一项代表什么否则参数填错了很难排查。我这里把关键的方程逻辑梳理一下。渗流场的基础是质量守恒方程COMSOL的达西接口一般写作ρ (∂(ε_p p)/∂t) ∇·(ρ u) Q_m其中u是达西速度ε_p是孔隙度Q_m是质量源项。在煤层气问题里源项不只是注入井的流量还包括甲烷从吸附态变成游离态的产出项这个需要用Langmuir吸附动力学或者平衡解吸模型来描述。动量层面用达西定律u -(k/μ) ∇p这里的渗透率k不是常数是随应力和温度动态变化的这就要和固体力学场联立起来。常用的经验公式是渗透率随有效应力的指数关系k k0 exp(-α(σ_eff - σ_eff0))α是应力敏感性系数煤层的α值通常在0.1~0.6 MPa⁻¹之间具体取决于煤阶和裂隙发育程度。固体力学场用的是孔隙弹性方程COMSOL内置的多孔介质弹性接口已经包含了孔隙压力引起的附加体积力项。对煤层这种裂隙发育的介质还需要在材料本构里加入塑性选项煤在拉压应力超过屈服点后会发生不可恢复的变形这时候需要定义屈服准则和硬化参数。COMSOL里常用Drucker-Prager准则来描述煤岩的塑性行为需要输入的参数包括内聚力、内摩擦角这些可以从煤岩力学实验中获取。温度场是经典的多孔介质传热方程(ρC_p)_eff ∂T/∂t ρ_f C_pf u·∇T ∇·(-k_eff ∇T) Q_T这里(ρC_p)_eff是骨架和流体的等效体积热容k_eff是等效导热系数。Q_T除了外部热源外还要把吸附/解吸热效应加进去煤对甲烷的吸附是一个放热过程解吸吸热这个热效应在注气驱替过程中不能忽略。三场的耦合路径归纳起来有四条孔隙压力变化通过有效应力原理影响应力场应力应变通过渗透率模型反馈到渗流场温度变化通过热膨胀影响应力场吸附/解吸通过源项和质量守恒影响渗流场同时对温度场产生热效应。COMSOL里做耦合时这些交叉项都在各自的系数里关联了所以你不用写耦合代码但要保证各物理场的变量名一致表达式引用正确。2.3 关键物性参数怎么定THM模型对参数极其敏感参数定不准模型再精细也是白搭。我做项目时总结了一套参数获取优先级煤岩力学参数优于热物性参数热物性优于流动参数流动参数优于解吸参数。力学参数里弹性模量、泊松比、内摩擦角、内聚力是必须的。弹性模量煤层一般1~4 GPa比普通砂岩低一个数量级泊松比0.25~0.4内摩擦角20°~40°内聚力0.5~2 MPa。热物性参数方面煤的导热系数0.2~0.5 W/(m·K)比热容1200~1700 J/(kg·K)热膨胀系数10⁻⁵ 1/K量级。流动参数的核心是初始渗透率和孔隙度裂隙煤层的渗透率跨度很大从0.01 mD到几十mD都见过要根据试井数据定。解吸参数用Langmuir体积和Langmuir压力描述这个通过等温吸附实验获取不同煤阶差异很大低阶煤的Langmuir体积能到30 m³/t以上高阶煤可能只有15 m³/t左右。参数太多的时候建议先做一个敏感性分析把影响最大的几个参数标定准其他用经验值就行。这个后面章节细说。3. COMSOL建模全流程从几何到求解器3.1 几何建模与网格划分COMSOL的几何建模思路和CAD软件不一样它更强调“够用就好”。煤层注气模型不需要把每个裂隙都画出来那是离散裂隙网络的事在宏观连续介质模型里裂隙的影响是通过渗透率张量体现的。我建议第一步先建立二维剖面的简化几何比如一个500 m × 50 m的长条形区域中间设置一口注气井两端设置生产井。二维模型跑通以后再扩展成三维这样可以大幅减少前期调试时间。几何建立有两种路径。一种是在COMSOL内直接画矩形、圆形再布尔运算适合规则几何另一种是从外部CAD导入适合实际矿区边界。外部导入时要注意COMSOL的CAD内核和外部软件之间的拓扑兼容性问题我在项目里就遇到过导入的复杂曲面直接被判定为不支持的拓扑最后只能用简化边界替代。这个问题的根源是外部CAD模型里存在微小缝隙或者重叠面COMSOL转换时无法自动修复处理方法是在外部CAD软件里先做几何清理删除小面、合并重合边再以STEP格式导入成功率会高很多。网格划分方面煤层模型需要重点考虑裂隙流动的特征。注气井附近压力梯度和速度梯度大必须加密生产井附近同理。流场变化平缓的区域可以适当放宽。我常用的是三角形网格加边界层在井筒圆周附近设置8~10层边界层网格用来捕捉近井地带的压力漏斗。网格数量控制在10万到50万之间比较合适太多了求解时间成倍增加太少了精度不够。这里要给个建议先跑一个粗网格模型看趋势再加密看关键位置的解有没有明显变化这叫网格无关性验证做完了才能放心用结果说话。3.2 物理场接口与耦合节点物理场接口选择上我的推荐组合是这样固体力学接口负责应力应变计算达西接口负责气体和水的渗流。这里要说一下COMSOL里有两种流动接口可选达西接口和多孔介质流动接口。达西接口更简洁适合单相或者两相简化问题多孔介质流动接口带相传递和毛细压力更贴近真实两相流。我在初始建模时用达西接口单相气体流动先把耦合逻辑跑通再根据需要升级为两相模型。传热用多孔介质传热接口它会自动计算固体和流体的等效热容和导热系数你只需要分别输入煤骨架和流体的热物性参数。如果注气温度和环境温度差异不大很多人会忽略温度场但实际工程中如果采用热气注入或者注气速率很高导致焦耳-汤姆逊效应明显温度场就必须保留。耦合节点的设置是COMSOL里最容易出错的地方。热-力耦合用“多孔介质热弹性”特征这个特征会自动把热应变和孔隙压力引起的应变加入应力方程你只需要激活多孔弹性材料模型并填入热膨胀系数和毕奥系数。流-力耦合主要通过“达西-固体力学耦合”节点实现它会在固体力学方程中加入孔隙压力引起的体积力项同时在达西方程中考虑体积应变引起的储集变化。吸附引起的煤基质变形和渗透率变化需要用变量表达式写入。我在模型里定义了一个渗透率修正因子表达式是基质收缩项和裂隙压缩项的乘积核心思路是当压力降低时基质解吸收缩导致裂隙增大渗透率上升同时有效应力增加导致裂隙压缩两者互相竞争。很多文献里用的Palmer-Mansoori模型就是这个思路你可以在COMSOL里用变量节点把它公式化。3.3 求解器配置与收敛控制求解器选型是COMSOL项目里最需要耐心的一步。THM耦合问题属于高度非线性问题求解器配置不当第一步迭代就发散毫不夸张。我的一般策略先做稳态计算得到初始应力场和压力场作为瞬态计算的初始条件。这一步很关键煤层原岩应力状态是地应力决定的如果直接从零应力开始算初始渗透率就和现场对不上。稳态计算用全耦合求解器依赖直接求解器PARDISO矩阵规模不大时速度很快。瞬态计算阶段建议改用分离式求解器。把固体力学和渗流场分开迭代每个物理场内部用Newton法物理场之间做几次外部迭代。这样做的原因是THM耦合问题中各个物理场的时间尺度差异很大压力传播快、应力调整慢、温度传播更慢分开求解可以让每个场在自己的时间步上稳定推进不容易出现刚性问题。时间步长用自适应步进初始步长取很小比如1e-5天然后慢慢放大。如果遇到迭代不收敛不要急着调求解器先检查初始条件和材料参数是否合理。COMSOL报错信息里有很多线索比如“在牛顿求解器中找不到收敛”往往是因为初始压力场和边界条件冲突“达西流动不收敛”则常和渗透率数值过小导致矩阵病态有关。调整方法一般是修改初始猜测值、增加阻尼因子或者缩小时间步长。4. 注气方案设计与参数标定4.1 单井注采模型设计现场注气方案千差万别有单井吞吐式、一注一采式、多井井网式建模思路不太一样。我做项目时率先尝试的是单井吞吐模型就是同一口井先注气焖井再开井产气。这种模式最简单也最容易跟数值模拟结果对比适合做机理研究。单井模型里的操作序列我用分段函数实现时间段1注入气体时间段2关闭焖井时间段3开井生产。COMSOL里可以用阶跃函数或者解析函数来定义井的流量或井底压力随时间的变化但这个阶跃点在数值上是一个间断容易引起压力振荡解决方案是在阶跃位置加一个平滑过渡比如用flc2hs函数或者atan函数做斜坡上升过渡时间取1~2小时。对于一注一采模式则需要在几何里加第二口井两口井的距离根据矿区井距设定一般100~300 m。注入井定流量或定压力生产井定井底流压通常取3~5 MPa。在这种模式下注气突破时间是一个重要的关注指标如果生产井的气体组分变化曲线显示注入气过早突破说明注采参数需要调整。4.2 注气参数敏感性分析参数敏感性分析的目的是找到影响增采效果的主导因素。我做这一块时重点考察了注气速率、注气压力、注气温度和焖井时间四个参数。每个参数取三个水平采用控制变量法逐一扫描记录累计产气量和注入气突破比例两个指标。注气速率的影响比较直观速率越高压力传播越快驱替前缘推进越迅速但速率过高会导致近井区域压力超过煤体损伤阈值对渗透率反而有损害。注气压力的选择要参考煤层的破裂压力一般控制在地层破裂压力的80%以下太接近破裂压力容易造成裂隙激活虽然短期渗透率升高但长期看地层稳定性下降水窜风险也增大。注气温度是一个经常被忽略的参数高温气体注入后煤基质受热膨胀如果膨胀受限就会产生较高的热应力对裂隙可能是一个闭合效应实测下来并非温度越高越好存在一个最优范围。焖井时间则需要在“气体充分扩散”和“生产时间被压缩”之间权衡时间太短置换不充分时间太长则产能浪费。扫描完参数后用云图或者散点图画出各参数对累计产气量的响应面通常能直观看出哪个参数最敏感。我项目里的结论是注气压力最敏感其次是焖井时间注气温度敏感度最低。这个结论可以指导现场把有限的资源用在最关键的环节上。4.3 与现场数据的匹配标定数值模型最终要和现场数据对上否则就是自娱自乐。标定的思路是反向拟合先调整渗透率和Langmuir参数让模拟的历史产气量曲线和实际产量曲线尽量重合再检查压力监测数据的拟合程度如果不匹配再回头修改有效应力系数等力学参数。这里有一个经验先标定稳态期的产量再标定瞬态期的产量。稳态期的产量主要受渗透率和边界条件控制瞬态期的产量则更多反映储容系数和吸附参数。逐段匹配可以有效避免多个参数同时调整导致的多解性问题。COMSOL本身不带反演工具我是通过参数化扫描加人工判断来完成的。如果需要系统化的历史拟合可以引入优化模块但大多数情况下人工调参已经能满足工程精度要求。标定完成后模型再外推用来预测各种注气方案的产出效果可靠性才有保障。5. 结果怎么看压力、温度、甲烷产量5.1 压力场与气体前缘演化压力场云图是最直观的结果。注气开始后压力波从注入井向外传播形成一个近似圆形的超压区域到达生产井后压力梯度趋稳流场进入拟稳态。通过COMSOL的动画功能你可以看到压力前锋的推移过程这一点非常有助于向非专业人士解释注气驱替的机制。我特别关注的是压力传播速度和生产井见气时间之间的对应关系。如果压力传播很快而见气时间很短说明优势通道发育注入气体沿高渗裂隙快速突进波及效率差如果压力传播慢且见气时间适中说明注入气体在基质中有效扩散驱替体积大。这种判断直接决定注采方案的调整方向。压力场结果还可以用来评估井间干扰程度。多井模型里如果注气井的压力影响应传到了邻近生产井而该井的产量没有明显上升就要考虑注入气体是否只是“短路”到了生产井而非置换出更多甲烷。5.2 温度变化的工程含义在THM模型里温度场的变化是解读耦合效应的重要窗口。注气温度高于地层温度时注气井附近会出现一个明显的温度升高区域这个区域的大小取决于注入速度和热容匹配。温度升高区域和压力传播区域并不重叠热前锋的推进速度通常远慢于压力前锋这是因为热容储能效应拖延了热的传播反映在注采指标上就是“热突破”比“气突破”滞后得多。温度对甲烷解吸的影响主要体现在两个方面一是热效应直接增强解吸速率二是热应力改变渗透率。项目模拟显示在注气温度高于地层温度15~20°C的条件下局部区域的甲烷解吸量可以提升5%~10%但这种提升会被热应力引起的渗透率降低部分抵消。最终增产幅度取决于哪个因素占上风这个通过模型里的两个输出量解吸量和渗透率变化量的对比就可以定量判断。有一个值得注意的现象是焦耳-汤姆逊效应。当气体在近井带减压膨胀时温度会降低COMSOL的多孔介质传热模型如果包含了压力对温度的影响项可以捕捉到这个降温现象。这个效应对热采方案影响不大但在低温注气场景下可能导致水合物生成需要评估。5.3 甲烷增产效果评估评估增采效果不能只看累计产量还要看产出气组分。设置一个探针在采出井处监测甲烷浓度和注入气浓度的变化据此绘制产出气体组分曲线。理想的曲线特征是生产初期甲烷浓度高98%以上中后期注入气浓度逐渐上升说明注入气开始突破甲烷浓度下降。生产停止的时机应该选在甲烷浓度降到经济下限之前。我用两组模型做对比一组是纯降压开采不含注气一组是注气开采两者采用完全相同的网格和参数只是边界条件不同。最终评价指标是“注气增采量”即注气条件下的累计产甲烷量减去降压条件下的累计产甲烷量。这个差值除以注气量还可以得到单位注气量对应的增采效率这在经济评价中是一个重要指标。另外THM模型相对传统渗流模型的优势要在结果呈现中突出出来。我会做一个对比案例开启动态渗透率更新THM耦合和关闭该更新固定渗透率两种结果往往能发现固定渗透率条件下预测的产量偏高因为现场实际情况中应力封闭效应会抑制渗透率的增长。对比结果可以用来论证THM耦合建模的必要性这个论证在评审或者汇报中非常有力。6. 常见问题与排查技巧实录6.1 不收敛的排查思路做COMSOL遇到不收敛太正常了关键是有一套系统的排查路径。我总结的顺序是先查模型问题再查数值问题最后查代码问题。模型层面的常见问题包括初始条件给得离谱比如初始压力设得和边界条件冲突导致第一步计算压力反转还有材料参数数量级错误比如渗透率填了国际单位制下的值但实际模型单位用了毫达西数值直接差了一千倍。数值层面最常见的是时间步长过大瞬态计算中几个场的时间尺度差异太大固定大时间步会导致压力求解失败。这种情况下把初始时间步长设置到足够小配合自适应步进通常能解决问题。在COMSOL的日志窗口里查看是哪一个物理场的牛顿迭代残差没有收敛非常关键。比如日志提示“固体力学中的牛顿求解器找不到收敛”那问题多半在弹塑性材料模型里常见原因是塑性参数内聚力、内摩擦角设置不合理或者初始应力状态位于屈服面之外。这时候的排查办法是暂时关闭塑性选项先跑通弹性模型再把塑性逐步加回去观察哪一个环节开始发散。6.2 几何导入与网格变形前面提过COMSOL转换CAD拓扑的问题这里具体说说我的处理办法。外部CAD模型导入后如果报“不支持拓扑”第一步是用修复功能自动清理第二步是手动删除细碎特征第三步如果还不行就在CAD软件里重新建模删掉小圆角、小倒角、薄壁只保留主要边界。煤层构造模型的几何一般并不复杂但如果是矿区地表起伏的曲面模型导入出问题的概率会大很多。移动网格是另一个容易出问题的点。在THM模型中如果考虑了煤层的大变形网格节点会随固体位移移动COMSOL需要启用移动网格接口。这个功能在生产井附近容易出现网格畸变尤其在近井区域煤岩发生较大剪切应变时网格质量退化导致求解发散。经验做法是把井周围区域设置较大的积分区域或者采用自适应网格重构。但说实话宏观煤层模型里应变一般不到几个百分点移动网格不是必须的我在大多数情况下会关闭移动网格用固定的欧拉描述来评估。6.3 材料非线性与塑性问题煤层THM模型里最容易出“隐性bug”的地方就是塑性设置。COMSOL的弹塑性材料需要输入屈服函数参数默认的Mises屈服面适合金属材料对煤岩不合适必须改成Drucker-Prager准则或者摩尔-库仑准则。摩尔-库仑准则有角点问题数值稳定性差Drucker-Prager是它在偏平面上的光滑近似工程上更常用。如果计算过程中提示“塑性应变变量在迭代未收敛请检查塑性参数”大概率是硬化模量给的太小或者屈服应力初始值不匹配。处理办法是先给一个较大的硬化模量让计算能推进再慢慢降低到实验值。另外在材料本构中加入了非关联流动法则后如果流动势的膨胀角设得太大体积膨胀会导致计算不稳定膨胀角一般取摩擦角的1/4到1/2不要直接取到很大。关于塑性还有一个容易踩的坑如果你在模型里设置了塑性但是只在局部区域应力超过屈服面而COMSOL默认是全局使用同一套弹塑性本构这会导致非塑性区域的刚度矩阵也变成非对称增加计算负担。更合理的做法是把塑性区域人工限定到井周围或者断层破碎带其余区域用弹性模型这样效率更高也更稳定。6.4 热物性与流场参数的调试心得最后分享一个调试参数的技巧。THM耦合模型的参数太多一旦结果不合理很难判断是哪个参数出问题。我的做法是把耦合逐级拆开先跑单渗流场检查压力分布和产量是否符合预期再加入应力场检查位移和渗透率变化是否合理最后加入温度场看温度分布和热耦合项是否正确。三级递进调试每一步都验证后再进入下一步能省下大量排查时间。做温度场调试时可以把注气温度直接设成地层温度取消热扰动检查结果是否和两级耦合模型一致。如果一致说明传热项和热力耦合项的表达式里有bug如果不一致就顺藤摸瓜找到误差来源。这个“对照组”思想是调试所有耦合仿真的通用方法。至于流体黏度随温度变化COMSOL里要特别注意单位。很多案例里用的黏度公式是基于工程单位制cp的但COMSOL默认的国际单位制是Pa·s换算时要小心。建议把所有自定义函数统一用COMSOL单位制或者在表达式中写清楚换算系数不然一个量级错误就能让整个模型的流动行为面目全非。做了这么多轮COMSOL仿真我个人最深的体会是THM耦合模型的价值不在于“算得有多精确”而在于它能帮我们把现场的模糊认识转变成可量化的判断。注气方案调整之前先在模型里跑几组对比哪些参数敏感、哪些环节有风险心里就有底了。现场实验成本和周期都高能靠仿真提前筛掉一批方案比什么都重要。如果后续想继续往深做可以考虑在这个模型框架里加入离散裂隙网络把宏观裂隙的局部非均质性考虑进去也可以结合地质力学套件做更加精细的应力场模拟还可以把经济评价模块直接嵌到后处理流程里。COMSOL的可扩展性足够支撑这些方向关键是先把THM耦合这个地基打好。
基于COMSOL的煤层气注气增采THM耦合数值仿真全流程解析
做煤层气强化开采的数值仿真最近几年几乎绕不开COMSOL Multiphysics。尤其是注气增采这一块单纯用黑箱模型已经很难满足现场需求了大家都开始往THM耦合这个方向走也就是热-流-力三场耦合。这个项目就是基于COMSOL平台把煤层注气过程中的温度场、渗流场和应力场统一装进一个模型里去评估注气对甲烷解吸、运移和产出的综合影响。如果你正在做煤层气、页岩气或者致密砂岩气的开发方案论证或者你手头有COMSOL但不知道THM耦合怎么做这篇内容可以给你一条相对完整的参考路径。先说清楚这套思路能解决什么问题现场注气强化甲烷开采本质上不是一个单一的渗流问题。注入气体比如氮气或者二氧化碳会改变孔隙压力孔隙压力变化会引起有效应力调整应力调整又会影响裂隙开度和渗透率渗透率变化反过来改变气体流动路径。如果注入的是热气或者注气后伴有温度变化温度场还会通过热膨胀、吸附/解吸热效应进一步扰动这个系统。这就是THM耦合要描述的核心场景。COMSOL在这里的优势很明显它不需要你从底层写有限元代码物理场接口和耦合项都内置了你要做的是把地质模型、边界条件和耦合机制配置对。下面我从项目思路、机理方程、建模实操到参数标定一条条展开讲。1. 为什么在这个时间点做煤层注气THM仿真1.1 煤层气强化开采的现实难题煤层气开采和常规天然气差别很大。甲烷在煤里主要以吸附态存在于微孔表面只有不到两成是游离气。常规排水降压能解吸一部分但越往后解吸速率越慢单井产量衰减很快。很多矿区的井打到中后期日产气连500方都保不住这时候如果还想提升采收率注气是公认的方向之一。原理上不复杂注入二氧化碳或者氮气可以降低甲烷在煤表面的分压通过竞争吸附或者浓度驱替把甲烷置换出来。但问题在于注气不是简单的“往里灌就行”注入气在裂隙网络中的波及范围、压力传播速度、对煤体渗透率的动态影响这些都会决定增采效果。更麻烦的是煤层本身是强应力敏感性的多孔介质。注气抬升孔隙压力后有效应力降低裂隙会张开渗透率随之增大这是有利的一面。但压力升高到一定程度后可能诱发局部塑性变形甚至破坏裂隙又会起变化。温度场参与进来以后情况更复杂煤对温度很敏感升温能促进解吸但同时热应力会改变裂隙开度。这一连串的耦合效应靠室内实验一条条测出来不现实靠现场试错成本太高数值仿真就成了重要的评估手段。1.2 为什么COMSOL适合干这个活COMSOL在THM耦合建模上优势非常明显。它的多物理场耦合机制不是简单的“数据传递”而是在同一个求解框架里把各个物理场的方程联立在一起这对强耦合问题很友好。煤层的THM问题恰恰属于强耦合温度变化影响应力应力变化影响渗透率渗透率变化影响压力分布压力变化又反馈到有效应力和吸附平衡这种多向交互用弱耦合方式计算很容易发散或者失真。再一个就是自定义灵活。煤储层和常规岩层不太一样它需要考虑Langmuir吸附、基质收缩、裂隙压缩等特殊机制。COMSOL里虽然预置了达西流动、多孔介质传热、固体力学这些接口但真正的煤层气模型通常需要在材料属性、源项或者耦合变量里加入自定义公式这一点COMSOL做得很顺手。项目里我把渗透率随应力的变化、气体吸附引起的煤基质变形都做成了自定义表达式整体集成度很高后处理出图也方便给甲方汇报的时候可视化效果比Fortran程序跑出来的结果直观得多。另外还要说一句THM耦合并不仅限于煤层气。多年冻土路基、地热开采、核废料处置这些都是典型的THM场景COMSOL在这些领域已经有大量案例可以借鉴。模型框架搭好之后改材料参数和边界条件换一个工程场景并不难这也是我一开始就选COMSOL而不是其他专业油藏模拟器的核心原因。2. THM耦合机理三场如何互相咬合2.1 热-流-力三场的物理图像理解THM耦合先得有清晰的物理图像。你可以把煤体想象成一块干燥的海绵里面浸满了水但这里的水换成了压缩气体海绵内部的骨架还会随着温度和压力的变化自己变形。温度场负责描述这块“海绵”内部温度的分布和传递热源来自注入气体的温度差异还来自吸附/解吸产生的热量。渗流场负责描述气体和水在裂隙网络里的流动流动的驱动力是压力梯度流动的阻力取决于渗透率和流体黏度。应力场负责描述煤骨架在外力和内部压力作用下的变形与破坏受力状态改变以后裂隙的张开程度会发生变化进而改变渗透特性。这三场不是各自独立跑一遍再拼在一起而是每一步都在互相影响。注气后压力升高孔隙压力变化直接改变有效应力有效应力变化改变裂隙宽度裂隙宽度改变渗透率渗透率改变流量和压力分布。温度升高时煤基质受热膨胀但裂隙边缘的煤块膨胀方向不同会挤压或者拉开裂隙面同时温度升高会加速解吸解吸又带走一部分热量。这种循环依赖关系就是THM耦合的本质。2.2 控制方程与耦合路径在COMSOL里搭THM模型底层方程并不需要你从零推但你必须理解每一项代表什么否则参数填错了很难排查。我这里把关键的方程逻辑梳理一下。渗流场的基础是质量守恒方程COMSOL的达西接口一般写作ρ (∂(ε_p p)/∂t) ∇·(ρ u) Q_m其中u是达西速度ε_p是孔隙度Q_m是质量源项。在煤层气问题里源项不只是注入井的流量还包括甲烷从吸附态变成游离态的产出项这个需要用Langmuir吸附动力学或者平衡解吸模型来描述。动量层面用达西定律u -(k/μ) ∇p这里的渗透率k不是常数是随应力和温度动态变化的这就要和固体力学场联立起来。常用的经验公式是渗透率随有效应力的指数关系k k0 exp(-α(σ_eff - σ_eff0))α是应力敏感性系数煤层的α值通常在0.1~0.6 MPa⁻¹之间具体取决于煤阶和裂隙发育程度。固体力学场用的是孔隙弹性方程COMSOL内置的多孔介质弹性接口已经包含了孔隙压力引起的附加体积力项。对煤层这种裂隙发育的介质还需要在材料本构里加入塑性选项煤在拉压应力超过屈服点后会发生不可恢复的变形这时候需要定义屈服准则和硬化参数。COMSOL里常用Drucker-Prager准则来描述煤岩的塑性行为需要输入的参数包括内聚力、内摩擦角这些可以从煤岩力学实验中获取。温度场是经典的多孔介质传热方程(ρC_p)_eff ∂T/∂t ρ_f C_pf u·∇T ∇·(-k_eff ∇T) Q_T这里(ρC_p)_eff是骨架和流体的等效体积热容k_eff是等效导热系数。Q_T除了外部热源外还要把吸附/解吸热效应加进去煤对甲烷的吸附是一个放热过程解吸吸热这个热效应在注气驱替过程中不能忽略。三场的耦合路径归纳起来有四条孔隙压力变化通过有效应力原理影响应力场应力应变通过渗透率模型反馈到渗流场温度变化通过热膨胀影响应力场吸附/解吸通过源项和质量守恒影响渗流场同时对温度场产生热效应。COMSOL里做耦合时这些交叉项都在各自的系数里关联了所以你不用写耦合代码但要保证各物理场的变量名一致表达式引用正确。2.3 关键物性参数怎么定THM模型对参数极其敏感参数定不准模型再精细也是白搭。我做项目时总结了一套参数获取优先级煤岩力学参数优于热物性参数热物性优于流动参数流动参数优于解吸参数。力学参数里弹性模量、泊松比、内摩擦角、内聚力是必须的。弹性模量煤层一般1~4 GPa比普通砂岩低一个数量级泊松比0.25~0.4内摩擦角20°~40°内聚力0.5~2 MPa。热物性参数方面煤的导热系数0.2~0.5 W/(m·K)比热容1200~1700 J/(kg·K)热膨胀系数10⁻⁵ 1/K量级。流动参数的核心是初始渗透率和孔隙度裂隙煤层的渗透率跨度很大从0.01 mD到几十mD都见过要根据试井数据定。解吸参数用Langmuir体积和Langmuir压力描述这个通过等温吸附实验获取不同煤阶差异很大低阶煤的Langmuir体积能到30 m³/t以上高阶煤可能只有15 m³/t左右。参数太多的时候建议先做一个敏感性分析把影响最大的几个参数标定准其他用经验值就行。这个后面章节细说。3. COMSOL建模全流程从几何到求解器3.1 几何建模与网格划分COMSOL的几何建模思路和CAD软件不一样它更强调“够用就好”。煤层注气模型不需要把每个裂隙都画出来那是离散裂隙网络的事在宏观连续介质模型里裂隙的影响是通过渗透率张量体现的。我建议第一步先建立二维剖面的简化几何比如一个500 m × 50 m的长条形区域中间设置一口注气井两端设置生产井。二维模型跑通以后再扩展成三维这样可以大幅减少前期调试时间。几何建立有两种路径。一种是在COMSOL内直接画矩形、圆形再布尔运算适合规则几何另一种是从外部CAD导入适合实际矿区边界。外部导入时要注意COMSOL的CAD内核和外部软件之间的拓扑兼容性问题我在项目里就遇到过导入的复杂曲面直接被判定为不支持的拓扑最后只能用简化边界替代。这个问题的根源是外部CAD模型里存在微小缝隙或者重叠面COMSOL转换时无法自动修复处理方法是在外部CAD软件里先做几何清理删除小面、合并重合边再以STEP格式导入成功率会高很多。网格划分方面煤层模型需要重点考虑裂隙流动的特征。注气井附近压力梯度和速度梯度大必须加密生产井附近同理。流场变化平缓的区域可以适当放宽。我常用的是三角形网格加边界层在井筒圆周附近设置8~10层边界层网格用来捕捉近井地带的压力漏斗。网格数量控制在10万到50万之间比较合适太多了求解时间成倍增加太少了精度不够。这里要给个建议先跑一个粗网格模型看趋势再加密看关键位置的解有没有明显变化这叫网格无关性验证做完了才能放心用结果说话。3.2 物理场接口与耦合节点物理场接口选择上我的推荐组合是这样固体力学接口负责应力应变计算达西接口负责气体和水的渗流。这里要说一下COMSOL里有两种流动接口可选达西接口和多孔介质流动接口。达西接口更简洁适合单相或者两相简化问题多孔介质流动接口带相传递和毛细压力更贴近真实两相流。我在初始建模时用达西接口单相气体流动先把耦合逻辑跑通再根据需要升级为两相模型。传热用多孔介质传热接口它会自动计算固体和流体的等效热容和导热系数你只需要分别输入煤骨架和流体的热物性参数。如果注气温度和环境温度差异不大很多人会忽略温度场但实际工程中如果采用热气注入或者注气速率很高导致焦耳-汤姆逊效应明显温度场就必须保留。耦合节点的设置是COMSOL里最容易出错的地方。热-力耦合用“多孔介质热弹性”特征这个特征会自动把热应变和孔隙压力引起的应变加入应力方程你只需要激活多孔弹性材料模型并填入热膨胀系数和毕奥系数。流-力耦合主要通过“达西-固体力学耦合”节点实现它会在固体力学方程中加入孔隙压力引起的体积力项同时在达西方程中考虑体积应变引起的储集变化。吸附引起的煤基质变形和渗透率变化需要用变量表达式写入。我在模型里定义了一个渗透率修正因子表达式是基质收缩项和裂隙压缩项的乘积核心思路是当压力降低时基质解吸收缩导致裂隙增大渗透率上升同时有效应力增加导致裂隙压缩两者互相竞争。很多文献里用的Palmer-Mansoori模型就是这个思路你可以在COMSOL里用变量节点把它公式化。3.3 求解器配置与收敛控制求解器选型是COMSOL项目里最需要耐心的一步。THM耦合问题属于高度非线性问题求解器配置不当第一步迭代就发散毫不夸张。我的一般策略先做稳态计算得到初始应力场和压力场作为瞬态计算的初始条件。这一步很关键煤层原岩应力状态是地应力决定的如果直接从零应力开始算初始渗透率就和现场对不上。稳态计算用全耦合求解器依赖直接求解器PARDISO矩阵规模不大时速度很快。瞬态计算阶段建议改用分离式求解器。把固体力学和渗流场分开迭代每个物理场内部用Newton法物理场之间做几次外部迭代。这样做的原因是THM耦合问题中各个物理场的时间尺度差异很大压力传播快、应力调整慢、温度传播更慢分开求解可以让每个场在自己的时间步上稳定推进不容易出现刚性问题。时间步长用自适应步进初始步长取很小比如1e-5天然后慢慢放大。如果遇到迭代不收敛不要急着调求解器先检查初始条件和材料参数是否合理。COMSOL报错信息里有很多线索比如“在牛顿求解器中找不到收敛”往往是因为初始压力场和边界条件冲突“达西流动不收敛”则常和渗透率数值过小导致矩阵病态有关。调整方法一般是修改初始猜测值、增加阻尼因子或者缩小时间步长。4. 注气方案设计与参数标定4.1 单井注采模型设计现场注气方案千差万别有单井吞吐式、一注一采式、多井井网式建模思路不太一样。我做项目时率先尝试的是单井吞吐模型就是同一口井先注气焖井再开井产气。这种模式最简单也最容易跟数值模拟结果对比适合做机理研究。单井模型里的操作序列我用分段函数实现时间段1注入气体时间段2关闭焖井时间段3开井生产。COMSOL里可以用阶跃函数或者解析函数来定义井的流量或井底压力随时间的变化但这个阶跃点在数值上是一个间断容易引起压力振荡解决方案是在阶跃位置加一个平滑过渡比如用flc2hs函数或者atan函数做斜坡上升过渡时间取1~2小时。对于一注一采模式则需要在几何里加第二口井两口井的距离根据矿区井距设定一般100~300 m。注入井定流量或定压力生产井定井底流压通常取3~5 MPa。在这种模式下注气突破时间是一个重要的关注指标如果生产井的气体组分变化曲线显示注入气过早突破说明注采参数需要调整。4.2 注气参数敏感性分析参数敏感性分析的目的是找到影响增采效果的主导因素。我做这一块时重点考察了注气速率、注气压力、注气温度和焖井时间四个参数。每个参数取三个水平采用控制变量法逐一扫描记录累计产气量和注入气突破比例两个指标。注气速率的影响比较直观速率越高压力传播越快驱替前缘推进越迅速但速率过高会导致近井区域压力超过煤体损伤阈值对渗透率反而有损害。注气压力的选择要参考煤层的破裂压力一般控制在地层破裂压力的80%以下太接近破裂压力容易造成裂隙激活虽然短期渗透率升高但长期看地层稳定性下降水窜风险也增大。注气温度是一个经常被忽略的参数高温气体注入后煤基质受热膨胀如果膨胀受限就会产生较高的热应力对裂隙可能是一个闭合效应实测下来并非温度越高越好存在一个最优范围。焖井时间则需要在“气体充分扩散”和“生产时间被压缩”之间权衡时间太短置换不充分时间太长则产能浪费。扫描完参数后用云图或者散点图画出各参数对累计产气量的响应面通常能直观看出哪个参数最敏感。我项目里的结论是注气压力最敏感其次是焖井时间注气温度敏感度最低。这个结论可以指导现场把有限的资源用在最关键的环节上。4.3 与现场数据的匹配标定数值模型最终要和现场数据对上否则就是自娱自乐。标定的思路是反向拟合先调整渗透率和Langmuir参数让模拟的历史产气量曲线和实际产量曲线尽量重合再检查压力监测数据的拟合程度如果不匹配再回头修改有效应力系数等力学参数。这里有一个经验先标定稳态期的产量再标定瞬态期的产量。稳态期的产量主要受渗透率和边界条件控制瞬态期的产量则更多反映储容系数和吸附参数。逐段匹配可以有效避免多个参数同时调整导致的多解性问题。COMSOL本身不带反演工具我是通过参数化扫描加人工判断来完成的。如果需要系统化的历史拟合可以引入优化模块但大多数情况下人工调参已经能满足工程精度要求。标定完成后模型再外推用来预测各种注气方案的产出效果可靠性才有保障。5. 结果怎么看压力、温度、甲烷产量5.1 压力场与气体前缘演化压力场云图是最直观的结果。注气开始后压力波从注入井向外传播形成一个近似圆形的超压区域到达生产井后压力梯度趋稳流场进入拟稳态。通过COMSOL的动画功能你可以看到压力前锋的推移过程这一点非常有助于向非专业人士解释注气驱替的机制。我特别关注的是压力传播速度和生产井见气时间之间的对应关系。如果压力传播很快而见气时间很短说明优势通道发育注入气体沿高渗裂隙快速突进波及效率差如果压力传播慢且见气时间适中说明注入气体在基质中有效扩散驱替体积大。这种判断直接决定注采方案的调整方向。压力场结果还可以用来评估井间干扰程度。多井模型里如果注气井的压力影响应传到了邻近生产井而该井的产量没有明显上升就要考虑注入气体是否只是“短路”到了生产井而非置换出更多甲烷。5.2 温度变化的工程含义在THM模型里温度场的变化是解读耦合效应的重要窗口。注气温度高于地层温度时注气井附近会出现一个明显的温度升高区域这个区域的大小取决于注入速度和热容匹配。温度升高区域和压力传播区域并不重叠热前锋的推进速度通常远慢于压力前锋这是因为热容储能效应拖延了热的传播反映在注采指标上就是“热突破”比“气突破”滞后得多。温度对甲烷解吸的影响主要体现在两个方面一是热效应直接增强解吸速率二是热应力改变渗透率。项目模拟显示在注气温度高于地层温度15~20°C的条件下局部区域的甲烷解吸量可以提升5%~10%但这种提升会被热应力引起的渗透率降低部分抵消。最终增产幅度取决于哪个因素占上风这个通过模型里的两个输出量解吸量和渗透率变化量的对比就可以定量判断。有一个值得注意的现象是焦耳-汤姆逊效应。当气体在近井带减压膨胀时温度会降低COMSOL的多孔介质传热模型如果包含了压力对温度的影响项可以捕捉到这个降温现象。这个效应对热采方案影响不大但在低温注气场景下可能导致水合物生成需要评估。5.3 甲烷增产效果评估评估增采效果不能只看累计产量还要看产出气组分。设置一个探针在采出井处监测甲烷浓度和注入气浓度的变化据此绘制产出气体组分曲线。理想的曲线特征是生产初期甲烷浓度高98%以上中后期注入气浓度逐渐上升说明注入气开始突破甲烷浓度下降。生产停止的时机应该选在甲烷浓度降到经济下限之前。我用两组模型做对比一组是纯降压开采不含注气一组是注气开采两者采用完全相同的网格和参数只是边界条件不同。最终评价指标是“注气增采量”即注气条件下的累计产甲烷量减去降压条件下的累计产甲烷量。这个差值除以注气量还可以得到单位注气量对应的增采效率这在经济评价中是一个重要指标。另外THM模型相对传统渗流模型的优势要在结果呈现中突出出来。我会做一个对比案例开启动态渗透率更新THM耦合和关闭该更新固定渗透率两种结果往往能发现固定渗透率条件下预测的产量偏高因为现场实际情况中应力封闭效应会抑制渗透率的增长。对比结果可以用来论证THM耦合建模的必要性这个论证在评审或者汇报中非常有力。6. 常见问题与排查技巧实录6.1 不收敛的排查思路做COMSOL遇到不收敛太正常了关键是有一套系统的排查路径。我总结的顺序是先查模型问题再查数值问题最后查代码问题。模型层面的常见问题包括初始条件给得离谱比如初始压力设得和边界条件冲突导致第一步计算压力反转还有材料参数数量级错误比如渗透率填了国际单位制下的值但实际模型单位用了毫达西数值直接差了一千倍。数值层面最常见的是时间步长过大瞬态计算中几个场的时间尺度差异太大固定大时间步会导致压力求解失败。这种情况下把初始时间步长设置到足够小配合自适应步进通常能解决问题。在COMSOL的日志窗口里查看是哪一个物理场的牛顿迭代残差没有收敛非常关键。比如日志提示“固体力学中的牛顿求解器找不到收敛”那问题多半在弹塑性材料模型里常见原因是塑性参数内聚力、内摩擦角设置不合理或者初始应力状态位于屈服面之外。这时候的排查办法是暂时关闭塑性选项先跑通弹性模型再把塑性逐步加回去观察哪一个环节开始发散。6.2 几何导入与网格变形前面提过COMSOL转换CAD拓扑的问题这里具体说说我的处理办法。外部CAD模型导入后如果报“不支持拓扑”第一步是用修复功能自动清理第二步是手动删除细碎特征第三步如果还不行就在CAD软件里重新建模删掉小圆角、小倒角、薄壁只保留主要边界。煤层构造模型的几何一般并不复杂但如果是矿区地表起伏的曲面模型导入出问题的概率会大很多。移动网格是另一个容易出问题的点。在THM模型中如果考虑了煤层的大变形网格节点会随固体位移移动COMSOL需要启用移动网格接口。这个功能在生产井附近容易出现网格畸变尤其在近井区域煤岩发生较大剪切应变时网格质量退化导致求解发散。经验做法是把井周围区域设置较大的积分区域或者采用自适应网格重构。但说实话宏观煤层模型里应变一般不到几个百分点移动网格不是必须的我在大多数情况下会关闭移动网格用固定的欧拉描述来评估。6.3 材料非线性与塑性问题煤层THM模型里最容易出“隐性bug”的地方就是塑性设置。COMSOL的弹塑性材料需要输入屈服函数参数默认的Mises屈服面适合金属材料对煤岩不合适必须改成Drucker-Prager准则或者摩尔-库仑准则。摩尔-库仑准则有角点问题数值稳定性差Drucker-Prager是它在偏平面上的光滑近似工程上更常用。如果计算过程中提示“塑性应变变量在迭代未收敛请检查塑性参数”大概率是硬化模量给的太小或者屈服应力初始值不匹配。处理办法是先给一个较大的硬化模量让计算能推进再慢慢降低到实验值。另外在材料本构中加入了非关联流动法则后如果流动势的膨胀角设得太大体积膨胀会导致计算不稳定膨胀角一般取摩擦角的1/4到1/2不要直接取到很大。关于塑性还有一个容易踩的坑如果你在模型里设置了塑性但是只在局部区域应力超过屈服面而COMSOL默认是全局使用同一套弹塑性本构这会导致非塑性区域的刚度矩阵也变成非对称增加计算负担。更合理的做法是把塑性区域人工限定到井周围或者断层破碎带其余区域用弹性模型这样效率更高也更稳定。6.4 热物性与流场参数的调试心得最后分享一个调试参数的技巧。THM耦合模型的参数太多一旦结果不合理很难判断是哪个参数出问题。我的做法是把耦合逐级拆开先跑单渗流场检查压力分布和产量是否符合预期再加入应力场检查位移和渗透率变化是否合理最后加入温度场看温度分布和热耦合项是否正确。三级递进调试每一步都验证后再进入下一步能省下大量排查时间。做温度场调试时可以把注气温度直接设成地层温度取消热扰动检查结果是否和两级耦合模型一致。如果一致说明传热项和热力耦合项的表达式里有bug如果不一致就顺藤摸瓜找到误差来源。这个“对照组”思想是调试所有耦合仿真的通用方法。至于流体黏度随温度变化COMSOL里要特别注意单位。很多案例里用的黏度公式是基于工程单位制cp的但COMSOL默认的国际单位制是Pa·s换算时要小心。建议把所有自定义函数统一用COMSOL单位制或者在表达式中写清楚换算系数不然一个量级错误就能让整个模型的流动行为面目全非。做了这么多轮COMSOL仿真我个人最深的体会是THM耦合模型的价值不在于“算得有多精确”而在于它能帮我们把现场的模糊认识转变成可量化的判断。注气方案调整之前先在模型里跑几组对比哪些参数敏感、哪些环节有风险心里就有底了。现场实验成本和周期都高能靠仿真提前筛掉一批方案比什么都重要。如果后续想继续往深做可以考虑在这个模型框架里加入离散裂隙网络把宏观裂隙的局部非均质性考虑进去也可以结合地质力学套件做更加精细的应力场模拟还可以把经济评价模块直接嵌到后处理流程里。COMSOL的可扩展性足够支撑这些方向关键是先把THM耦合这个地基打好。
