近年来，国内外页岩气勘探开发取得了较快发展^[1-3],但页岩气藏渗流机理复杂^[4]，不仅涉及解吸、扩散和对流机理，而且伴随着气体滑脱和应力敏感效应。目前，关于滑脱和应力敏感效应对气藏储层物性和产能影响的研究已经取得了一定进展和定性认识，主要集中在致密低渗透气藏及煤层气藏的滑脱效应和页岩气渗流理论方面^[5-9]：J.E.Warren等人^[7]认为储层孔隙压力较低时，气体滑脱效应较强；V.Shabro和V.Swami等人^[10-11]认为滑脱效应使页岩气井生产时储层的视渗透率增大，从而使页岩气井产能升高；P.N.Azom 等人^[12-14]认为建立页岩气渗流模型时，应该考虑滑脱效应。此外，夏阳等人^[15]的研究结果表明，页岩存在强应力敏感效应，导致页岩气井产能降低。由此可见，滑脱和应力敏感效应是影响页岩气藏开发动态的2个重要因素，亟需开展其对页岩气开采影响的定量研究。

目前，页岩气井产能的定量研究通常采用解析法^[16-19]和数值模拟方法^{[6, 7, 20-22]}，由于解析法计算产能时无法考虑气层非均质性等复杂因素的影响，因而，数值模拟方法成为研究页岩气井产能的主要手段。目前，页岩气井产能数值模拟主要采用双重孔隙介质、多重孔隙介质和等效孔隙介质模型^[20-21]。E.S.Carlson等人^[22]的研究结果表明，采用双重孔隙介质模型预测页岩气井长期产能时的结果偏低。虽然采用多重孔隙介质模型描述渗流介质在理论上是合理的，但是存在介质类型划分标准不统一和参数确定困难等问题，因此，采用等效孔隙介质模型预测页岩气井产能是一种比较可行的手段^[23]。文献调研表明，国内外尚无成熟可靠的页岩气藏数值模拟软件，为了定量研究页岩气藏开发动态和产能，笔者运用等效孔隙介质模型，建立了同时考虑滑脱和应力敏感效应的页岩气藏渗流数学模型，采用有限差分数值离散方法进行求解，并采用Fortran计算机编程语言开发了页岩气藏数值模拟器。利用该模拟器模拟研究了滑脱与应力敏感效应对页岩气井产能的影响规律，得出了页岩气井产气量的变化规律，并与生产井实际生产数据进行对比，验证了该模拟器的可靠性。

1 滑脱和应力敏感效应 1.1 页岩气的滑脱效应

因气体在孔道壁表面的切向速度不为零,气测视渗透率随着压力升高而降低的现象称为气体滑脱效应，Klinkenberg建立了考虑滑脱效应的克氏渗透率计算公式^[24]：

(1)

式中：K_a为气体的渗透率，D；K₀为孔隙介质的本征渗透率，D；p为目前气藏的压力，Pa；b_k为滑脱系数，Pa。

1.2 页岩的应力敏感效应

页岩气开采过程中，气藏压力降低使岩石骨架受到的有效应力增大，从而导致储层的渗透率、孔隙度降低的现象称为页岩应力敏感效应。页岩的渗透率、孔隙度与有效应力间的关系式为^[25-26]：

(2)

(3)

式中：p₀为气藏原始压力，Pa；K(p)为压力p下的气藏渗透率，D；φ和φ₀分别为p和p₀对应的孔隙度；α_K和α_φ分别为渗透率和孔隙度变化系数，Pa^-1。

2 页岩气藏渗流数学模型的建立

目前，人们在建立页岩气藏渗流数学模型时一般选择单一气相^{[4, 27]}，但实际上，由于页岩气藏渗流机理复杂，地层水可能参与流动而出现气水两相渗流。为此，笔者利用等效孔隙介质模型，同时考虑滑脱和应力敏感效应建立了页岩气藏气水两相渗流模型。主要假设条件包括： 三维气、水两相等温渗流，考虑气体滑脱及应力敏感效应，岩石为可压缩孔隙介质，流体为可压流体，考虑重力和毛细管压力。

对于三维气、水两相渗流，水相和气相的流动过程符合守恒定律，可以分别表示为^[28]：

(4)

(5)

其中，式(4)为水相渗流方程，式(5)为气相渗流方程。除了式(4)和式(5)，还需要联立辅助方程和约束方程，以减少变量的数量。

考虑重力作用下水相的总压强为：

(6)

综合考虑重力和气水毛细管压力作用下气相的总压强为：

