煤层气藏全流固耦合数学模型

未志杰 ，康晓东 ，刘玉洋 ，曾 杨

（1.海洋石油高效开发国家重点实验室，北京100028；2.中海油研究总院有限责任公司，北京100028）

0 引言

煤层属于天然裂缝性储层，力学强度低，具有显著的应力敏感性。流固耦合作用，又称地质力学效应。研究表明，在数值模拟过程中是否考虑流固耦合作用，对于准确预测煤层气产能至关重要［1-2］。相比其他裂缝性储层，煤层的流固耦合作用更为复杂，不仅存在有效应力效应，还包括基质膨胀或收缩作用。基质因煤层气吸附或解吸现象产生形变，吸附使基质膨胀，导致有效渗流孔道和渗透率均减小；解吸使基质收缩，导致有效渗流孔道和渗透率均增大。

为刻画煤层流固耦合作用，研究人员［3-6］曾提出包括ARI模型、Palmer&Mansoori模型和Shi&Durucan模型等在内的多种模型。ARI模型［3］为经验公式，没有地质力学理论基础，基质膨胀或收缩作用引起的应变量与煤层气吸附量成正比；Palmer模型［4］基于地质力学理论，认为煤层是均质各向同性线弹性孔隙介质，并将基质膨胀或收缩应变等效类比为热膨胀应变；Shi模型［5］不同于 Palmer模型，渗透率与水平有效应力呈对数关系。这些模型均属于解析或经验流固耦合模型，具有形式简洁直观、便于与商业煤层气模拟软件结合的优势，但缺点是须引入较多强假设，如固定上覆应力与单轴向应变假设，导致渗透率计算结果失真，影响产能预测精度。由此通过引入煤岩形变本构方程来准确刻画煤层地质力学效应，即建立全流固耦合数学模型，以期获得更准确的储层物性参数及产量预测结果。

1 全流固耦合数学模型

煤层气藏全流固耦合数学模型，包括多过程流动方程和多孔介质地质力学方程。前者采用三孔双渗模型来描述煤层孔隙结构，可进一步发展目前的双孔单渗模型；后者通过计算煤层应力、应变等，可获得更准确的孔隙度和渗透率。

1.1 多过程流动方程

煤层属于非常规天然气储层，既是生气岩又是储气岩，原始地层条件下煤层气的主要赋存状态为吸附态而非游离态。目前常用双孔单渗模型来描述煤层［6］，即煤层由连续分布的裂缝系统以及离散分布的基质系统构成，相应的流动过程为地层孔隙压力随排水下降，煤层气由基质解吸并扩散进入裂缝，而后以渗流方式进入井筒。在矿场应用中，这一模型的计算结果与实际生产数据存在较大偏差（煤层气产量被高估、产水量被低估、气体突破时间比实际早）［7-9］，有必要对煤层多过程物质运移和煤层孔隙结构进行细分再认识，以得到合理准确的产量计算结果，为此进一步发展出三孔双渗模型。3套孔隙系统是指基质、中孔隙系统和裂缝系统等。基质离散分布，是煤层气主要的吸附和存储空间；中孔隙系统与裂缝系统虽分布连续，但物性差异显著，中孔隙系统孔隙度大而渗透率小，裂缝系统孔隙度小而渗透率大，中孔隙系统提供地层水存储空间，而裂缝系统提供主要的渗流通道。三孔双渗模型的流动过程为：首先进行排水降压，使煤层气解吸，基质内煤层气在浓度梯度驱动下扩散进入中孔隙系统；然后以窜流方式进入裂缝；最后以渗流方式进入井筒并采出。

1.1 .1 质量守恒方程

忽略煤层气在水中的溶解和水分挥发，建立质量守恒方程，由连续性方程和达西渗流方程构成。

连续性方程为

达西渗流方程为

式中：ρ为流体密度，kg/m3；下标l为流体相、η为孔隙系统、s为基质；Φ为孔隙度，%；S为流体饱和度，%；V为运动速度，m/s；q为源汇项，kg/（m3·s）；t为时间，s；K为渗透率mD，Kr为相对渗透率；μ为黏度，mPa·s；P为流体压力，Pa为梯度，·为散度。

将式（2）代入式（1），可得

将式（3）右边第二项展开，可得

同时根据定义，固相运动速度Vs与固相位移u的关系满足

煤层气开采过程中固相位移与固相运动速度相对较小，整理式（4）与式（5）后得

式中：u为位移，m；e为煤岩体应变。对于三孔双渗模型，质量守恒方程具体形式为

式中：下标w为水相、g为气相、f为裂缝、m为中孔隙系统；qgms为煤层气解吸量，kg/（m3·s）；qgmf和qwmf分别为煤层气和水在中孔隙系统与裂缝之间的窜流量，kg/（m3·s）；qgwell和qwwell分别为煤层气和水的源汇项，一般指开发井的生产量或注入量，kg/（m3·s）。

1.1.2 吸附及扩散方程

煤层气吸附采用Langmuir等温吸附曲线描述，基质表面的平衡吸附量为

