企业官网建设流程全解析

做煤层气开发或者瓦斯治理工程的人，这两年应该对“THM”这个缩写不陌生。THM 三个字母分别对应热（Thermal）、流体（Hydraulic）、力学（Mechanical），放在煤层注气增强甲烷开采的场景里，就是研究注气过程中温度场、渗流场和应力应变场之间怎么互相影响。以前很多人建模只考虑气水两相渗流，把煤体当成刚性骨架，早期做工程估算还行，但真要在注气方案上做决策，误差往往大得离谱。COMSOL Multiphysics 这类多物理场仿真软件，恰好提供了一个把热力流三套方程装进同一个模型中迭代求解的环境，所以这几年在非常规天然气、煤层气方向的应用频率明显上升。

这篇文章想聊的，是我在 COMSOL 里搭建煤层注气 THM 耦合模型、模拟增强甲烷开采过程的一些流程和心得。理论推导可以写得很复杂，但真正上手操作时，边界条件怎么设、渗透率用哪个演化模型、求解器怎么调，每一步都可能把你卡在原地。文章适合刚开始接触多场耦合的研究生，也适合在矿区做瓦斯治理和增产方案的技术人员，跟着思路重新搭一套自己的模型，完全没有问题。

1. 煤层注气和 THM 耦合的项目背景

1.1 渗流、温度和应力到底怎么互相纠缠

煤层注气的物理过程，一句话可以概括成：通过井筒向煤层注入二氧化碳、氮气或者混合烟道气，把原本吸附在煤基质微孔表面的甲烷置换出来，再通过产气井采出。听起来是纯渗流问题，但煤层的特殊性在于，煤是一种双重孔隙介质，基质里有大量微孔，天然割理和裂隙则构成宏观渗流通道。甲烷主要以吸附态存在，想要让它解吸，就必须打破原有的吸附平衡。

打破吸附平衡的手段不止一种。降低孔隙压力可以让甲烷从吸附态转为游离态，这对应常规压降开采；注入二氧化碳则利用煤对 CO2 更强的吸附亲和力，把甲烷从吸附位上“挤”下来，这属于置换驱替。但问题也随之而来：煤基质吸附 CO2 后会发生明显的基质膨胀，基质膨胀会挤压割理裂隙，导致渗透率下降；反过来，甲烷解吸又会让基质收缩，孔隙通道扩大，渗透率上升。膨胀和收缩这对矛盾，直接受温度和流体压力的影响，而流体压力和温度又随注气时间和注气参数变化，于是“热-流-力”三个物理场就彻底绑在了一起。

用单物理场模型处理这个链条，最典型的错误是只算压力场变化，不管应力场，得出的渗透率往往是恒定值，或者只按压力做简单修正。可实测数据经常显示，同样的注气压力，不同温度、不同围压条件下的增产效果差异很大，根源就在煤体的力学响应。所以只要问题涉及注气方案优化、注入压力优选、煤层渗透率动态演化，就必须把 THM 耦合纳入考虑。

1.2 COMSOL 做这类问题有哪些优势

COMSOL 的优势不在于哪一物理场算得最精细，而在于“耦合”这件事做得足够顺手。它有独立的“达西接口”“固体力学接口”“流体传热接口”，每个接口背后都是成熟的偏微分方程框架，然后通过多物理场节点的预设耦合项把三者串起来，比如把体积应变和孔隙度变化引入达西方程，把流体压力和温度变化作为体积力引入力学方程。整个过程可以逐步搭建、逐项验证，比从零写有限元程序效率高很多。

另外，COMSOL 的变量定义很灵活，大多数商业软件里参数只能填固定系数，这里可以直接写渗透率随孔隙度、吸附应变变化的表达式，甚至可以用全局 ODE 或者偏微分方程接口再补一套自定义本构关系。做科研需要反复修改模型假设时，这种开放性非常关键。

当然，COMSOL 也不是没有短板。它对超大规模三维网格的并行计算效率，和专门的油藏数值模拟软件相比并不占优，而且煤层气藏中的裂缝网络非常复杂，纯连续介质假设容易摸不到真实储层的脾气。所以比较合理的定位是：用 COMSOL 做注气机理研究和参数敏感性分析，把规律先吃透，再用更专业的油藏软件做矿区级规模的历史拟合和产能预测。

2. 物理模型、控制方程与关键参数设置

2.1 热流力三场怎么装入同一套方程