(7)

对于饱和孔隙介质，水、气两相的饱和度之和为1，则约束方程为：

(8)

式中：S为饱和度；K为地层渗透率，D；K_r为相对渗透率；μ为黏度，mPa·s；R_s为溶解气油比；z为标高，m；B_w,B_g分别为水和气的体积系数；φ为地层孔隙度；p_c为毛细管压力，Pa；γ为流体重度，N/m³；q为源汇项，从单位体积储层岩石中采出(取负号)或向单位体积储层岩石中注入(取正号)的地面条件下的流体量，m³/d；η为单位换算系数；Φ_w为水相的压强与位置水头产生的压强之和，Pa；Φ_g为气相的压强与位置水头产生的压强之和，Pa；下标w，o和 g分别表示水，油和气。

式(1)—式(8)构成了页岩气的渗流数学模型。其中，式(4)和式(5)中的渗透率(K)和孔隙度(φ)受气体滑脱和应力敏感效应的影响，计算时要将式(1)—式(3)耦合到式(4)和式(5)中。

3 页岩气藏渗流模型数值求解

笔者采用有限差分方式对页岩气藏渗流数学模型进行差分离散，差分求解过程见参考文献[29]。在笛卡尔坐标系中，三维气、水两相偏微分方程的差分离散结果为:

(9)

(10)

其中，式(9)为气相差分方程，式(10)为水相差分方程。式中：T_g和T_w分别为气和水的传导系数,m²/s；D为网格中心的深度，m；V_i,j,k为第i,j,k个网格的体积,m³；q_wi,j,k和q_gi,j,k分别为单位体积储层岩石中采出水和气的速度，s^-1(注：这两项是源汇项，若网格内没有注采井，这两项均为零)。

采用隐式求压力、显式求饱和度的方法(即IMPES方法)，合并式(9)和式(10)，得到一个只含有水相压力一个未知量的方程。通过有限差分方法将非线性偏微分方程转化成关于压力的线性方程组：

(11)

式中：A为某个网格压力方程的系数；E为三维体系中与某个网格相邻的6个网格的压力方程系数与单位时间内气藏综合压缩性引起孔隙体积变化率之和；R_si为原始溶解气油比；B为压力方程的右端项；C_t为综合压缩系数，Pa^-1;Δt为时间步差，s。

用逐次线松弛迭代技术(line successive over relaxation，LSOR)对方程组进行求解^[28]。与常规油藏数值模拟不同的是，笔者建立的页岩气藏渗流模型考虑了滑脱和应力敏感效应对孔隙度和渗透率的影响，每个迭代步必须更新2个效应引起的孔渗参数变化，求解过程见图 1。笔者应用Fortran 90 在Powerstation 4.0 的编译环境下编程开发了页岩气藏数值模拟器。

4 算例模拟 4.1 定产量开采页岩气算例

该算例的目的是研究定产量衰竭式开发方式引起的页岩气藏压力变化及其对渗透率的影响，应用了三维地质模型，模拟条件是模型外边界封闭，4口页岩气生产井置于模型的4个角上，每口井以1×10⁴ m³/d的等产量衰竭方式开采。该算例中主要的地质和流体物性参数如下：储层深度为2 537.46 m，原始气藏压力为32.96 MPa,网格节点在x、y和z方向上的划分为10×10×3，平面上网格尺寸dx=dy=304.8 m，地层水黏度为0.50 mPa·s，气藏条件下页岩气黏度为0.03 mPa·s，束缚水和原始含气饱和度分别为0.20和0.80。表 1给出了各模拟层的厚度、孔隙度和3个方向上的渗透率的差异，体现了地质模型的非均质特性。

图 2为利用页岩气藏数值模拟器得到的页岩气藏压力在xy平面上随着开发时间的变化情况。由图 2可以看出：随着开发时间增长，页岩气藏的压力逐步下降；开发时间为2，5和8年时，页岩气藏局部位置的最大压力分别为32.1，29.9和27.6 MPa，生产井所射开地层(仅射开了最上面的一层)位置的压力分别为29.9，27.7和25.4 MPa；在任一开发时刻，4口生产井所在位置的压力较低；压力等值线并不完全对称，这是由于平面x和y方向上地层的渗透率不同所致。

图 3为利用页岩气藏数值模拟器得到的页岩气藏开采5年后3个方向上的渗透率分布特征。由图 3可以看出，页岩气藏各个部位3个方向上的渗透率均低于原始渗透率，而且生产井附近储层的渗透率最低。分析认为，由于生产井附近储层压力降低，导致储层受到的有效地应力增大，使储层受到的压实作用增大，从而引起渗透率急剧降低。

