针对深层页岩气在生产过程中存在较为严重的气水两相渗流问题，建立了深层页岩气压裂水平井气水两相流动模型(图1)，并对模型作出如下假设：(1)矩形储层，厚度相等，上下均为封闭边界；(2)温度在整个生产过程中不发生变化；(3)将气体在页岩基质中的吸附解析以及气体的压缩性等非线性影响因素考虑到模型当中；(4)气、水两相相互独立，气体不溶于水；(5)储层流体为二维流动，裂缝为一维流动，忽略流体重力影响；(6)不考虑压裂液注入过程，即不考虑生产过程中由于压裂液和地层水造成的气水两相流动。

图  1  深层页岩气压裂水平井物理模型

Figure 1.  Physical model of a fractured horizontal well in the deep shale gas reservoir

储层基质系统是指未改造区的页岩基质，根据四川盆地某页岩气井页岩原始含水饱和度测定结果，基质系统含水饱和度高于束缚水饱和度，因此考虑基质系统气水两相流动。

基质孔隙压缩系数

气体压缩系数

超临界气体解吸压缩系数^[13]

基质综合压缩系数

式中，c_m为基质孔隙压缩系数；ϕ_m为基质区孔隙度，小数；p_m为基质区压力，MPa；c_g为气体压缩系数；c_d为超临界气体解吸压缩系数；p_sc为标准状况下压力，MPa；Z为压缩因子；T_sc为标准状况下温度，K；Z_sc为标准条件下压缩因子；T为地层温度，K；V_L为兰氏体积，m³/kg；p_L为兰氏压力，MPa；ρ_g为气相密度，kg/m³；ρ_a为吸附相密度，kg/m³；c_t为基质综合压缩系数。

根据质量守恒定律，考虑气体在页岩基质中的渗流、扩散、滑脱和超临界吸附作用，得到基质中气相渗流方程为

其中　　　　　　　${k}_{\text{a}}=k{{\rm{exp}}\left[{{a}}\left ({{{p}}}_{\text{mg}}-{{{p}}}_{\text{e}}\right)\right]}$

页岩基质表观渗透率为

式中，k_a为考虑了应力敏感性和表观渗透率的渗透率，10⁻³μm²；p_mg为基质区气相压力，MPa；µ_g为气相黏度，mPa·s；S_mg为基质区气相饱和度，小数；k为表观渗透率，10⁻³μm²；k_m0为页岩基质渗透率，10⁻³μm²；p_e为初始压力，MPa；p为压力，MPa；D_k为克努森扩散系数，mm²/s；a为应力敏感系数；c为稀薄系数；d为滑脱系数。

同理得到基质中水相渗流方程为

式中，k_mw为基质区水相相对渗透率；p_mw为基质区水相压力，MPa；S_mw为基质区水相饱和度，小数；µ_w为水相黏度，mPa·s；B_w为水相体积系数。

根据质量守恒定律，得到改造区中气相渗流方程为

式中，k_f为考虑应力敏感性后得到的渗透率，10⁻³μm²；k_fg为改造区气相相对渗透率；p_fg为改造区气相压力，MPa；ϕ_f为改造区孔隙度，小数；S_fg为改造区气相饱和度，小数；B_g为气相体积系数；b为应力敏感系数；k_f0为页岩改造区渗透率，10⁻³μm²。

同理可得到SRV中水相渗流方程为

式中，k_fw为改造区水相相对渗透率；p_fw为改造区水相压力，MPa；S_fw为改造区水相饱和度，小数。

对于人工主裂缝，由于其改造性好，渗透率大，考虑气体高速运动产生的非达西现象，根据Forchheimer方程，可以得到非达西校正系数为

其中　　　　$\beta = \dfrac{{6.92 \times {{10}^{10}}}}{{{{\left ( {{k_{\text{F}}}{k_{{\text{Fg}}}}} \right)}^{0.5}}{{\left[ {{\phi _{\text{F}}}\left ( {1 - {S_{{\text{F}}w}}} \right)} \right]}^{1.5}}}}$

