适用于深层页岩气井生产模拟和预测的自扩散流热耦合模型

夏 阳 韦世明 金 衍 陈康平

1. 油气资源与探测国家重点实验室·中国石油大学（北京） 2. 美国亚利桑那州立大学

0 引言

四川盆地及其周缘埋深介于3 500 ～4 500 m 的深层页岩气资源丰富，占该区域页岩气总资源量的85%以上[1]。深层页岩气储层具有高温、高压、高地应力的特点[2]，准确预测深层页岩气井产能可以为气藏开发方案的制订提供重要的依据。

页岩气井生产数据统计结果表明，采用达西渗流模型计算得到的气井产量往往低于实际气产量[3]。现有的页岩气渗流模型则在达西渗流模型的基础上，增加考虑滑脱、扩散等因素的影响[4-7]，对达西渗流模型进行修正。Klinkenberg[8]最早通过岩心渗透率测试实验，发现同一块岩心的气测渗透率大于其液测渗透率，并且气体压力越低，气测渗透率越高；气测渗透率与气体压力倒数近似成正比，此即为气体滑脱效应，又称为Klinkenberg 效应。2007 年，Javadpour等[9]将努森扩散引入页岩气渗流模型中，并给出页岩孔隙中努森扩散系数的具体形式。2015 年，吴克柳等[10]采用努森数定义了页岩气流动过程中努森扩散和黏性流的权重系数，发现努森扩散在低压条件下对页岩气渗流的贡献不可忽略。由于滑脱效应仅发生在稀薄气体中，而实测气体滑脱因子介于0.1 ～1.0 MPa[11-13]，对于深层页岩气储层而言，由于地层压力较高，滑脱效应对气体流动影响很小。Patzek[14]、Chen 等[15]指出，在深部致密气储层中气体滑脱效应和努森扩散可以忽略。

页岩储层基质渗透率极低，应力敏感效应使得生产过程中渗透率将进一步降低[16-17]。由于深层页岩储层压力、温度高，在气井生产过程中气体流动极大地依赖其自身膨胀，国内外大多数学者将自扩散效应作为页岩孔隙中气体的主要流动方式，并建立了初步的气体自扩散流动模型[18-21]。Larry[22]研究了溶解气驱油流动机理，发现气体膨胀是原油流动的主要驱动力。Kolesar 等[23]，Carlson 等[24]、Amann-Hildenbrand 等[25]随后提出扩散效应是致密气流动的主要方式。Wan 等[20]基于气体扩散系数计算方法，模拟页岩油吞吐生产过程中的油气产量，对实验结果进行了很好地拟合。Cronin 等[21]基于Fick 定律，建立了不同气体在页岩油吞吐生产过程中的质量扩散方程，计算的油气产量与现场生产数据吻合程度高，但模型中扩散系数是基于对生产数据进行拟合，尚未建立普适的自扩散流动模型。Hoteit[26]提出基于浓度梯度的Fick 定律不能描述致密储层中的气体扩散行为。Jin等[27]从可压缩流体的N-S 方程出发，在Klainerman等[28]提出的小马赫数流动理论基础上，建立了致密气一次开采气体自扩散流动控制方程，并基于现场实例对自扩散流动模型的特征和适用性进行了研究。

深层页岩储层具有高温、高压的特征，CH4的体积黏度远大于其剪切黏度且受温度的影响较大[29]。在页岩气井生产过程中，由于压裂液冷却、气体膨胀使得近井温度降低，温度的变化将对气体流动产生影响，进而影响深层页岩气井的产能[30]。截至目前，国内外关于流体自扩散流动模型的研究都未考虑高温、高压的影响，尤其是温度变化对气体物理化学性质及传质速度的影响规律尚未见到阐明。

为此，笔者基于Jin 等[27]提出的自扩散流动模型，建立了考虑渗流场和温度场耦合的自扩散流动模型（以下简称自扩散流热耦合模型），然后，将该模型与以达西定律和努森扩散为基础的渗流模型（以下简称修正达西模型）计算结果进行对比；基于自扩散流热耦合模型，分析考虑温度场变化的深层页岩气自扩散流动行为，进而探讨了井底温度变化对深层页岩气采出程度的影响；在此基础上，将自扩散流热耦合模型应用于四川盆地长宁区块上奥陶统五峰组—下志留统龙马溪组一口页岩气水平井的生产模拟，以期为页岩气井产能的准确预测提供支撑。

1 自扩散流热耦合模型的建立

在地层高温高压条件下，气体状态方程为：

式中ρG表示气体密度，kg/m3；p 表示气体压力，Pa；M 表示气体分子摩尔质量，kg/mol；Z 表示气体偏差因子；R 表示普适气体常数，取值为8.314 J/(mol·K)；T 表示气体温度，K。