建立 THM 模型之前，首先得把煤体简化为连续介质，用孔隙弹性理论描述煤的力学行为。这样可以写出第一个核心方程，固体变形平衡方程。在 COMSOL 固体力学接口中，默认解的是动量和应力平衡，关键在于把有效应力表达式写对：σ' = σ - αp，其中 α 为 Biot 有效应力系数，p 为孔隙压力。温度引入膨胀项，所以总应力本构里还要加一项与温差成正比的热应变。

第二个核心方程，渗流方程。气体和水在割理裂隙中流动满足达西定律，COMSOL 的达西接口可以处理多孔介质中单相或多相流动。做注气模拟时以气相为主，也可以考虑气水两相，但两相设置会让方程非线性更强，收敛难度明显上升。如果前期以机理研究为主，先做单相气体流动，重点把渗透率演化效应做进去。

第三个核心方程，传热方程。煤炭是一种导热性能一般的多孔介质，但注气气体温度、地层初始温度、吸附/解吸过程的热效应，都会引起温度场变化。吸附热一定要考虑，因为甲烷解吸是吸热过程，二氧化碳吸附是放热过程，这两者产生的热量变化直接影响近井地带温度，进而改变煤体热应变和吸附平衡，最终反馈到渗透率上。

三个方程通过三个桥梁变量耦合：一是孔隙度 φ，变化后直接改变渗透率和储容系数；二是有效应力，它由骨架应力和孔隙压力组合，决定煤体的压缩或膨胀状态；三是吸附膨胀应变，通常表示成 Langmuir 形式的附加应变项，与压力、温度均有关系。

2.2 吸附模型与动态渗透率演化模型

吸附模型我通常选择 Langmuir 等温吸附的扩展形式。甲烷组分和二氧化碳组分在煤表面存在竞争吸附，二元混合气体的吸附量可以表达为：

V_i = (V_i,max · b_i · p_i)/(1 + b_1 p_1 + b_2 p_2)

其中 V_i 是组分 i 的吸附量，V_i,max 是 Langmuir 体积，b_i 是 Langmuir 常数，p_i 是组分分压。温度对吸附的影响不可忽略，通常把 Langmuir 常数写成温度和吸附热相关的指数形式，比如 b_i = b0 · exp(ΔH_i/(R·T))。这一步很关键，假如忽略温度修正，温度场再怎么算也都是“陪着跑”，对渗透率演化起不到实际作用。

渗透率演化模型是 THM 模型里最容易引起争议的地方。Palmer-Mansoori 模型、Shi-Durucan 模型等各有各的假设，选择时要根据煤体的力学参数和现场应力条件来定。Palmer-Mansoori 模型把渗透率表示为割理孔隙度、有效应力、基质收缩应变共同作用的结果，表达式相对常见，也比较容易在 COMSOL 中作为变量写入。Shi-Durucan 模型则更多考虑水平渗透率与水平应力的关系，适用于水平割理发育的矿区。

以 Palmer-Mansoori 模型为例，渗透率比值可以直接写成与孔隙度比值相关的三次方形式，而孔隙度变化由压力变化和吸附应变共同控制：

k/k0 = (φ/φ0)^3

φ/φ0 = 1 + (c_f/φ0)((p-p0)/K) - (ε_l/φ0)(K/K_M - 1)(p/p+p_ε) ...

这类公式看起来参数多，但在 COMSOL 里就是一组变量定义，写一次之后所有物理接口都能调用。重要的是把各符号代表的物理量搞清楚，别在单位上翻车。

2.3 参数取值参考与边界条件概览

不同矿区煤层的力学性质差异非常大，没有通用的万能参数。我这里列一组我国典型中阶烟煤的常见数量级参数，方便刚开始接触的同学做参考：

表格：

物理量 | 常用单位 | 常见取值范围 | 备注 孔隙度 | % | 2%～8% | 割理孔隙度，基质孔隙占比较小 Langmuir 体积 | m3/t | 15～30 | 换算体积按标准状态气体 Langmuir 压力 | MPa | 1～3 | 吸附难度相关 Biot 系数 | 无量纲 | 0.6～1.0 | 煤体裂隙发育程度高时取大值 弹性模量 | GPa | 2～6 | 割理发育煤体偏低 泊松比 | 无量纲 | 0.25～0.4 | 煤体塑性较强 最大吸附膨胀应变 | 无量纲 | 0.01～0.03 | 取决于气体组分，CO2 吸附时偏大 渗透率 | mD | 0.1～10 | 原煤渗透率现场离散很大