4.2 定压开采页岩气算例

该算例的目的是在定压衰竭式开采条件下，研究滑脱和应力敏感效应对页岩气井产量的影响。模拟条件是模型外边界封闭，4口页岩气井置于模型的4个角上，采用20 MPa的井底流压开采页岩气。模型的其他参数与定产量开采页岩气算例相同。

笔者利用页岩气藏数值模拟器模拟了滑脱和应力敏感效应对页岩气藏定压开采5年后压力分布的影响，结果如图 4所示。由图 4可以看出：滑脱和应力敏感效应均对储层压力分布产生影响，但滑脱效应对储层压力分布的影响不大，这与前人的观点吻合^[30-31]；应力敏感效应对储层压力分布的影响较大，特别是对于生产井附近的储层，应力敏感效应导致其压力等值线密集分布。

图 5为利用页岩气藏数值模拟器得到的滑脱和应力敏感效应对4口井总日产气量的影响。由图 5可以看出：投产第1年，页岩气井产气量迅速下降，但滑脱和应力敏感效应对产气量的影响极小；生产500 d后，滑脱和应力敏感效应对产气量的影响开始凸显，滑脱效应使产气量升高，应力敏感效应导致产气量降低；生产到3 500 d时，不考虑和考虑滑脱与应力敏感效应条件下的产气量分别为4.08×10⁴和3.65×10⁴ m³/d，滑脱和应力敏感效应导致产气量降低10.5%；与滑脱效应相比，应力敏感效应对页岩气井产气量的影响更大。

图 6为滑脱和应力敏感效应对4口页岩气井总累计产气量的影响。

由图 6可以看出，同一时刻累计产气量由高到低的顺序为：不考虑应力敏感只考虑滑脱时的产气量>不考虑滑脱和应力敏感时的产气量>考虑滑脱和应力敏感时的产气量>考虑应力敏感不考虑滑脱时的产气量，应力敏感比滑脱效应对累计产气量的影响大得多。

为了研究页岩气井产气量的衰减规律，利用页岩气藏数值模拟器模拟得到的数据，对4口井的总产气量与生产时间的关系进行了分段回归，结果见图 7。由图 7可以看出：总产气量递减规律可以近似地分为2段，在开采初期的730 d内，总产气量与生产时间的关系呈近似二次多项式规律递减；2年之后，总产气量与生产时间的关系呈近似自然对数规律递减。

4.3 计算结果与实际生产数据的对比

为了验证页岩气藏数值模拟器的可靠性，应用其计算的产气量与美国Haynesville页岩气藏W3井^[32]的实际生产数据进行对比，分析产气量递减趋势及吻合程度。数值模拟中的平面网格尺寸为130 m×130 m，孔隙度变化系数为0.08 Pa^-1，3个方向上的平均渗透率变化系数为0.62 Pa^-1；其他主要参数见定产量开采页岩气算例和表 1，模拟开始气井的产量定为1.6×10⁴ m³/d，与该井实际产气量保持一致。该井生产850 d的模拟计算结果和实际生产数据的对比见图 8。由图 8可以看出，模拟结果与生产数据的衰减趋势吻合度较高，从而验证了页岩气藏数值模拟器的可靠性。

5 结 论

1) 在综合考虑滑脱和应力敏感效应对页岩气藏渗流影响的条件下，建立了三维两相渗流数学模型，对模型进行了数值求解并开发了页岩气藏数值模拟器。

2) 数值模拟结果表明：页岩气藏投产后储层3个方向上的渗透率均低于储层的初始渗透率，生产井附近储层的渗透率最低；应力敏感效应导致页岩气藏储层渗透率下降，造成页岩气井产能降低。

3) 与滑脱效应相比，应力敏感效应对产气量的影响更大，不考虑应力敏感只考虑滑脱情况下的产气量最高，考虑应力敏感不考虑滑脱情况下的产气量最低。

4) 模拟结果表明，产气量递减速度呈现先快后慢的特征，与实际生产数据一致。产气量曲线可以近似分为两段，开采初期产气量与生产时间的关系呈近似二次多项式规律递减，后期近似呈自然对数规律递减。

5) 页岩气藏渗流影响因素多，机理复杂，建议将来重点开展页岩气藏渗流影响因素的试验研究，特别是压裂对导流能力的影响，进一步完善页岩气藏数值模拟器的功能。