气体偏差因子和剪切黏度随气体组成、压力和温度发生变化，笔者采用Heidaryan 等[31]提出的经验公式，即式中ppr表示气体对比压力，取值为1；Tpr表示气体对比温度，取值为1；A1—A11表示由实验数据拟合得到的常系数，依次取值为3.252 838、-1.306 424× 10-1、-6.444 919 4×10-1、-1.518 028、-5.391 019、 -1.379 588×10-2、6.600 633×10-2、6.120 783×10-1、2.317 431、1.632 223×10-1、5.660 595×10-1；μ表示气体剪切黏度，mPa·s；B1—B10表示由实验数据拟合得到的常系数，依次取值为1.022 872、 -1.651 432、5.757 386、-7.389 282×10-2、8.389 065× 10-2、2.977 476 10-1、-1.451 318、4.682 506、1.918 239、 -9.944 968×10-2。

当考虑温度场的变化时，气体体积流量与热通量可分别表示为[32]：

式中qG表示通过岩石单位面积的气体体积流量，m/s；Dj表示自扩散系数，m2/s；βT表示热渗透系数，m2/(s·K)；qT表示通过岩石单位面积的热通量，W/m2；CG表示气体比热容，kJ/(kg·K)；φ 表示基质孔隙度；kG、ks分别表示气体、岩石的热传导系数，W/(K·m)。

Dj计算式为[27]：

式中μb表示气体体积黏度，mPa·s。

图1 甲烷与T 关系曲线图[29]

多孔介质中流体的质量守恒方程为：

含流体的多孔介质热平衡方程为[30-31]：

式中ρs表示岩石密度，kg/m3；Cs表示岩石比热容，J/(kg·K)；αG表示气体热膨胀系数，K-1；KG表示气体体积模量，Pa。

将式（1）、（4）代入式（7），则气体自扩散流动控制方程为：

将式（4）、（5）代入式（8）中，得到温度场控制方程[30-31]，即

其中

将式（9）、（10）结合，即为自扩散流热耦合模型，其中待求解变量为p 与T。该模型仅将游离气的膨胀考虑为气体流动的驱动力，暂未考虑吸附气的影响。

2 模型分析

笔者针对压裂井模拟区域建立了一维、二维模型（图2），模型参数如表1 所示。一维模型的左端代表裂缝面，其截面积为1 m2，边界条件为左端定压定温、右端封闭。笔者采用有限元软件COMSOL中的数学模块对式（9）、（10）的偏微分方程进行求解。通过对比自扩散流热耦合模型（考虑温度场的影响）与修正达西渗流模型的计算结果，明确自扩散流热耦合模型的正确性；在此基础上，分析温度场对气体自扩散流动行为的影响，进而研究温度场变化对深层页岩气井产量的影响。

图2 压裂井模拟区域划分示意图

表1 模型参数表

2.1 计算结果对比

常规渗流模型认为气体流动受到压力梯度和浓度梯度的影响。基于达西定律，增加考虑努森扩散的影响，流量计算式为[6-7]：

式中J 表示通过单位岩石截面积的气体质量流量，kg/(m2·s)；K 表示岩石渗透率，mD；D 表示努森扩散系数，m2/s；M 表示气体摩尔质量，g/mol；C 表示气体浓度，mol/m3。

仅使用式（11）代替式（4）右侧自扩散流量时，就可以得到修正达西模型，之后可以采用有限元数值模拟软件COMSOL 中的数学模块进行求解。本节基于一维模型（图2）进行相关计算。

图3 自扩散流热耦合模型与修正达西模型计算孔隙压力剖面对比图

图4 自扩散流热耦合模型与修正达西模型计算气体密度剖面对比图

通过对比自扩散流热耦合模型与修正达西模型计算的压力剖面，发现在同一时间，自扩散流热耦合模型计算的压力扩散速度更快（图3）。如图3、4 所示，采用修正达西模型计算得到压力剖面和气体密度剖面，压力、气体密度均在裂缝面处最小。由于气体密度、浓度与压力正相关，因而图3 中基于自扩散流热耦合模型计算的压力剖面与修正达西模型计算得到的压力剖面变化趋势一致。然而，由于温度场的影响，自扩散流热耦合模型计算得到的气体密度剖面在基质内部出现最小值，且该最小值对应的位置随着生产时间延长向基质内部移动（图4），可以看出，气体自扩散并非仅受气体浓度梯度的影响[33]。如图5 所示，采用自扩散流热耦合模型计算的气体累计产出质量远大于修正达西模型计算值，说明基于气体膨胀驱动计算的产气量大于基于气体压力和浓度梯度驱动计算的产气量，从而解释了页岩气井生产初期气产量通常大于修正达西模型计算值的原因[3]。