式中，$\varepsilon _{\text{n}} $为气相非达西校正系数；k_F为主裂缝区渗透率，10⁻³μm²；k_Fg为主裂缝区气相相对渗透率；$\beta $为非达西系数；v为达西流速，cm/s；ϕ_F为主裂缝区孔隙度，小数；S_Fw为主裂缝中水相饱和度，小数。

忽略主裂缝区应力敏感效应及毛管力，即气相压力与水相压力相等，得到主裂缝中气相渗流方程为

式中，p_F为主裂缝中压力，MPa；Q_Fg为地面标准状况下产气量，m³/s；S_Fg为主裂缝区气相饱和度，小数。

同理可得到主裂缝中水相渗流方程为

式中，Q_Fw为地面标准状况下的产水量，m³/s；k_Fw为主裂缝区水相相对渗透率。

饱和度约束方程

由于主裂缝渗透性好，忽略主裂缝中的毛管力，毛管压力方程为

式中，S为饱和度，小数；p为压力，MPa；下标m代表基质区，下标f代表改造区，下标F代表主裂缝区。

对三区渗流方程进行有限元离散，得到气水流动控制方程。

基质系统气水流动控制方程

改造区系统气水流动控制方程

主裂缝系统气水流动控制方程

式中，S_g0为初始含气饱和度，小数；c_w为地层水压缩系数；η为导压系数；S_me为基质有效含水饱和度，小数；p_mw为基质区水相压力，MPa；p_mc为基质区毛管力，MPa；p_fc为改造区毛管力，MPa；R为摩尔气体常数，J/(mol · K)。

先将整体空间划分为基质与主裂缝两个空间，运用伽辽金法，得到基质系统的等效积分弱形式为

网格划分选择的是三角形网格，因此其基函数为

式中，A为模型单元面积。

主裂缝系统的等效积分弱形式为

其中　　　　　$M_{\mathrm{p}}^{\prime}=\dfrac{S_{\mathrm{g0}}}{p_{\mathrm{F}}} C_{\mathrm{wF}}$

用COMSOL Multiphysics对深层页岩气水两相流动模型进行求解。通过COMSOL构建页岩压裂水平井模型并进行网格划分，用COMSOL中的系数形式偏微分方程(PDE)模块将气水两相控制方程及等效积分“弱”形式输入到模型中，对深层页岩气水两相流动模型求解，得到压力及含水饱和度变化场图。

采用三角形网格进行划分，为了使模拟结果更精确，将主裂缝周围的网格进行加密处理，模型网格包含752个域单元和99个边界元(图2)。

图  2  网格剖分结果

Figure 2.  Meshing of the model

模型基础参数取值参考四川盆地某页岩气井，原始地层压力为50 MPa，原始地层温度为403.15 K，初始含水饱和度为30%，兰氏体积为3 m³，兰氏压力为6.5 MPa，基质孔隙度为4%，基质渗透率为1×10⁻⁸ μm²，SRV区孔隙度为2%，SRV区渗透率为2×10⁻³ μm²，主裂缝半长150 m，主裂缝宽度1 mm，主裂缝渗透率为200×10⁻³ μm²。

模拟运行后，不同时刻地层压力变化如图3所示，由于主裂缝渗透率较大，因此随着生产进行，主裂缝处压力下降最快，并逐渐向外传播，由改造区传播至基质。不同时刻的含水饱和度变化如图4所示，主裂缝处流体最先流向井筒，随着生产进行，流体流向主裂缝并最终汇入井筒，动用区域也从主裂缝开始向外逐步扩大。

图  3  地层压力随累计生产时间变化场图

Figure 3.  Contour plot of variation of formation pressure with cumulative production time

图  4  含水饱和度随累计生产时间变化场图

Figure 4.  Contour plot of variation of water saturation with cumulative production time

