1. 煤储层渗透率为何"先降后升":一切建模工作的事实起点
从常规天然气储层转到煤层气项目的人,第一个认知冲击往往不是解吸机理,而是渗透率的变化规律。常规储层在衰竭开采中渗透率通常缓慢单调下降(孔压降低导致骨架压实),但煤储层完全不同。我跟踪过两口井的长期试井数据,初始渗透率大约5 mD,排采一年半内掉到了0.8 mD左右,之后又开始慢慢回升,两年后回到接近初始值的水平。这种非单调演化让很多刚开始做煤层气数值模拟的人一头雾水,甚至以为试井解释出了问题。
其实这个现象背后是两个互相竞争的物理机制在"拔河":
- 有效应力压缩机制:排采过程中孔隙压力持续降低,上覆岩层压力不变,作用在煤体骨架上的有效应力相应增大,导致天然裂隙(割理)开度被压缩,渗透率下降。在排采初期,这个效应主导。
- 基质收缩机制:甲烷从煤基质微孔表面解吸后,基质体会产生收缩变形,相当于给割理系统"松绑",裂隙开度反而增大,渗透率回升。当解吸量积累到一定程度后,这个效应开始压过应力压缩效应。
这两个机制的变化速率、响应幅度决定了现场渗透率曲线的具体形态。更麻烦的是,如果煤层埋深大(超过1000米)、地温高,温度场的变化还会通过改变吸附平衡和煤体热应变进一步影响渗透率演化。这正是必须把热(Thermal)、水(Hydraulic)、力(Mechanical)三个物理场耦合起来建模的根本原因——如果你只用单场或双场模型,根本描述不了完整的排采过程。
这篇文章我打算把煤层气运移THM模型的骨架、渗透率-孔隙度模型的选型逻辑、数值实现中的关键细节,以及我自己踩过的坑完整梳理一遍。适合正在做储层数值模拟的研究生、刚接触煤层气开发的工程师,以及想从"用软件"进阶到"理解软件"的同行参考。
2. THM三场耦合的物理骨架:温度、渗流、应力如何互相咬合
THM耦合听起来高端,本质上就是描述三个物理场之间的相互作用。煤储层和常规储层最大的区别在于:固体骨架(煤基质)对流体压力和温度的变化高度敏感,而这种变形又反过来改变流体的流动通道。三个场谁都不能单独拎出来说"我已经考虑得很完备了"。
2.1 力学场:有效应力是连接流体压力和煤体变形的那根轴
力学场遵循经典的连续介质力学框架。对煤体这种裂隙发育的介质,通常采用双重孔隙或等效连续介质模型来描述。核心变量是应力张量σ、应变张量ε、位移u。本构关系一般采用线弹性或弹塑性模型,但对于工程尺度的排采模拟,线弹性假设(配合适当的模量)往往已经够用,非要上非线性本构反而容易在参数标定上陷入泥潭。
这里面的关键枢纽是Biot有效应力原理:
σ' = σ - αP
其中P是裂隙系统中的孔隙压力,α是Biot系数(裂隙煤体通常接近1)。孔隙压力下降,有效应力升高,煤体被压缩。但需要注意的是,煤层在排采过程中并不是"自由压缩"——上覆岩层是恒定的应力边界,而平面方向又受到周围煤体的约束,所以实际变形状态接近单轴应变条件(水平位移为零,垂向自由沉降)。这个约束条件直接决定了理论模型中应力应变的换算关系,也决定了不同渗透率模型之间那点微妙的差别。
力学场还耦合了温度效应。温度升高时煤基质产生热膨胀,等效为附加的应力/应变源项;同时吸附态甲烷解吸时基质收缩,也产生一个"负膨胀"效应。这两项在深层煤层气开发中分量不小,尤其是热采方案设计时。
2.2 渗流场:从基质解吸到裂隙流动的三步接力
渗流场描述的是气体和水在煤储层中的运移过程。煤储层是典型的双重孔隙介质:基质体内部有大量微孔(孔径从纳米到微米级),比表面积巨大,甲烷绝大部分以吸附态储存在这一部分;基质之间发育着面割理和端割理,构成天然裂隙网络,是气水流动的主要通道。
运移过程可以概括为三步接力:
- 解吸:压力降到临界解吸压力以下后,吸附态甲烷从微孔表面脱离,遵循Langmuir等温吸附方程。V = V_L·P/(P_L+P),V_L是Langmuir体积(最大吸附量),P_L是Langmuir压力(吸附量达到一半时的压力)。
- 扩散:解吸后的甲烷在基质微孔中通过浓度差驱动向裂隙系统扩散,服从Fick定律。基质块尺寸越小、扩散系数越大,这一步越快。
- 渗流:甲烷进入裂隙系统后,与裂隙中的水一起按达西定律向井筒流动。这里涉及气水两相相对渗透率、毛管压力、Klinkenberg效应(低压下气体滑脱)等一系列问题。
我把这三个环节在模型里统一处理成"基质-裂隙质量交换项",基质内用扩散传质,裂隙内用两相达西流。这样做的好处是物理过程清晰,每个参数都有明确的实验测定方法;代价是计算量增加,尤其当基质块划分较细时。
2.3 温度场:吸附热和热应变的双重角色
温度场在THM模型中往往被当成"配角",但深层煤层气和热采方案研究时它是实打实的主角。
首先是温度对吸附平衡的影响。甲烷在煤表面的吸附是放热过程,温度升高会导致Langmuir体积略微下降、Langmuir压力明显上升——换句话说,温度越高,同样压力下的吸附量越少,解吸越容易。我在拟合高温高压吸附实验数据时发现,如果忽略温度对Langmuir参数的影响,在60°C、10 MPa条件下吸附量的计算误差可以达到15%~20%,这在工程上已经不容忽视了。
其次是吸附/解吸本身伴随的热效应。解吸吸热,会让局部温度下降,反过来影响解吸速率;注热增产(如热蒸汽驱)则是主动利用这个效应。这些过程需要耦合能量守恒方程,考虑热传导、对流换热以及吸附热源项。
最后是热应变。温度变化引起煤基质膨胀或收缩,直接作用于力学场中的应变项。深层煤层地温梯度大约3°C/100米,1500米埋深的地层温度就有50°C上下,和室内实验25°C的温差产生的热应变不可忽略。
三个场的耦合关系可以大致归纳为:压力场通过有效应力影响变形,变形通过渗透率变化反馈到流动;温度场通过Langmuir参数影响气源和解吸,通过应变影响变形,同时流动过程又搬运热量。任何一个方向的忽略,在特定工况下都可能让模型预测偏离实际。
3. 渗透率-孔隙度模型的选型逻辑:从立方定律到主流模型的差异
如果说THM框架是"骨",那渗透率-孔隙度模型就是"血"。模型选得对不对、参数定得准不准,直接决定模拟结果能不能反映现场。
3.1 基础桥梁:孔隙度与渗透率的立方定律关系
煤储层渗透率的主体贡献来自裂隙网络,而裂隙渗透率和裂隙开度之间存在经典的立方定律关系:单条裂隙的渗透率正比于开度的三次方。假设裂隙开度的变化正比于孔隙度的变化,就可以推导出常用的关系式:
k/k0 = (φ/φ0)^3
这个式子简洁到让人觉得"太理想了",但在工程上它出奇地好用。理由其实不复杂:煤储层裂隙开度只有几十微米量级,孔隙度每变化零点几个百分点,对应的开度变化已经足以让渗透率产生数量级上的改变。立方定律把这种非线性放大效应表达得刚好到位。我在实际项目中通常以它为基础框架,再根据具体储层特征选择修正形式。
孔隙度的演化方程则来自力学场:总应变等于孔弹性应变、吸附/解吸应变、热应变的叠加,通过体积应变的变化来更新孔隙度。这就是所谓"动态孔隙度"的内核——把力学方程算出来的变形反馈到流动方程的存储项和渗透率项中。
3.2 三个主流模型对比:P-M、S-D、C-B到底差在哪
目前工程和科研中用得最多的三个渗透率动态模型是Palmer-Mansoori模型(P-M)、Shi-Durucan模型(S-D)和Cui-Bustin模型(C-B)。它们都试图回答同一个问题:压力变化+解吸应变作用下,渗透率如何演化?但各自的理论切入点和适用边界差异明显。
| 模型 | 理论基础 | 关键假设 | 核心表达 | 优势 | 局限 |
|---|---|---|---|---|---|
| P-M(Palmer-Mansoori, 1998) | 孔弹性+吸附应变,结合立方定律 | 单轴应变、恒定上覆应力、线弹性 | k/k0 = [1 + A·ΔP + B·Δε_sorption]³ | 形式简洁,参数物理意义明确,工程应用最广 | 未显式考虑温度项;对水平应力比变化敏感 |
| S-D(Shi-Durucan, 2004) | 应力-渗透率指数关系 | 单轴应变,应力变化中含基质收缩贡献 | k/k0 = exp(-3c_f·Δσ_eff) | 能较好匹配现场排采数据,指数形式对有效应力变化敏感 | c_f取值往往需要历史拟合反推,物理意义弱化 |
| C-B(Cui-Bustin, 2005) | 各向同性线孔弹性完整推导 | 既支持单轴应变也支持恒定围压条件 | 基于Biot孔弹性理论完整求解应变场后代入立方定律 | 理论基础最严谨,适用范围广 | 方程结构复杂,参数需求多,工程推广阻力大 |
从我的项目经验看,选型建议如下:
- 常规煤层气排采模拟、历史拟合为主:优先P-M模型。它的参数(割理压缩系数c_f、初始孔隙度φ0、最大收缩应变ε_L、Langmuir压力P_L)都能从常规实验获得,工程师友好度高。
- 重点研究渗透率急剧变化、近井压降带分析:S-D模型往往拟合效果更好,因为它用指数关系放大了有效应力变化的影响。但要注意c_f的取值范围,超出0.05~0.6 MPa⁻¹就要警惕物理合理性。
- 研究尺度复杂(如三维应力状态、应力反转):上C-B模型。不过要做好参数标定工作量翻倍的准备。
3.3 模型参数获取的工程路径
模型的可靠性完全取决于参数质量。我整理了一套参数获取的推荐路径:
- Langmuir参数V_L和P_L:通过平衡条件下不同压力点的等温吸附实验获取。但必须注意:实验室用的是干燥样,现场是含水条件;而且温度不同结果差异明显。有条件就做原位温度下的吸附实验,否则需要用热力学方法对参数做温度修正。
- 最大收缩应变ε_L:通过煤粒或煤块的解吸应变实验测得,一般在0.005~0.04之间。这个参数对渗透率回升幅度影响极大,实测时要注意吸附平衡时间,测试周期通常需要两周以上。
- 割理压缩系数c_f:由变围压渗透率实验拟合得到。煤样在加载过程中存在不可逆损伤,所以加载路径要和现场应力路径尽量一致。
- 弹性模量E、泊松比ν:通过三轴压缩或声波测井获取。煤的非均质性很强,同一层段不同钻孔的E可能差一倍,取均值时要和地质统计配合。
- 初始渗透率k0和初始孔隙度φ0:初始渗透率用试井资料;孔隙度用测井解释或岩心分析。这两个参数是模型标定的"锚点",优先级最高。
参数获取有一条总原则:能用实测的不用拟合的,能用试井的不用岩心的,能用原位条件的不用实验室条件的。
4. 数值实现:耦合策略、边界条件与参数标定的实际工程经验
模型选好、参数备齐之后,照理说就是"开算"了。但数值实现环节才是真正拉开差距的地方。同一套THM方程,不同人去实现,结果可能天壤之别。
4.1 耦合策略:顺序耦合还是全耦合,以及折中路线
THM三个场的数值求解有两种经典思路:
- 全耦合(Monolithic):把力学方程、流动方程、能量方程整合成一个大刚度矩阵,使用Newton-Raphson迭代同步求解。好处是稳定性好、无时间滞后误差,代价是内存占用大、Jacobian矩阵结构复杂、调试困难。对常规的单井排采模拟,性价比不高。
- 顺序耦合(Staggered/Sequential):每个时间步内先解渗流场,得到压力分布后传给力学场求变形,再由变形更新渗透率和孔隙度,接着进入下一时间步。实现简单、模块化程度高,可以用成熟的力学求解器和流动求解器拼接。缺点是时间步长过大时可能出现渗流-变形"不同步"导致的数值振荡。
我实际用过多种组合方案,一个折中路线是迭代耦合:每个时间步内顺序求解两到三轮,每轮检查渗透率和孔隙度的收敛残差,残差过大就重复迭代,直到满足容差。这个策略在COMSOL和FLAC3D接口方案里都验证过,稳定性接近全耦合,计算量却只有全耦合的一半左右。
软件选型方面:COMSOL Multiphysics做THM耦合很方便,自带固体力学、流体流动和传热模块,内置耦合接口;TOUGH+FLAC3D组合是地热和天然气水合物领域的老牌方案,但对工程师的代码能力有要求;Abaqus的孔压单元也能做,只是对吸附应变这类材料行为的二次开发比较繁琐。如果只是想快速验证机理,MATLAB里的偏微分方程工具箱也能搞定简化版THM模型。
4.2 初始应力场与边界条件:90%的模型误差从这里来
很多THM模型跑出来结果离谱,根子不在方程,而在初始条件和边界条件设错了。
第一个高频错误是初始应力场设置过于简化。煤储层通常处于三向应力不等的状态:垂向应力和埋深相关,水平应力受构造运动影响,往往两个水平主应力方向差异明显。模拟时不能只给一个"各向同性"的初始应力,否则渗透率各向异性根本体现不出来。我处理这类问题会先做一次应力平衡(也称"地应力回归"),让模型在无扰动条件下达到稳定的初始地应力状态,再开始排采模拟。
第二个高频错误是边界条件类型选错。储层模型常规的顶部边界应该是应力边界(相当于上覆岩层加载),而不是位移边界。用位移固定边界会导致上方岩层不能沉降,人为制造了过大的水平应力,渗透率被严重低估。侧向边界在对称条件下可以用辊轴约束,但底部边界要考虑热效应时就不能简单用绝热边界。
第三个问题是井筒网格尺寸。煤层渗透率低、压降带集中,井筒附近需要对数加密网格。我在一个项目中把近井网格从10米加密到1米后,累计产气量预测变化了15%,而计算时间才增加了几分钟。这个投入产出比太划算了。
4.3 参数标定:历史拟合的完整工作流
历史拟合不是"调参数调到曲线对上",而是有章法地缩小不确定性。我的标准工作流大致是:
- 先找锚点参数:初始渗透率、初始孔隙度、初始压力梯度、解吸压力,这些有实测依据的参数一律锁死,不允许参与调参。
- 敏感性排序:对候选参数做一次全局敏感性分析,用拉丁超立方采样跑几十个案例,识别出对产量曲线影响最大的前五个参数。通常包括ε_L、c_f、基质块尺寸、相对渗透率曲线形态、表皮因子。
- 分阶段拟合:先拟合产水阶段(确定渗透率和表皮因子),再拟合产气启动阶段(确定解吸参数),最后拟合稳产和递减阶段(确定收缩参数和边界影响)。不要一上来就全方位拟合。
- 多目标验证:历史拟合不只是匹配产量曲线,还要验证压力点、饱和度剖面和渗透率随时间的变化趋势。如果压力和产量不能同时匹配,说明模型结构有问题而不是参数问题。
这套流程走下来,模型的不确定性通常能控制在工程可接受范围内。
5. 踩坑实录:从模型算到现场验证的几个教训
再好的理论模型放到现场都会有落差。这几个坑我踩过不止一次,写下来给后来者省点时间。
5.1 实验室渗透率与现场渗透率的"系统偏差"
室内岩心实验测出来的渗透率经常比现场试井解释值低一个数量级,甚至更多。原因有几个:取心过程中煤样卸压导致割理不可逆张开损伤;实验室测试通常在高围压下进行,和原位应力状态不一致;煤样经过干燥处理,含水饱和度的变化改变了有效渗透率。
我踩过的坑是:拿实验室渗透率直接当作模型的初始渗透率,结果历史拟合怎么也调不上去。后来意识到该用试井数据作为锚点,实验室数据只用来确定渗透率变化趋势(比如应力敏感系数),而不是绝对值。这个思路调整之后,拟合收敛速度明显加快,预测结果也靠谱得多。
5.2 忽略温度项的代价:深层项目中的"系统性偏差"
有一次我负责一个埋深1200米项目的机理建模,为了"简化起见"用了等温假设,把地温固定在45°C。结果模型预测的累计产气量比实际单井数据偏低20%左右。排查了很久,最终定位到问题出在等温假设:实际排采过程中,近井区压力快速下降诱发大量解吸,解吸吸热导致局部温度下降3~5°C,而温度降低反而增强了吸附能力(Langmuir体积回升),形成"负反馈",最终解吸速度比等温模型预测的快。之后我把温度场完整加上,这个偏差消失了大半。
深层煤层气项目的温度效应远不是"加一个热膨胀系数"那么简单,它是通过Langmuir参数温度修正、吸附热源项、热膨胀应变三条路径同时作用的。
5.3 判断模型结果合理性的快速检查清单
经过几个项目的打磨,我总结了一套"算完先自查"的流程,分享出来:
- 质量守恒检查:累计产气量是否与模型质量平衡项一致,偏差大于2%就要查数值误差。
- 渗透率变化幅度:模拟的k/k0变化范围是否落在0.1~10倍区间内。如果出现渗透率放大几十倍的情况,多半是参数(通常是ε_L或c_f)取值越界了。
- 压力梯度合理性:近井压降漏斗形态是否与试井解释结果一致,远场压力是否保持稳定。
- 饱和度分布逻辑:含水饱和度剖面在排采过程中是否出现违反相对渗透率物理规律的突变。
- 无外力振荡:产量曲线和压力曲线在时间序列上是否平滑,出现锯齿状多半是时间步长或迭代容差设置不合理。
这套检查清单虽然不能完全替代严格的验证工作,但可以拦截掉大多数"模型算完感觉不对"的问题。
6. 下一步我打算尝试的方向
THM耦合在煤层气领域还有不少值得深入研究的方向,简单说说我个人的计划。
一是多场耦合向化学场延伸。煤气层水中的溶解气(CO2、H2S)以及煤岩与水长期作用下的矿物溶解—沉淀过程,都会改变裂隙系统的渗流能力。THM升级为THMC(加入Chemistry)是必然趋势,虽然实现复杂度和参数需求会明显增加。
二是机器学习辅助参数反演。历史拟合本质上是一个高维优化问题,传统方法容易陷入局部最优。我最近在尝试用神经网络代理模型替代昂贵的THM正演计算,把敏感性分析和历史拟合的速度提升一个量级。初步结果不错,代理模型的预测误差控制在5%以内,而单次正演时间从几小时压缩到几十毫秒。
三是多尺度联合建模。微米尺度的纳米CT孔隙结构、岩心尺度的渗流实验、数十米尺度的试井解释,再到千米尺度的储层模拟,每个尺度的参数怎么传递和升级,是当前模型不确定性最主要的来源。现在常用的做法是经验性的"等效参数"升级,但这个路子在复杂储层的适用性越来越吃力。
最后分享一个实际工作中的小技巧:无论用多么完善的THM模型,每次模拟前都先跑一个极度简化的解析解或半解析模型做交叉验证。比如单相微可压缩流体的产量递减解析解、简单的渗透率-压力衰减关系式。结果对得上,再进入全耦合模拟。这个习惯帮我拦截了无数次由网格问题、边界问题引起的"隐性错误"——这些错误在参数曲面和产量曲线上往往隐藏得很深,肉眼根本看不出来。