图5 自扩散流热耦合模型与修正达西模型模拟 气体累计产出质量对比图

2.2 温度场变化对气体自扩散流动行为的影响

下面分析温度场变化对深层页岩气自扩散流动行为和产量的影响，本节计算采用的模型参数、边界条件与2.1 部分相同，对于恒温地层，则考虑地层温度保持不变。

如图6 所示，若不考虑温度场变化的影响，自扩散系数沿裂缝面向基质内部逐渐减小；考虑温度场变化后，自扩散系数沿裂缝面向基质内部呈先增大后减小的趋势，且随着生产进行，自扩散系数峰值所在位置逐渐向基质内部移动。使用气体状态方程，将自扩散系数计算式变形为：

式中bg表示气体体积黏度与剪切黏度的比值。

由式（12）可以看出，Dj随着T 升高而增大，随着p 升高而减小。在裂缝面温度低于储层初始温度时，从裂缝面向基质内部，地层温度和孔隙压力均不断升高，自扩散系数呈先增大后减小的趋势。

不考虑温度场变化时，自扩散流热耦合模型计算得到的气体密度剖面（图7）与修正达西模型计算得到的气体密度剖面形状相同（图4）。如图7 所示，随着地层温度降低，从裂缝面向基质内部，气体密度先降低后增大，其原因是井底温度的降低使得自扩散系数降低，造成气体流动速度在裂缝面附近减小，引起气体在近裂缝面区域发生质量堆积，即气体自扩散的近井阻塞；随着生产时间进行，低温区域范围扩大，发生阻塞的区域范围增大，气体密度谷值所在位置向基质内部移动。

如图8 所示，考虑井底温度低于地层温度时，预测的气体累计产出质量低于恒温地层。由此可知，近井低温将使页岩气井产量降低。

3 井底温度变化对深层页岩气采出程度的影响探讨

基于第2 节的模型参数，在裂缝面设定不同的温度，分别低于初始地层温度10 ～80 ℃，计算生产3 年末的气采出程度（累计产气量与模拟区域内页岩气储量的比值）。如图9 所示，井底温度越低，即温差越大，气采出程度越低；随温差增大，采出程度近直线下降，温差从10 ℃增至80 ℃，采出程度降低2.3%。在压裂过程中，压裂液会冷却裂缝面，并且随压裂液注入量增大和焖井时间延长，裂缝面温度还会进一步降低。由于深层页岩气储层具有高温、高应力特征，裂缝宽度更窄，在页岩气开采过程中，上述两方面的影响会更明显，因而对于深层页岩气的开采，应该考虑温度场变化所产生的影响。

图9 温度差对气井生产3 年末气采出程度影响曲线图

4 实例井应用

选取四川盆地长宁区块五峰组—龙马溪组一口页岩气水平井进行生产模拟，该井基础参数如表2 所示，岩石和气体热力学参数如表1 所示。建立压裂井二维模型（图2），假设水力裂缝为无穷导流能力裂缝，则水力裂缝处边界条件为定压定温，其中水力裂缝处的压力则为由井口油压计算得到的井底压力[34]。

表2 实例井基础参数表

如图10 所示，考虑温度场变化影响的模拟日产气量与实际生产数据更加吻合。刚开始生产时，若不考虑温度场的变化，近裂缝处自扩散系数计算值将大于实际值，从而使模拟日产气量远高于实际结果。2014 年6 月18 日—7 月1 日，关井使得该井地层压力得以恢复；由于井筒储存效应的影响，在开井复产初期，实际日产气量高于模拟结果。2014 年8 月12日以后，关井压力恢复对产量模拟产生的影响消失，且考虑温度场变化模拟的日产气量与实际日产气量的吻合度更好。在压裂过程中，压裂液将冷却裂缝面，并且裂缝面温度与压裂液注入量、压裂液浸泡时间均相关，而且裂缝面温度与初始地层温度相差越大，日气产量越低。因此，对于采用大液量压裂的深层页岩气井来说，必须考虑温度场变化的影响才能准确预测单井气产量。

5 结论

1）在相同参数下，采用自扩散流热耦合模型计算的页岩气产量大于修正达西模型计算值。

2）考虑温度场变化时，自扩散系数剖面出现峰值，气体密度剖面出现谷值，并且峰值/谷值对应数据点同步向地层内部移动，说明自扩散系数影响着气体传质速率，且随生产进行，近井低温对气体自扩散的影响范围不断增大。

图10 长宁区块五峰组—龙马溪组页岩气水平井生产模拟结果图

3）当井底温度低于地层温度时，近井气体自扩散系数减小，发生气体近井阻塞现象，从而使气井产量降低。

4）根据现场实例井的生产模拟结果，考虑温度场的变化，自扩散流热耦合模型能够更准确地模拟深层页岩气井产量。