初始条件下，地层孔隙压力等于原始地应力对应的有效应力状态，温度取煤层原位温度；注气井位置上，给定注入压力和注气温度；产气井给定井底流压。边界上最外侧一般取远场边界，压力保持原始地层压力不变，位移边界为约束法向位移。

注意，参数不是越精确越好。在机理研究阶段，很多参数只求数量级正确，重点是观察各物理场之间的相对变化趋势。到具体矿区应用时，弹性模量、渗透率、Langmuir 参数这类敏感参数，最好通过实验室岩心试验、现场试井数据来标定。

3. COMSOL 建模实操：从模块装配到求解收敛

3.1 模块选择和多物理场装配

COMSOL 的操作界面里需要装配的核心节点包括：流体传热、达西定律、固体力学。建模时建议直接从三维组件里建立二维轴对称模型，也就是在平面内表示注气井周围的径向剖面。二维轴对称的好处非常明显，几何建模快、网格量少、求解速度快，而且注气井问题的物理场本身就具有轴对称性，先用二维跑通全流程，再扩展到三维检验边界效应，是性价比最高的路径。

实际操作中，我会先在模型向导中选择三维空间维度，然后添加物理场接口，分别勾选固体力学、流体传热、达西定律，在后续设置里把研究类型选为瞬态。注意，COMSOL 的二维轴对称几何必须用 x 轴代表半径方向，y 轴代表深度或高度方向，这一点后期导出结果时特别容易搞混。

多物理场耦合节点建议这样添加：第一个耦合节点把孔隙压力和温度变化作为体载荷加载到固体力学中，公式表现为增加的应力项；第二个耦合节点把固体力学计算得到的体积应变传递到达西接口中，用于更新孔隙度和渗透率；第三个耦合节点考虑吸附热作为传热方程的热源项。添加顺序并不重要，重要的是把每个耦合变量都检查一遍，确保单位和方向正确。

3.2 注气井和产气井的边界条件处理

在 COMSOL 中，井筒通常被简化成边界条件而不是实体几何。比如注气井位置可以画成一条边界线，然后设置成给定压力边界，也就是 Dirichlet 边界，值取注入压力；远场边界设置成压力恒定，模拟无限大地层。产气井如果和注气井不在同一模型中，可以单独设置另一组边界条件，也可以通过多口井模型共同求解，几何相对复杂一点，但思路相同。

这里有一个容易被忽略的问题：注入气体温度影响吸附过程和煤体温度场，但在近井地带，注气井入口气体温度和地层温度之间可能存在几十度的温差，初始时刻传热现象非常剧烈。如果初始条件设置成整个地层温度均匀，注气井温度直接跳到目标值，计算早期往往出现剧烈的数值振荡。解决办法是给注气压力或注气温度加一个缓慢升高的瞬态过渡函数，比如前 0.5 天从初始值缓慢过渡到目标值，既能避免数值不收敛，也更贴近现场开井时逐渐提压提注的实际操作。

力学边界条件方面，模拟注气引起的地面抬升或沉降时，表面设为自由边界，底部和远端设为法向位移约束。如果要模拟储层处于一定围压条件下，则需要在外部边界上加地应力面载荷。

3.3 网格划分与求解器设置的经验

THM 耦合模型的网格精度问题，我踩过不少次坑。最初我把网格画得很密，觉得算得更准，结果瞬态计算一跑就是十几个小时，到后期还频繁出现不收敛。后来经验是：近井地带网格加密，因为压力梯度和温度梯度大，远井区域网格渐变放粗，固定在 2～5 米，整体网格量控制在几万量级以内。如果做完网格无关性检验，再增加网格数量对结果影响低于 1%，就没必要继续加网格了。

求解器设置上，COMSOL 的瞬态研究默认会用自适应时间步，但 THM 耦合方程刚性强，建议手动改成 BDF 算法，最大阶数为 2，初始步长设得非常小，比如 1e-3 秒级别，然后逐渐放大。同时要仔细设置每个物理场的相对容差，不要把容差放松，一般相对容差取 1e-4 或更小，否则渗透率演化容易出现台阶式突变。求解过程中如果一直报“找不到一致的初始值”，常见原因是方程中出现了除以孔隙度的表达式，而初始孔隙度给成了 0 或者极小值。这种问题排查起来很费时间，但本质上就是检查每个变量的初始值是否落在方程定义域内。