深层页岩气水运动受到多种因素的影响，分析基质渗透率、改造区面积和主裂缝渗透率等因素对气水运动的影响程度，由于含水饱和度是最为直观的气水运动参数，因此引入改造区含水饱和度变化率和含水饱和度场评价各因素对气水运动的影响，同时选择累积产气量作定量评价。

式中，S_wm为变化前的含水饱和度，小数；S_wn为变化后的含水饱和度，小数；$S^\prime $为含水饱和度变化率，%。

对比基质渗透率为1×10⁻⁸、2×10⁻⁸、4×10⁻⁸、6×10⁻⁸ μm² 等4种情况下的含水饱和度变化率、累计产气量的变化，如图5、图6所示，可以看出，基质渗透率越大，越有利于流体流动，改造区边缘的含水饱和度越高，说明基质中的水补充进入改造区，从而使得含水饱和度变化率降低；更多的水进入改造区后，推动气体采出，使得产量增高，但增加幅度越来越小。

图  5  基质渗透率与含水饱和度变化率、累产量关系

Figure 5.  Correlations of water saturation variation rate and cumulative gas production with matrix permeability

图  6  改造区含水饱和度随基质渗透率变化场图

Figure 6.  Contour plot of variation of the reconstruction zone water saturation with matrix permeability

模拟不同改造区面积条件下的含水饱和度变化率、累计产气量的变化，分析改造面积对气水运动的影响程度，由图7、图8可以看出，改造面积越大，累计产量也越大，当改造面积增大到一定程度时，累计产气量增加幅度减小，即当改造面积增大到一定程度后，对产量的贡献降低。改造面积越大，改造区含水饱和度变化率越小，分析认为，增大改造面积会使得改造区边缘改造不充分，导致生产过程中有裂缝闭合，改造区边缘流体难以流入主裂缝，即随着改造面积增大，主裂缝含气饱和度变小，且由于改造面积增大，有更多的压裂液滞留，从而导致含水饱和度变化率降低。

图  7  改造区面积与含水饱和度变化率、累计产气量关系

Figure 7.  Correlations of water saturation variation rate and cumulative gas production with the reconstruction zone area

图  8  改造区含水饱和度随改造面积变化场图

Figure 8.  Contour plot of variation of water saturation with the reconstruction zone area

由于基质中的流体是先流入改造区，再进入主裂缝，然后从主裂缝流入水平井筒，因此主裂缝渗透率会直接影响产量。从图9可以看出，主裂缝渗透率越大，产气量越大，含水饱和度变化率越快，但增加幅度越来越小。这是因为主裂缝内的流体会先流入井筒，改造区中的流体对主裂缝进行补充，最后再由基质对改造区进行补充，因此前期主要受主裂缝渗透率影响，后期主要受基质渗透率影响。从图10可以看出，主裂缝渗透率越大，流体向主裂缝流动，使得裂缝之间的含水饱和度越小。

图  9  主裂缝渗透率与含水饱和度变化率、累计产气量关系

Figure 9.  Correlations of water saturation variation rate and cumulative gas production with the main fracture permeability

图  10  改造区含水饱和度随主裂缝渗透率变化场图

Figure 10.  Contour plot of variation of water saturation with the main fracture permeability in the reconstruction zone

(1)建立了考虑超临界吸附和非线性流动机理的深层页岩气压裂水平井气水两相流动数学模型，基于有限元方法，采用COMSOL Multiphysics软件求解得到累计产气量及含水饱和度变化场图，并进行气水运动规律分析。

(2)压裂设计需要进行综合考量，并不是改造面积越大越好，改造区面积越大，使得改造区边缘改造不充分，导致生产过程中有裂缝闭合，改造区边缘流体难以流入主裂缝。

(3)开采过程中，主裂缝内的流体会先流入井筒，改造区中的流体对主裂缝进行补充，最后再由基质对改造区进行补充，因此前期主要受主裂缝渗透率影响，后期主要受基质渗透率影响。