煤层气THM耦合模拟：渗透率孔隙度模型与现场应用要点

煤层气项目的现场数据，是最会骗人的。我前几年在某区块做开发方案时，实验室测出来的岩心渗透率是0.5毫达西，按这个数建模型，预测的单井产气量能到每天两千方。结果打下去一压裂一排采，实际稳定产气量连零头都不到，日产一百多方。当时整个团队都懵了，反复查压裂施工记录、排采制度，始终找不到原因。后来把问题归到储层物性动态变化上，才意识到我们一直用“静态渗透率”去算一个“动态过程”的账——煤层气开采本身就是一场压力、应力、温度交互作用的拉锯战。这也是我今天想聊透这个话题的原因：煤层气运移的THM耦合模型，以及与之配套的渗透率孔隙度模型，到底是怎么一步步搭起来、又怎么用到实际项目里的。

这套东西不是什么前沿黑科技，就是一套把温度场、渗流场、应力场放在同一个方程框架里联立求解的研究思路。它的核心价值在于，能把“排采降压→有效应力增加→裂隙压实闭合”“解吸→基质收缩→裂隙重新张开”这一对矛盾机制同时反映出来，再叠加温度变化对吸附解吸和岩石骨架的扰动。适合同行里正在做数值模拟、做开发动态分析、或者正被现场渗透率异常折磨的工程师和科研人员参考。如果你只是需要做一个静态储量计算，这篇内容可能对你帮助不大；但只要你的模型需要对历史产气曲线做拟合、要做开发指标预测，渗透率不能当常数，那你就是这篇文章的目标读者。

1. 从现场数据对不上说起：为什么煤层气运移必须看三场耦合

1.1 单场模型失灵的两个案例

先聊两个我自己遇到或者身边同事遇到过的真实案例。第一个是典型的“降压反而不增产”：有一口排采井，套压降到很低了，产气量却没有像教科书说的那样稳定上升，反而在某个阶段下跌。单看流动方程的话，压力降低意味着解吸驱动力增大，产气量应该涨。但事实是，大压降带来了明显的孔压衰减，有效应力急剧上升，裂隙被压实，渗透率掉了好几个数量级，产气通道被“掐住”了，这个问题靠单纯的渗流模型无解。

第二个案例是注入实验和数值结果严重错位。做CO₂驱替煤层气先导试验时，注入井端的渗透率应该是增加的——因为CO₂吸附导致煤基质膨胀，理论上会挤压裂隙，反而降低渗透率？这里逻辑很绕，实际数据出来以后，预测和实测对不上。单纯用等温吸附曲线去算储层膨胀，忽略了注入过程带来的温度扰动和应力扰动，结果偏差自然大。这两个案例指向同一个结论：煤储层不是刚性骨架，压力和温度一变，储层本身也在变，必须把流动和变形放到同一框架里看。

1.2 THM耦合的三个场到底在耦合什么

THM是Thermo-Hydro-Mechanical的缩写，对应温度场、渗流场、力学场。要理解耦合，不能把三个场分开看，要关注它们之间的交叉作用。

渗流场（H）和力学场（M）的耦合是最容易理解的。孔隙压力变化会改变有效应力，有效应力公式很简单：σ′ = σ − αP。煤体受压变形，裂隙开度改变，反过来渗透率就变了。而渗透率一变，渗流阻力随之变化，压力分布又不一样，这是一个闭环。

温度场（T）的加入是很多人一上来容易忽略的。温度升高，煤基质受热膨胀，裂隙受压有闭合趋势；与此同时，温度升高会加速气体解吸，解吸导致的基质收缩又会抵消掉一部分热膨胀。煤体的热膨胀系数和吸附应变这两个参数往往被单独拿出来测，实际问题里它们同时在起作用，数学上就体现为温度场对力学场和渗流场都贡献了一个源项。

三个场联立以后，方程组会变得非常硬：流动方程系数的变化剧烈，力学方程涉及塑性或弹塑性行为，热传导方程又有长时间尺度特征。这也是为什么很多商业模拟器单独跑一个模块都很顺，一旦耦合就各种发散、不收敛。因为时间步长、网格尺寸、非线性迭代策略都得跟着换。

1.3 温度场不是配角：深部开采与注热增产的现实需求

早些年做煤层气模拟，很多人不把温度当回事，温度项直接设成恒温。这在浅部开采、温度变化不大时确实够用。但这两年深部煤炭资源开发、注热增产的研究越来越多，温度效应再也绕不过去。