4. 结果分析：透过云图和数据曲线看增强甲烷开采规律

4.1 温度、压力、应力场的动态演化关系

结果输出之后，第一步不是直接看甲烷产量，而是把三类云图放在一起对照。

注气初期，近井地带压力最先升高，压力波以井筒为中心向外扩散，温度则随着注入流体进入逐渐改变。如果注入的是低温二氧化碳，温度下降的区域和压力升高的区域高度重合。温度下降会引起煤体收缩，有一点增加渗透率的正面作用，但此时压力升高引起的有效应力增大，又会让割理被压缩，渗透率下降。两方面作用互相抵消，最后反映到实际产量上，可能增产也可能减产，完全取决于参数配比。建模时如果不耦合温度场，这类竞争机制就完全看不见。

应力场的变化更值得注意。煤体在注气初期受到压力波挤压，近井地带的体应变表现为压缩状态；当注入气体扩散到足够范围，吸附作用逐渐占主导，基质体积开始膨胀，但膨胀受到外边界约束后，会产生环形拉应力区。环形应力区的存在，一方面会改变裂隙开度，另一方面在工程上还可能诱发微震或者裂隙扩展。COMSOL 里查看应力分布云图时，要特别留意应力集中带的移动规律。

4.2 渗透率演化曲线怎么看

渗透率演化曲线是 THM 模型最核心的输出。常规做法是把渗透率对初始渗透率的比值画成随时间变化的曲线，同时把孔隙度演化曲线作为辅助。结果通常不是单调变化，而是呈现先下降后上升或者先上升后下降的趋势，关键看吸附膨胀和有效应力哪个先占上风。

如果模型起点是原始煤层，注气压力较高时，初期压力效应占主导，渗透率快速下降；随着时间推移，吸附膨胀效应积累，渗透率可能继续下降，也可能因甲烷解吸引起的基质收缩而出现回升。现场试验也观察到类似的非单调变化规律，这正好验证了 THM 耦合的必要性。很多工程人员拿到渗透率曲线后，只关注终值，其实波动过程中的极值点、拐点反而隐藏着更多关于储层参数的信息。

4.3 从模拟结果反推注气方案

模拟的最终目的是指导工程决策。我常用的分析手段是参数敏感性扫描：把注入压力、注入温度、注气持续时间作为变量，设置 5～10 组数据，对比累计甲烷产量和注气量比值。以注入温度为例，当注入气体温度较高时，甲烷解吸速率提高，但同时对煤体的热应力和吸附平衡的影响很复杂；低温注气能增强吸附气体的封存能力，但可能降低甲烷解吸速率。通过扫描曲线找到拐点，就能给现场提出合理的注入温度范围，而不是摸着石头过河。

另外，模型结果还能回答“该不该提高注气压力”这类现实问题。注气压力升高能加快驱替，增强甲烷产量，但带来的副作用是有效应力增大、渗透率降低，而且可能增加煤岩破裂风险。在 COMSOL 里把不同井底压力条件下的应力场调出来，检查任何点是否超过煤体拉压强度，就能大致划定安全注气压力区间。

5. 实际建模中的常见问题与实用排查技巧

5.1 如何对付不收敛

THM 模型不收敛的原因，排在第一位的是初始条件设置不当，尤其是孔隙压力与应力场的初始条件没有执行“初始平衡”。COMSOL 中求解稳态应力场时，如果初始位移全部为零，而初始孔隙压力分布与实际地应力不匹配，模型就会报“求解器未收敛”。我的做法是先单独跑一遍只有地应力和孔隙压力的稳态求解，得到初始应力场，然后把结果作为瞬态模型的初始值，这样能避免很多早期数值震荡。

排在第二位的是变量表达式出现断点。比如 Darcy 速度方程里用了隙宽、渗透率的三次方，当渗透率趋近零时，方程会变得极刚性。检查全局变量定义里是否有除法、幂运算中存在潜在零值，把表达式改写成数值稳定的形式，比如添加一个极小正常数作为保护项。

排在第三位的是时间步长太激进。COMSOL 自适应步长在强非线性问题中很容易自动放大步长，造成下一步无法收敛。这时候不要硬等，直接暂停，把研究设置里的最大时间步长主动限制在合理范围内，再从头接着算。

5.2 单位、网格与几何简化的细节