式中：C（Pgm）为煤层气平衡吸附体积浓度，m3/m3；Pgm为中孔隙系统的煤层气压力，Pa；VL为基质对煤层气的Langmuir吸附量，m3/m3；PL为煤层气Langmuir吸附压力，Pa。

对于基质中煤层气的扩散运移，一般采用拟稳态非平衡模型进行描述，根据菲克第一定律，推导得出煤层气扩散方程

式中：C为基质中煤层气的平均浓度，m3/m3；τ是解吸附时间，s。

1.1.3 窜流方程

中孔隙系统与裂缝系统存在势差，在此作用下产生窜流或物质交换，中孔隙系统与裂缝系统之间的窜流量qgmf和qwmf的计算方程［9］分别为

式中：δ为形状因子，m-2。

1.2 地质力学方程

1.2.1 煤岩本构方程

根据Biot线弹性理论，对于各向同性弹性的煤层，本构方程为

式中：T 为应力张量，Pa；Te为有效应力张量，Pa；Γ为应变张量；I为单位张量；G为剪切模量，Pa；λ为平均轴向模量，Pa；es是基质体应变。

应变与位移的关系为

动量距平衡方程为

联立式（15）—（18），可得

对式（19）取散度，最终将转化成如下简洁形式?

煤层形变受制于有效应力与基质膨胀或收缩效应。实验表明，基质体应变量与气体压力的关系呈现Langmuir曲线特征［4］

式中：εL为基质Langmuir体应变量。

1.2.2 孔隙度与渗透率方程

煤层渗透率的影响因素包括有效应力效应、基质膨胀/收缩作用和气体滑脱效应等。有效应力升高导致渗透率降低［10-11］；煤层气吸附或解吸导致基质膨胀或收缩，同样可使渗透率发生变化［3］；此外，煤层渗透率低，气体滑脱效应同样不可忽略［12］。

通过数学推导，裂缝与中孔隙系统的孔隙度为

式中：Vpη为孔隙系统 η的孔隙体积，m3；Vb为煤岩总体积，m3。中孔隙系统与裂缝的孔隙体积变化方程［13］为

考虑到煤岩体应变dVb/Vb=de和平均主应力，孔隙度方程为

式中：cpn为孔隙体积压缩系数，Pa-1；β为有效应力系数；Kb为煤岩体积模量，Pa。

裂缝一般采用火柴棒束模型来描述，这样渗透率与孔隙度呈三次方关系

式中：Pkf为裂缝Klinkenberg系数，Pa。上式右端第一项表征有效应力和基质收缩效应，第二项表征气体滑脱效应。

对于中孔隙系统，火柴棒束模型已不再适用，采用Kozeny-Carman模型来表征渗透率与孔隙度的关系［14-16］，可表达为

式中：Pkm为中孔隙系统Klinkenberg系数，Pa。根据实验结果［17-20］，Pkf和 Pkm分别取值 0.28 MPa 和0.84 MPa。

2 数值求解

根据全流固耦合数学模型，控制方程包括质量守恒方程［式（7）～（10）］、扩散方程［式（12）］和地质力学方程［式（20）］等，主要未知量包括水相压力（Pwf和 Pwm）、水相饱和度（Swf和 Swm）、基质煤层气平均吸附浓度（C）和煤岩体应变（e）等。采用全隐式有限体积法离散质量守恒方程［21-24］，得

将式（27）展开，可得

对于扩散方程式（12）和地质力学方程式（20），有

采用 Newton-Raphson 迭代全隐式方程［25-28］求解上述离散方程，可获得三孔双渗全流固耦合煤层气模拟器。主要未知量 Pwf，Pwm，Swf，Swm，C，e 被同时解出，据此同步更新煤层孔隙度/渗透率、流体PVT属性、流度和传导系数等，进而获得气、油产量，然后进入下一个时间步计算。

3 模型及算法验证

通过与商业模拟软件Coalgas进行对比分析，来验证三孔双渗全流固耦合数学模型及模拟器的准确性。因二者所采用的流固耦合模型不同，本次研究为全流固耦合数学模型，而Colagas使用的是ARI模型，所以为使计算结果具有可比性，可忽略煤层地质力学作用，即使孔隙度和渗透率保持不变。

煤层概念模型的输入参数为：煤层长度为360 m、宽度为360 m、厚度为5 m；网格平面划分为51×51，几何结构对称，生产井位于中心并保持定压生产；初始地层条件下裂缝系统孔隙度为0.5%，渗透率为5.0 mD，中孔隙系统孔隙度为5.0%，渗透率为0.05 mD；煤层气吸附参数、流体PVT属性参数、煤岩力学参数、裂缝系统与中孔隙系统物性参数如表1所列；煤层初始状态是非饱和吸附，初始孔隙压力为6.895 MPa，临界解吸压力为5.17 MPa。

煤层气与水的产量模拟结果如图1所示，模拟器计算结果与Coalgas的模拟结果一致，初步验证了全流固耦合数学模型的准确性和可靠性。

4 流固耦合效应的影响

流固耦合效应包括有效应力效应和基质收缩作用，下面分别剖析这2种作用对渗透率和煤层气生产动态的影响。

4.1 有效应力效应