深部煤层温度一般超过40摄氏度，有些区块地温梯度高，储层温度能到60度以上。温度影响吸附等温线的形态，Langmuir压力参数和温度强相关，温度升高10度，吸附量可以从2.5降到1.8左右，这个变化对产气潜力是数量级影响。注热增产的思路更直接：通过热流体注入把煤储层温度抬起来，利用热解吸效应抢出一部分吸附气。但温度升高同时带来热应力、煤体膨胀、有效应力的再分布，处理不好就会把裂隙胀裂或者压实，增产变成减产。这种工程场景，不考虑THM耦合就是空中楼阁。

2. 渗透率孔隙度模型拆解：从P&M公式到工程修正

2.1 P&M模型的推导逻辑：两个机制在抢渗透率

做煤层气数值模拟的人，对P&M模型这个名字应该不陌生。P&M模型的核心逻辑说起来其实很朴素：把煤体想象成一块被均匀裂隙切割的方糖，裂隙开度决定渗透率高低。开采过程中，裂隙开度的变化由两项竞争：

一项是有效应力压实项。孔隙压力降了，有效应力升了，煤骨架被压缩，裂隙开度变小。这一项永远站在渗透率的对立面。另一项是基质收缩项。压力降到解吸临界点以下，吸附气开始脱离煤基质表面，基质本身发生收缩，裂隙被“拉开”，渗透率增加。这一项又是抬升渗透率的。

两个机制互相拉扯，最终渗透率走向取决于谁占上风。P&M模型把它公式化成渗透率比：

k/k₀ = (1 + 3/b₀ · (Δb_stress + Δb_shrinkage + Δb_thermal))³

这里的b₀是初始裂隙开度，三个Δb分别代表应力变化引起的开度变化、基质收缩引起的开度变化、温度变化引起的开度变化。开度的总变化量占初始开度的比例，然后取三次方，放大成渗透率的变化。因为渗透率对裂隙宽度的依赖是立方级别的，裂隙开度即使只变化百分之几，渗透率也可能翻倍或者腰斩。这也就解释了为什么现场渗透率常常“测不准，变太快”。

2.2 参数到底怎么取：Langmuir应变、裂隙压缩系数、初始孔隙度

公式好写，参数难取。P&M模型真正难用的地方在参数标定。我吃过亏的地方，重点说三个参数。

第一个是裂隙压缩系数C_f。它描述有效应力每增加一个单位，裂隙孔隙度衰减多少，单位通常是MPa⁻¹。这个参数在不同煤阶里的差异极大：高阶煤刚性更强，C_f可能只有0.05左右；低阶煤松散，C_f可以到0.3以上。更麻烦的是，室内实验用的小煤样往往裂隙发育不完全，压出来的C_f和原位储层出入很大。我见过有的文献直接用0.1作为通用值，这在实际项目中风险极高。你不做一次三轴加载渗透率实验，完全不知道自己手里的煤样C_f到底落在哪个区间。

第二个是Langmuir应变参数ε_max。这个参数的物理意义是：吸附量达到饱和时，煤基质的极限体积应变。它一般是正值，表示膨胀；解吸时则是负的收缩量。ε_max的取值受煤阶、含水率、温度共同影响，一般在0.002到0.01之间波动。注意，这是应变，不是百分比。换算成裂隙开度变化时，还要乘上一个基质块体尺寸和裂隙间距的比例系数，很多人在这里搞混。

第三个是初始裂隙孔隙度φ₀。P&M模型的早期形式会把C_f/φ₀作为一个组合参数出现，因为裂隙孔隙度本身就是开度除以间距的比值。问题是，裂隙孔隙度在地质上是极难测准的，测井解释给的是一个范围，岩心分析又代表不了整个储层。实操中我会花大力气先做一次单孔压降试井，反演初始渗透率和初始裂隙孔隙度，而不是直接拿室内孔隙度值替代。

2.3 模型对比：P&M与Shi-Durucan谁更适合你的场景

除了P&M，业内常用的还有Shi-Durucan模型，两者常常被拿来对比。它们最大的区别在于：P&M把基质收缩直接写成裂隙开度的几何变化，Shi-Durucan则把基质收缩折算成有效应力的一部分，然后套指数衰减关系。

Shi-Durucan的核心表达式是：

k/k₀ = exp(−3C_f · Δσ_eff)

其中Δσ_eff不单是外应力变化，还扣掉解吸引起的“内应力”变化。这个模型的好处在于形式简洁，工程实现容易；缺点是“内应力”的物理图像比较牵强，参数拟合痕迹重。