单位是 COMSOL 初学者最容易忽略的隐形陷阱。几何可以用米，也可以用毫米，但所有物理量单位必须统一到同一套国际制单位。特别是渗透率，现场习惯用 mD，但 COMSOL 达西接口默认单位是 m2，1 mD 约等于 9.869233e-16 m2，经常出现换算错误导致结果偏大几个数量级。建议在做任何后处理分析前，先检查单位设置，并手工计算一遍基础工况，确保量级没有跑偏。

网格无关性验证一定要做。曾经我有一组模型，网格粗细变化一倍，渗透率演化曲线居然差了将近 30%，问题出在局部流动速度梯度太大，粗网格没有捕捉到近井压力骤降。如果网格数量增加后曲线基本重合，才可以放心继续正算。另外，如果有注气井和产气井同时存在，井边界附近网格至少要加密到局部尺寸的 1/10，否则井底流压相关结果偏差很大。

5.3 模拟结果为什么和现场数据对不上

很多工程人员拿到模拟结果后第一反应是和现场实测数据对比，发现偏差不小，就认为模型不可靠。其实这不一定说明模型有问题，而要先检查现场数据本身的误差。比如渗透率现场值往往通过注入压降测试获得，测试期间煤体处于非平衡状态，和模型里稳定初始条件的背景并不同频。

更常见的是模型简化过头了。单孔、单井模型忽略了煤层气的天然裂缝网络、局部断层、水相分布，这些因素在地质条件复杂的矿区是决定性的。若现场数据来自多煤组互层，建模时就要考虑层间窜层；若现场测试段内含水率高，优先考虑气水两相模型，而不是单纯的气相达西模型。我的建议是，先用模拟结果和现场数据之间的差值反推哪些关键参数需要修正，把这个过程理解为参数标定而不是模型错误。经过几轮标定后，模型预测能力和可靠性会显著提升。

6. 建模思路的进一步扩展

6.1 从单孔模型到双重介质模型

前面提到的连续介质模型把煤体看作等效多孔介质，这在宏观尺度上是合理的，但煤的双重介质特性决定了基质微孔和割理裂隙的交换阻力会影响注气前缘的推进速度。如果模拟结果和实验数据之间的偏差始终存在，可以考虑把模型升级为双重介质双重渗透率模型。COMSOL 可以用两个达西接口分别对应基质量和裂隙量，再用交换项连接两者。这个升级会让方程数量翻倍，收敛难度加大，但能捕捉到注气过程中“裂隙先响应、基质后响应”的动态特征，在分析注气焖井时间、压力恢复曲线时价值很大。

6.2 引入煤层气抽采与注气协同的多井场景

实际工程中很少是单井注气、单井采气，更常见的是井网部署。注气井阵列和采气井阵列之间会形成压力波叠加、气体突破和窜流。COMSOL 里可以按实际井网布置建立平面二维模型，把每口井都设为边界条件，研究不同井距、不同注采比下的增产效果。多井场景下的计算量比单井高出不少，但规律也更贴近现场，能够顺带回答“井距缩到几米以后发生气体突破”这类工程关键问题。

6.3 力学本构的拓展方向

煤作为一种强度较低、塑性特征明显的沉积岩，弹性本构只是第一步。当注气压力较高、接近地层破裂压力时，塑性变形不可忽略，此时可以考虑在固体力学接口中引入塑性本构或损伤本构。COMSOL 提供了多种材料非线性选项，也可以基于用户自定义的失效准则进行二次开发。前期先把 THM 弹性模型做扎实，再逐步增加非线性效应，这样即使后边模型复杂了，也清楚每一步变化的来源是什么。

一点个人体会

做煤层注气 THM 耦合模拟这几年，我最大的感受是：这个方向真正的难点不在软件操作，而在物理过程的理解和简化。COMSOL 把方程接口做得很友好，但模型里放了哪些耦合项、哪些被忽略，直接决定结果靠不靠谱。每加一个物理过程，都应该问自己一句：这个耦合在这个尺度下真的关键吗？它对结果的改变是真实机制，还是只是模型自我安慰式的复杂化？

如果让我给刚开始接触这个方向的人一个建议，那就是不要一上来就追求大而全的三维双孔双渗模型，先把二维轴对称的 THM 模型跑通，把渗透率演化、吸附热、温度场这些机制逐项验证，再逐步扩展。模型不是越复杂越高级，能够解释现场现象并指导下一步实验和工程决策，才算真正用到位。