为了评估有效应力效应的影响，可暂且忽略基质收缩作用，这样煤层在力学上等效于普通裂缝性储层。有效应力效应作用程度依赖于煤岩力学强度，而杨氏模量是衡量煤岩力学强度的关键参数。本次研究选择了3种不同力学强度的煤层，对应杨氏模量E值分别为 1.25 GPa，2.50 GPa和 5.00 GPa，同时选取刚性储层作为对照。不同力学强度储层的煤层气产量如图2所示。可知弹性储层的煤层气产量与刚性储层存在明显差异，弹性储层的煤层气高峰日产量较低，且高峰出现时机较晚，这是有效应力效应作用的结果。杨氏模量越大，煤层气高峰产量越大，且出现时机越早。

表1 煤层气模拟输入参数Table 1 Simulated input data of coalbed methane

图1 全流固耦合模拟算法与Coalgas模拟结果对比Fig.1 Comparison of simulation results between the fully coupled fluid flow and geomechanics model and Coalgas

图2 不同杨氏模量煤层的产气量对比Fig.2 Comparison of gas rate for reservoirs with different Young's modulus

随着煤层气的开采，地层孔隙压力减小，有效应力（压应力）增大，导致孔隙收缩，孔隙度持续下降，而渗透率因气体滑脱作用的影响而呈现不同的变化趋势。以裂缝渗透率为例，不同杨氏模量下的渗透率随孔隙压力的变化情况如图3所示，图中对渗透率进行了归一化处理。对于刚性煤层，当孔隙压力由初始值（6.895 MPa）降低至临界解吸压力（5.170 MPa）之前，孔隙中无游离气，渗透率保持不变；之后，煤层气由基质中解吸，渗透率受气体滑脱作用影响而不断升高。对于弹性煤层，渗透率的变化受有效应力和气体滑脱作用双重作用影响，且二者作用方向相反，前者导致渗透率减小、后者导致渗透率增大：有效应力效应在开采初期占主导，渗透率减小，但后期气体滑脱作用增大并逐渐占优，渗透率可止降为升。

图3 不同杨氏模量煤层井网格的裂缝渗透率变化情况对比Fig.3 Comparison of normalized fracture permeability of well block for reservoirs with different Young's modulus

图4 不同Langmuir体应变量煤层的产气量对比Fig.4 Comparison of gas rate for reservoirs with different Langmuir volumetric strains

4.2 基质收缩作用

为评估基质收缩作用的影响，保持煤岩杨氏模量不变，通过改变基质形变强度，相应Langmuir体应变量 εL分别取 0，0.006，0.012与 0.024等，煤岩杨氏模量为2.50 GPa。εL越大，单位煤层气吸附量引起的基质收缩量越大，对渗透率与孔隙度的影响也越大，εL=0即不考虑基质应变。

不同基质形变强度下的煤层气产量如图4所示：Langmuir体应变量εL对煤层气产量的动态影响显著，εL越大，峰值产量越大，且出现时间越早。可见对煤层气产量的影响而言，基质收缩作用与有效应力效应作用方向相反。

不同基质形变强度下裂缝系统与中孔隙系统的渗透率变化情况如图5所示，当εL≠0时，渗透率变化受有效应力效应、基质收缩作用和气体滑脱效应等3种作用影响，其中有效应力效应会使渗透率减小，但基质收缩作用和气体滑脱效应却使渗透率增大。当孔隙压力高于临界解吸压力时，煤层气无解吸，仅有效应力效应发挥作用，故渗透率不变；随着煤层气解吸，诱发基质收缩作用和气体滑脱效应，在三者共同作用下，裂缝渗透率先减小后增大，εL=0.024时，最终渗透率甚至可达初始值的2.2倍。中孔隙系统渗透率仅在初期微弱下降，而后快速上升至初始值的2.3倍，这主要是气体滑脱作用的影响。裂缝渗透率受基质Langmuir体应变量εL的影响显著，但对中孔隙系统的影响几乎忽略不计，这是因为中孔隙系统的力学强度大于裂缝系统（cpf＞cpm）。

图5 不同Langmuir体应变量煤层井网格的渗透率变化情况对比Fig.5 Comparison of normalized permeability of well block for reservoirs with different Langmuir volumetric strains

5 结论

（1）构建了煤层气藏全流固耦合数学模型，包括多孔介质地质力学方程和三孔双渗流动方程，前者定量刻画了煤岩形变，获得了更准确的孔渗参数；后者进一步完善了目前的双孔单渗模型，实现了煤层气多过程物质运移的准确表征。

（2）开发了基于全隐式有限体积法的全流固耦合、三孔双渗煤层气数值模拟器，通过与商业模拟软件Coalgas进行对比，验证了其准确性与有效性。

（3）地质力学效应对煤层气生产动态预测至关重要。有效应力效应与基质收缩作用均可显著影响孔渗参数及煤层气产量，前者降低渗透率，后者增大渗透率，二者作用方向相反，但通常不可抵消。有效应力效应的作用强度在开发初期大于基质收缩作用，后期逆转，导致裂缝渗透率呈现独特的“U”型变化，最终值甚至可达初始值的数倍。