我自己的选择标准是这样的：做机理研究、需要解释渗透率反弹现象时，优先用P&M；做工程预测、现场数据有限时，Shi-Durucan的指数形式说不定更稳。两个模型都是对真实行为的近似，别迷信任何一方。

2.4 氧化损伤-微裂隙扩展：渗透率模型的新拓展方向

传统的P&M和Shi-Durucan把渗透率归因于“应力+解吸+温度”三条线，但有些人忽略了另一条重要路径：氧化损伤导致的微孔微裂隙扩展。

这个方向是从煤岩流体化学反应延伸出来的。煤体在氧气或氧化性流体作用下，基质表面发生氧化反应，生成CO₂和有机小分子，煤骨架局部被“蚀刻”，微孔连通性增强，微裂隙沿力学薄弱面扩展。表现出来，就是渗透率在不改变外载条件的情况下，随时间缓慢上升。这个机制在煤层自燃风险区、空气驱替、原位氧化增产等场景里非常明显。

我在自己的模型里给这条路径加了一个工程近似项：

k_aug = k_base · (1 + γ·D_ox)

D_ox是氧化损伤变量，取0到1，γ是裂隙连通修正因子。这个写法不是从某一篇文献直接抄的，而是综合了几个氧化实验的孔隙度变化趋势做的经验拟合。它最大的作用是提醒你：渗透率演化不只是纯力学过程，化学反应也可以在长时标里改动物性。做长期数值模拟的时候，完全不考虑这一项，历史拟合到后期就会系统性偏移。

3. 实操环节：这套模型怎么落到具体项目里

3.1 实验设计与参数获取

要拿到可靠的模型参数，关键在试验方案设计，而不只是设备好坏。我建议按下面这个流程走一遍：

岩样准备阶段，钻取标准圆柱样，直径25毫米或50毫米均可，但一定要注明取心方向。煤岩层理发育较强，水平渗透率和垂直渗透率相差可以达到五到十倍，建模之前先搞清楚自己模拟的是哪个方向的流动。

加载阶段分三步走：第一步，固定围压，逐步改变孔压，测不同有效应力下的渗透率；第二步，固定有效应力，逐步升高温度到储层温度附近，测温度-渗透率响应；第三步，把气压降低到解吸临界以下，全程监测应变和渗透率变化，捕捉基质收缩的贡献。第三步是最容易出问题的——很多实验只测到“渗透率下降段”，没有把“反弹段”测出来，后续拟合P&M模型时，收缩参数完全无法约束。

数据处理时，把每一阶段的应变数据、渗透率数据、压力数据对齐到同一时间轴上。计算有效应力时要注意，孔压采用的是平均孔压还是端面压力，直接决定了C_f的数值，务必保持一致口径。

3.2 数值实现路径与软件选型

模型选型和参数标定一样重要。市面常见的实现方式有三类：

第一类是COMSOL Multiphysics这类多物理场平台。固体力学模块、达西流动模块、固体传热模块直接组合，耦合方式清晰。适合机理研究和方案对比。缺点是非线性耦合时计算量偏大，一个中等规模的三维模型跑一个工况可能就需要十几个小时。

第二类是FLAC3D这类岩土力学软件。它力学求解强，内置的裂隙单元和渗流模块能处理煤岩破坏和裂隙扩展，做压裂、开采扰动研究很合适。缺点是温度模块相对基础，如果你想做精细的热-流-固耦合，需要二次开发。

第三类是TOUGH系列。它在热-流耦合方面非常成熟，尤其适合CO₂注入、热注入这类涉及多组分多相流的场景。但它和力学变形的耦合是弱耦合，通常需要借助外部接口把应力场计算结果导回来。有人开玩笑说TOUGH是“算流动的利器，算力学的半残”，话糙理不糙。

实际项目中我倾向于两步走：用COMSOL或FLAC3D做机理阶段的小规模模拟，确定主控因素；再借助商业油藏模拟器做全区块历史拟合和指标预测。商业模拟器虽然没有单独的THM模块，但通常提供了渗透率随压力动态变化的表格式输入，你把P&M模型计算的k(p, T)结果导入进去，效果已经很接近真正的耦合模拟，而且稳定性和效率都好得多。

3.3 历史拟合与参数校准

历史拟合不是“调参数调到对”，而是“找到哪些参数值得调”。我在项目里有一个习惯：先做敏感性分析，把模型里10个参数挨个上下浮动50%，看产气量曲线变化幅度。变化大的参数列为重点校准目标，变化小的参数直接用经验值固定，不再折腾。

