我在做华北某区块煤层气排采试验井的模拟时第一版模型怎么都匹配不上试井数据后来发现问题不在井筒而在少算了温度。正是那次经历让我把煤层气开采里的热流固三场耦合从“听说过”变成了“一条总能对上的逻辑链”。这篇内容就围绕煤层气数值模拟中热流固三场耦合的机理、建模路线、参数取值和踩坑记录展开适合正在做煤储层模拟、从事非常规天然气开发或者刚开始接触多物理场耦合方程的工程师。1. 三场之间到底在“耦合”什么单场和两场模型为什么不够1.1 吸附态为主煤层气与常规气藏的本质差异常规天然气主要是游离气储集在孔隙和裂缝里开采时压降驱动气体流向井筒流固耦合的影响当然存在但很多时候可以简化处理。煤层气不一样绝大多数甲烷以吸附态附着在煤基质微孔内表面一般占总含气量的80%到95%。想采出这部分气必须先把压力降下来让吸附态甲烷解吸成游离气再经过基质扩散进入裂隙系统最后通过裂隙渗流到井筒。这是煤层气开发最核心的逻辑链排水—降压—解吸—扩散—渗流。链条每个环节都不是孤立的。降压会改变煤储层的有效应力有效应力变化又让裂隙张闭进而改变渗透率压力降低导致甲烷解吸煤基质收缩裂隙重新张开解吸本身是吸热过程温度场随即变化温度变化又会反过来影响吸附平衡和煤岩力学性质。单看任何一个场都能说通但放在一起就说不清了。所以“热流固三场耦合”不是学术上炫技的提法而是煤层气开采中实际存在的物理过程。温度场看似温和——储层温度变化往往只有几度——但它参与反馈的方式是持续的长期排采下影响会累积到不可忽略的程度。1.2 排采过程中的三重连锁反应我在实际建模时会把热流固三场之间的相互作用拆成三条主线去理解。第一条是流固耦合排水降压后储层孔隙压力降低上覆岩层压力基本不变煤骨架承受的有效应力增加。煤是裂隙性多孔介质裂隙对有效应力极其敏感有效应力一增大裂隙就被压得更紧渗透率下降。这在排采初期尤其明显也是很多低渗煤层气井产量迟迟起不来的原因之一。第二条是基质收缩降压到临界解吸压力以下后甲烷开始解吸煤基质随之收缩裂隙宽度反而增大渗透率回升。基质收缩的强弱取决于煤阶和含气量高煤阶煤的收缩应变可以达到1%左右低煤阶相对小一些。这个效应与有效应力效应是竞争的一个让渗透率下降一个让渗透率上升最终曲线呈现“U形”或“先降后升”的特征。第三条是温度反馈甲烷解吸吸热井筒周围温度下降温度降低让煤基质发生热收缩裂隙再变化同时低温会增强煤对甲烷的吸附能力等于给解吸过程“踩了刹车”。如果储层温度较高吸附能力本来就弱解吸更容易一旦局部温度降下来同一压力下的平衡吸附量反而上升产气速率会变慢。这三条线交织在一起就是标题里“变形记”的含义煤储层在排采过程中不是静态的骨架而是随压力、温度、吸附状态不断变形演化的活体。1.3 “变形记”到底指什么从力学角度看煤储层的变形来源至少有三个有效应力变化引起的压缩变形、温度变化引起的热应变、吸附/解吸引起的基质膨胀或收缩应变。三者同时发生方向还不一致所以不能简单线性叠加。我经常打一个比方把煤储层想象成一块泡发过的海绵。用力挤它水被挤出来海绵变薄——这是有效应力把海绵放进热环境里它膨胀——这是热应变把海绵里原本吸附的某种“溶胀剂”抽走它收缩——这是吸附应变。煤层气排采就是一边挤压、一边降温、一边抽走吸附质这三个过程在同一时间轴上进行海绵最终变成什么样取决于三者谁占上风。数值模拟要做的就是把这三者的竞争关系用方程表达清楚并在时间步进中保持它们的一致性。这也是后面所有参数设置和求解策略的出发点。2. 控制方程里的耦合通道有效应力、孔隙度与渗透率的勾稽关系2.1 应力场方程中多出来的三个力煤储层中的应力平衡方程形式上比普通弹性力学方程多出三个耦合项。第一个是孔隙压力项。煤体承受的总应力一部分由煤骨架承担一部分由孔隙流体承担这就是有效应力原理。用修正的Terzaghi有效应力表达有效应力等于总应力减去Biot系数与孔隙压力的乘积。Biot系数反映骨架和孔隙的刚度关系煤的Biot系数一般接近1但也有不少学者取0.8到0.9这个参数不能随便设要结合声波测试结果标定。第二个是热应力项。温度变化产生热应变在约束条件下转化为热应力。煤的线膨胀系数大约在1×10⁻⁵到3×10⁻⁵ /K量级不算大但在三场耦合里足够影响裂隙开度。排采井周围温度下降几度产生的热应变虽然小于基质收缩应变但在中后期也不容忽视。第三个是吸附应变项。这是煤储层特有的。吸附甲烷时煤基质膨胀解吸时收缩常用Langmuir型公式来描述吸附应变与压力相关达到极限值ε_l。不同煤阶的ε_l差异很大低阶煤可能只有0.2%到0.5%中高阶煤可以达到1%以上。这个参数直接决定基质收缩效应能多大程度抵消有效应力效应。把这三个项写进本构关系后应力场方程才真正与渗流场、温度场形成了双向耦合。只保留孔隙压力项的模型本质上就是常规流固耦合只有三项全包括才能算热流固三场耦合。2.2 渗透率动态模型PM与SD的取舍流动方程的核心是渗透率如何随生产动态变化。目前业界用得最多的两类经典模型是PM模型和SD模型。PM模型Palmer-Mansoori从裂隙孔隙度出发认为渗透率变化受有效应力压缩和基质收缩共同控制。裂隙孔隙度满足φ_f / φ_f0 1 (c_m / φ_f0)(p - p0) (ε_l / φ_f0)(K_m / M - 1)[ (b p) / (1 b p) - (b p0) / (1 b p0) ]其中c_m是基质压缩系数K_m是基质体积模量M是约束轴向模量b是Langmuir型压力常数。渗透率假设与孔隙度满足立方关系k / k0 (φ_f / φ_f0)^3这个模型的优点是物理意义清晰能直观看到有效应力项和基质收缩项之间的竞争。缺点是它对裂隙几何做了很多简化参数多每项都需要实验标定。SD模型Shi-Durucan走的是另一条路它认为渗透率与有效应力满足指数关系k / k0 exp[-3 c_f (σ_eff - σ_eff0)]σ_eff随孔隙压力和基质收缩变化所以同样能模拟出渗透率先降后升。这个模型更贴合裂隙应力敏感强的储层但参数c_f裂隙压缩系数和后期的渗透率回升幅度绑定在大需要谨慎。我的做法是粗算阶段用PM模型因为它对细观机理交代得清楚历史拟合阶段换SD模型或者把两套模型都跑一遍用现场产量数据倒逼参数。没有哪个模型是普适的关键看你手里的数据能标定到哪一层。2.3 温度场不是陪跑热源项与解吸吸热的量级估算温度场方程相对直接多孔介质的等效热容乘温度变化率加上流体的热对流项减去有效导热系数控制的扩散项等于源汇项。麻烦的是源汇项的种类和量级。煤层气温度场的主要源汇有三个。一是解吸吸热。甲烷从煤表面解吸要吸收热量这个热量大约在15到25 kJ/mol之间我一般取18 kJ/mol做初步估算。18 kJ/mol换算成单位质量大约是1.1 MJ/kg甲烷。假设1吨煤含气量25 m³标准状态下约17.9 kg甲烷完全解吸需要吸收热量约20 MJ。煤的比热容大约在1到1.5 kJ/(kg·K)单算煤体绝热温降可以超过10℃。当然实际排采中有水流和地层热补给温度不会降这么夸张但井周区域降温3到5℃是常见的模拟结果。这就足以影响吸附平衡和热应力了。二是焦耳-汤姆逊效应。甲烷在降压膨胀过程中温度会降低在地面管网里这种现象非常明显储层内同样存在。甲烷的焦汤系数在储层温压条件下大约0.4到0.5 K/MPa一口井从原始压力降到井底流压压降好几兆帕气体膨胀降温贡献1到2℃并不奇怪。三是粘性耗散。流体在裂隙中流动时摩擦生热这个量级通常很小在常规排采模拟中我一般忽略只有在注热增产或压裂液返排这类注入工况下才把它放回去。温度场的源项处理直接决定了温度模拟结果是否可信。我在第一次做三场耦合时把解吸热写成了常数源结果温度场出现明显漂移后来调整成与解吸速率联动的形式整个能量守恒才闭合。这是数值实现里最容易忽略、也最容易出错的地方。3. 从地质参数到可收敛模型热流固三场耦合的建模流程3.1 几何建模与网格井筒附近加密的尺度控制我做概念模型时通常建一个500m×500m范围、储层厚6m的单井模型模拟区域太大网格数量会失控太小又体现不了远井区的补给效应。网格划分时井筒附近必须加密因为压降漏斗、温度降和应力集中都集中在这个区域。网格加密的经验是井筒径向从内向外按等比数列扩张相邻网格尺寸比控制在1.1到1.2之间。这样做的好处是既保证井壁附近几十厘米到几米范围内的分辨率又不让远井区的网格数量爆炸。如果网格尺寸从井壁的0.5m直接跳到煤储层边界的50m相邻网格体积差太大有限元插值误差会成倍放大求解器很容易振荡。单元类型上也有讲究。煤岩泊松比往往在0.3以上有些情况下接近0.4接近不可压缩。使用一阶线性单元做力学计算可能会出现体积自锁应力结果失真。我一般对力学场采用二阶单元或者采用混合单元方案流体压力用低阶插值位移用高阶插值这种搭配能显著改善收敛性。另外如果模型中包含人工裂缝裂缝周围的网格还要再加密一个级别。但第一版模型我不建议急着把裂缝细节放进去先把热流固三条反馈跑通确认边界和参数都没问题再加裂缝否则出了数值问题很难定位是裂缝网格引起的还是耦合算法引起的。3.2 耦合求解策略全耦合与顺序耦合的组合用法三场耦合的数值实现本质上是如何组织三个物理场的求解顺序。主要有两条路线。全耦合monolithic是把位移、压力、温度、饱和度等所有未知量放进同一个非线性方程组一起Newton迭代求解。优点是物理一致性最好每一步得到的都是三场同时平衡的解缺点是矩阵规模大、内存消耗高而且非线性强时收敛窗口很窄网格稍微差一点就容易发散。顺序耦合staggered/operator splitting是先解某个场再用结果去更新其他场循环迭代直到收敛。优点是灵活、省内存可以复用成熟的单场求解器缺点是如果时间步长取得很大场与场之间的信息交换存在滞后可能引入数值误差。我在实际项目中用组合方案流固两场采用全耦合因为它们之间的反馈最快、最强温度场用顺序耦合接入在每一个时间步内先算流固场再算温度场更新吸附参数和温度应变后回头迭代。这样既保证了最关键的流固反馈不脱节又避免了把温度场强行拉进全局矩阵带来的收敛压力。时间步长的控制也直接影响成功与否。排采初期压降传播快解吸尚未大规模启动时间步需要压得很小我一般从0.01天开始进入平稳排采期后逐步放大到0.5天甚至1天。自适应时间步长功能要打开当某一步Newton迭代不收敛时自动减半重算这是最实用的兜底策略。3.3 边界条件与初始场不能凭感觉给的开采参数边界条件看起来简单实际是错误高发区。力学边界方面模型顶部施加上覆岩层自重对应的垂向应力侧向根据侧压系数取0.8到1.2倍垂向应力底部固定垂向位移。如果模型直接从地表建起顶部直接加自由面边界如果是箱体截断模型顶部应力值必须按实际埋深计算不能随便给一个数。我见过不少模型上覆应力取少了导致有效应力偏低模拟出来的渗透率一直上去产气量虚高。渗流边界方面外边界采用定压条件数值取原始储层压力井底采用排采制度控制最常用的是定井底流压或定产水量。排采制度要写成时间函数比如前90天以每5天降0.2 MPa的速度从原始压力降到4 MPa而不是一步到位。实际工程中为了控制煤粉产出和裂隙闭合降压速率比这更慢模拟时也应当对应。温度边界方面外边界取地层原始温度井底不能简单定温更合理的是对流换热边界让流体与井壁之间按换热系数交换热量。忽略这一点井周温度会被定得过低或过高三场耦合的效果就被歪曲了。初始场设定最关键的是初始有效应力平衡。给定了孔隙压力场和重力载荷后初始位移必须由静力平衡方程解出来不能直接置零。否则模型一开始就处于不平衡状态会“假变形”这个变形量会叠加在排采引起的变形上导致渗透率演化曲线失真。初始含水饱和度也不能随便给要与储层深度、相渗曲线端点匹配否则气水两相渗流在头几个时间步内剧烈振荡。4. 360天排采模拟渗透率、井周温度与产气曲线的联动4.1 典型储层参数与排采制度设计为了让内容具体我放一组典型华北中高阶煤储层参数这套参数来自我做过的一个概念模型用来研究三场耦合对产能预测的影响参数取值备注埋深800 m中深部煤层储层压力8 MPa近似静水压力储层温度35℃地温梯度约2.5℃/100m煤岩弹性模量3.2 GPa煤岩常见范围2~5 GPa泊松比0.35中高阶煤典型值初始裂隙孔隙度0.04裂隙孔隙度远小于基质孔隙度初始渗透率0.5 mD低渗储层典型值Langmuir压力2.5 MPa等温吸附实验拟合Langmuir体积25 m³/t对应含气量中高最大基质收缩应变0.01中高阶煤偏大煤体导热系数0.3 W/(m·K)干煤偏低饱和水后更高解吸热18 kJ/mol换算约1.1 MJ/kg甲烷排采制度按阶梯降压设计第1到30天井底流压从8 MPa降到6 MPa第30到90天降到4 MPa第90到360天再缓慢降到2 MPa。全程控制压降速率避免有效应力骤增对裂隙造成不可逆伤害。4.2 渗透率“先降后升”的经典演化路径这套参数下的模拟结果很有代表性。前30天井周有效应力迅速增加裂隙被压缩渗透率从0.5 mD下降最低会到0.32 mD左右降幅接近40%。这个阶段对应井筒附近压降漏斗的快速扩张煤基质收缩还没有充分启动有效应力效应完全占主导。进入第90天后储层压力进一步降到临界解吸压力以下甲烷解吸范围扩大基质收缩逐渐成为渗透率的主导因素。渗透率开始回升到第200天左右恢复到初始值第360天可以达到0.75 mD左右。整条曲线呈现明显的“U形”或“先降后升”形态。这个形态对现场的意义很大。如果只做流固耦合、不考虑基质收缩渗透率会随有效应力持续下降模拟出来的产气量会明显偏低但如果只强调基质收缩前期渗透率下降被弱化又会高估早期产量。三场耦合里的“热”在这里的作用是调节基质收缩的启动速度——温度降低抑制解吸基质收缩来得比纯两场模型更慢、更平缓。所以三场耦合预测的渗透率回升点会比两场模型晚一点最终稳态渗透率也低一些。4.3 井周温度降对产气峰值的“刹车”作用温度模拟结果显示井筒周围50m范围内降温最明显模拟360天时井壁附近温度比初始值低3.5℃左右。这个降温由解吸吸热主导辅以气体膨胀的焦汤效应。离井筒越远温度降越小100m处只有零点几度因为地层热传导和流体热对流在持续补充热量。温度降的反馈体现在吸附平衡上。甲烷吸附是放热过程温度降低会增强煤对甲烷的吸附能力。我用温度修正后的Langmuir参数重算发现每降低1℃相同压力下的平衡吸附量大约增加2%到3%。井周温度降低3.5℃相当于在局部给了解吸过程一个明显的反向阻力。对比试验更有说服力同一套参数下忽略温度场仅做流固耦合时第250天左右达到产气峰值峰值日产气约3200 m³/d考虑完整热流固三场耦合后峰值时间推迟到第270天附近峰值降到2800 m³/d左右降幅约12%。长期累计产气量也有约8%的差距。这个结果表明温度效应虽然“温和”但对产能预测的修正量完全在工程决策需要考虑的范围之内。5. 收敛性、单位制与实际校准不跑一遍不会知道的坑5.1 数值不收敛时的四个调优旋钮三场耦合模型出数值问题是常态不收敛的原因千奇百怪但我排下来发现80%集中在四个方面。第一个是网格尺寸跳跃。相邻单元尺寸比超过3局部插值误差就会显性化表现为压力或温度场出现“棋盘式振荡”。解决办法很简单网格过渡区域多设两层过渡单元让尺寸连续变化。第二个是初始时间步长过大。排采刚开始时井底压力突变非线性方程组需要很小的步长才能找到平衡。把初始步长压到0.01天前几个步长跑稳定后再逐步放大能解决大部分“第一步就发散”的问题。第三个是物性参数剧烈跳变。比如渗透率模型里基质收缩参数取太大渗透率在解吸启动瞬间从一个量级跳到另一个量级方程刚度剧增。我的处理是把基质收缩参数的敏感性分开测试先不激活温度场单独看渗透率曲线是否连续如果不连续就是参数范围不合理。第四个是强非线性材料的单元退化。当煤岩接近不可压缩时线性单元可能出现体积自锁表现为位移场严重振荡。升级为二阶单元或者混合单元后这个问题基本消失。如果求解器不支持混合单元也可以把泊松比稍微调低一档做试探但这样会牺牲物理精度只适合排查。5.2 单位制与吸附/热参数的一致性检查单位制混乱是三场耦合里最隐蔽的坑。有限元软件本身没有单位概念它只认数值大小。同一个方程里压力用Pa还是MPa、渗透率用m²还是mD、时间用秒还是天如果混用矩阵条件数会差好几个数量级结果自然一塌糊涂。我习惯统一采用SI基本单位长度m、时间s、质量kg、温度K、压力Pa、渗透率m²。煤层气行业常用的mD和MPa需要换算1 mD约等于9.869×10⁻¹⁶ m²工程上粗略记为10⁻¹⁵ m²1 MPa是10⁶ Pa。吸附参数也要同步转换Langmuir压力2.5 MPa在SI里是2.5×10⁶ PaLangmuir体积25 m³/t换成kg/kg需要乘以甲烷密度0.717 kg/m³再除以1000。最容易出问题的是解吸热和吸附量的单位配套。如果解吸热用J/mol而吸附量用m³/t两者不在一个计量体系能量源项就会差一个数量级。我在第一个三场模型里就栽过这个跟头温度场整体偏低十几度折腾了两天才发现是解吸热单位换算少了分子量这一步。后来我把所有涉及能量的物理量全部转换成“J/kg甲烷”吸附量也换算成“kg甲烷/kg煤”模型才稳定下来。另外气体性质不能当常数处理。甲烷压缩因子随着压力下降明显偏离1尤其在低压阶段如果不做修正产气量会被高估。我倾向于在模型中加入气体偏差因子与压力的关联式或者在每个时间步用闪蒸计算更新气体密度这样产气曲线更贴近现场。5.3 用试井和产气数据反推模型参数参数再精确模拟结果也必须经过现场数据校准这是整个流程的最后一道关。我的校准顺序是先用试井压降数据反演初始渗透率因为排采初期的压力传播速度对渗透率最敏感再用压力恢复曲线约束裂隙压缩系数这段数据能反映渗透率随有效应力变化的敏感程度然后用产气曲线的峰值和峰期位置匹配Langmuir参数和解吸热最后才轮到热物性参数。实际操作中我很少追求所有参数一步到位。先固定力学参数手调渗透率渗透率对上了再放开吸附参数吸附参数跑通了最后才动温度相关参数。这种逐层开放的策略能避免多参数同时调整导致的“参数补偿效应”——两个错误参数互相掩盖问题表面拟合很好但模型一点预测能力都没有。我在做历史拟合时以现场月产气量为目标要求模拟结果与实测误差控制在15%以内同时井底压力的下降趋势也要基本吻合。如果产气对但压力不对说明渗透率参数可能错了如果压力对但产气不对多半是吸附参数或解吸热的设置有问题。这种交叉验证比单纯对一条曲线要可靠得多。最后说几个跑完这个项目后的个人习惯。我拿到一个新模型一定会先做“温度场关闭/开启”的对比试验如果温度对结果影响小于2%说明储层温度系统本身就不敏感可以简化处理一旦超过5%热流固三场耦合就没得商量必须做。还有个小心得是每次调整参数后都把渗透率、温度和产气量三条曲线存一份快照跑完几十组工况后回看这些快照很容易发现参数之间的弱关联规律。三场耦合的计算量确实比两场耦合大不少但这些前期工作做扎实了后期反而更省时间。