以我那个区块为例，敏感性分析结果显示，前三敏感的依次是Langmuir压力P_L、裂隙压缩系数C_f、初始渗透率k₀。这三个参数直接决定了产气量峰值和衰减节奏。随后我结合室内实验结果初设参数，再用现场压降试井和早期产气数据反演修正。拟合时不要一上来就调全曲线，优先锁产水阶段和早期产气段，因为这两个阶段对应的是应力敏感与解吸启动，正是渗透率模型信息量最丰富的时期。等模型在这些阶段拟合满意了，再去匹配后期平台期和中低产期，基本就水到渠成。

4. 现场踩坑实录与排查速查表

4.1 三个踩过高频坑

第一坑：渗透率算成负值。有些工况压降特别大，基质收缩项加上热应力项把裂隙开度算到了负值，渗透率变成负的，模型直接崩掉。解决办法是给裂隙开度设下限，比如最小开度不小于初始开度的百分之五。并检查力学计算里是否加了接触应力约束，物理上裂隙压实到一定程度后刚度是急剧上升的，这部分需要一个非线性修正。

第二坑：参数不敏感性误判。我一开始用商业模拟器做参数拟合，发现改变C_f对产气量没有明显影响，以为这个参数不重要。后来才发现，模拟器的渗透率压力表默认用指数形式，且压力区间没覆盖到解吸段，参数压根没激活。这个坑提醒我，先确认模型软件里渗透率模型的形式和你的公式假设一致，再做敏感性分析，不然就是在一个假命题里折腾。

第三坑：室内实验数据和现场数据对接不上。实验室的煤样通常已经经历了卸压和干燥，裂隙和原始地应力状态与地下完全不同，测出来的渗透率往往偏高。我吃过的亏是拿着室内值直接当初始渗透率代入模型，导致预测产量整体偏高。修正办法是用现场试井解释值做基准，把室内实验数据主要用于拟合各参数的相对变化趋势，而不是绝对数值。

4.2 问题排查速查表

排查时还有一个经验：最先怀疑的一定是输入参数异常，而不是方程写错。模型报错，十次有九次是给了一个物理上不合理的参数组合。比如初始孔隙度取0.05但C_f取0.3，裂隙力学上根本不可能同时成立。先做一轮数学一致性检查，把参数范围框到合理区间，再动方程结构。

关于参数范围，我给你一个实操参考区间：C_f通常在0.02到0.3 MPa⁻¹之间，ε_max在0.002到0.01之间，初始裂隙孔隙度φ₀在0.005到0.05之间。超出这个区间的取值，不要直接使用，先回过头检查实验或者试井解释结果。

数值模拟方面还有一个稳定性的土办法。三场耦合的非线性很强，直接上生产工况很容易发散。我先用单场模型预热几百步，把压力场和应力场算出一个相对平滑的初值，再开启耦合迭代。这跟开车前先热车一个道理，看似多花了时间，实际上省掉后面反复试错的时间。迭代时加一个压力场的阻尼因子，让每步压力变化不超过上一步的20%，也能有效压制振荡。

关于现场数据的问题，再补充一条：排采初期压力下降快的时候，渗透率可能是先降后升的，现场上往往表现为产气量出现一段“酝酿期”之后才上升。这个并不是井筒积液或压裂效果差，很可能就是煤储层的P&M渗透率响应。做生产分析时要把这段“酝酿期”和真正的低效生产区分开，否则容易误判压裂效果，过早调整工作制度。我的做法是先把渗透率模型跑出来，看理论上的酝酿期大概多长，再对照现场数据判断是否有其他因素叠加。

最后分享一个实际项目中的应用思路。我们当时需要在方案阶段对不同井距方案做比选，直接跑全耦合模型工作量大到爆炸。后来我用P&M模型先把“渗透率随压力路径的变化关系”做成一个查表文件，嵌入到相对简单的商业化生产模型中，计算效率和全局耦合模型差不多，但速度却快了两个数量级。这种“机理模型离线标定，工程模型在线调用”的思路，在时间紧、又不愿意牺牲准确度的方案论证阶段，非常实用。我自己做完这个项目后，最大的体会是：THM耦合模型和渗透率模型不应该被当作两套独立的东西，它们是一层递进的关系——前者提供宏观框架，后者提供局部物性演化的内核。把这个关系理顺了，不管是做机理研究还是做工程预测，心里都会踏实很多。