致密砂岩气藏是一种重要的非常规油气资源，在油气田开发过程中具有重要地位^[1]，主要分布于鄂尔多斯盆地、大牛地气田等区域^[2-4]。致密砂岩气藏由于储层致密（渗透率小于0.1 mD），无法实现商业开采价值，因此，通常采用直井或水平井大规模压裂或二次重复压裂实现经济开发，从而提高致密砂岩气藏采收率^[5-6]。致密气储层大规模体积压裂可导致复杂裂缝网络和非均质压裂改造区域的形成，大量实验表明流体在致密储层中的流动存在明显的非线性渗流特征^[7]。为了更加准确地对致密砂岩气藏进行参数评价，部分学者建立了对应的数学模型进行参数反演和曲线分析。

井底压力动态模型求解方法包括（半）解析解法和数值解法。对于气藏而言，压裂井复杂裂缝（半）解析模型大多是线性或近似线性化求得压力动态分析模型^[8-11]，这种处理方法忽略了关井或生产过程中天然气的物性变化，对于拟启动压力梯度影响的（半）解析求解方法只有近似处理才可求解。基于以上（半）解析模型所存在的不足，相关学者开展了数值试井模型研究，对于复杂裂缝模型数值模型的研究采用嵌入式离散裂缝进行模拟^[12-18]，忽略了启动压力梯度的影响。Li等^[19]建立了考虑启动压力梯度影响的直井数值试井模型，分析了拟启动压力系数对双对数曲线的影响，但该模型未考虑复杂缝网、多井干扰和储层非均质性等的影响。

由于致密气藏多采用水力压裂技术，为了提高致密气采收率，其井距通常较小，从而会产生多井干扰。为了分析井间干扰对井底压力动态特征曲线的影响，闫正和等^[20]基于数值试井对海上干扰试井进行模拟设计及方案优化，部分学者基于（半）解析方法对压裂复杂裂缝多井干扰井底压力动态进行分析与研究，但上述模型忽略了储层非均质性和流体非线性渗流特征的影响^[21-23]。

基于上述文献调研结果，针对致密气非线性渗流特征，建立了综合考虑拟启动压力梯度和储层非均质性影响的压裂井复杂裂缝数值试井数学模型。该模型可修正非线性渗流体条件下传导率的计算公式，利用PEBI网格剖分方法对研究区进行网格剖分，采用有限体积和差分相结合的方法对模型进行求解，并通过简化模型与Saphir数值解对比验证模型的准确性，同时分析储层的非均质性、裂缝属性等参数对井底压力动态特征曲线的影响，以期为致密砂岩气藏大规模体积压裂提供一定理论基础。

致密储层非线性渗流特征导致致密砂岩气藏储层流体渗流不再满足达西渗流规律。对于致密气藏而言，是否启动压力梯度存在争议，姜瑞忠等^[24-26]提出了非线性渗流微分方程可准确描述非达西渗流与达西渗流的过程。根据Li等^[27]的研究成果，将反映启动压力梯度和非线性渗流的特征用一个参数来表示模拟非达西渗流与达西渗流过渡渗流的规律。表达式为

由于离散裂缝和压裂裂缝渗透率高，依然采用线性达西渗流公式进行描述

式（1）—式（2）中：k_m为基质系统渗透率，mD；k_F为裂缝系统渗透率，mD；p_m为基质系统压力，0.1MPa；P_F为裂缝系统压力，0.1MPa；μ为气体黏度，mPa·s；v_m为基质流体渗流速度，cm/s；v_F为裂缝流体渗流速度，cm/s；c为非线性系数，0.1MPa/cm。

将式（1）代入不考虑源汇项的基质系统渗流微分方程

将式（2）代入裂缝系统渗流微分方程

考虑井储和表皮影响时，内边界条件为

其中：$ W I=\frac{\theta k_{\mathrm{m}} h}{\mu \ln \left(d / r_{\mathrm{w}}\right)}$（井在基质系统）；

外边界条件包括封闭边界和定压边界2种

式（3）—式（7）中：ρ为气体密度，kg/m³；φ_m为基质系统孔隙度，%；φ_f为裂缝系统孔隙度，%；WI为井产能指数，m³/（s·0.1MPa）；p_wf为井底压力，0.1MPa；p_wb为井筒所在网格压力，0.1MPa；C为井储系数，m³/0.1MPa；n为外边界法向量；h为储层厚度，m；w_F为裂缝宽度，m；r_w为井筒半径，m；d为两网格中心距离，m；V_F为网格体积，m³；B为体积系数，m³/m³。q_mf和q_fm为基质和裂缝的交换流量，kg/s；

对基质和裂缝系统渗流微分方程[式（3）—式（5）]在空间上采用有限体积法、时间上采用差分法进行离散，基质、裂缝系统及内边界的离散微分方程分别如式（8）—式（10）。

其中：$ {{v^{n + 1}} = - \frac{k}{\mu }{\mathop{\rm grad}\nolimits} \left( {{p^{n + 1}}} \right);{\mathop{\rm grad}\nolimits} \left( {{p^{n + 1}}} \right) = p_i^{n + 1} - p_j^{n + 1}}$; $ {{\mathop{\rm div}\nolimits} {{(\rho v)}^{n + 1}} = - \frac{1}{{\mu B}}{T_{ij}}{\mathop{\rm grad}\nolimits} \left( {{p^{n + 1}}} \right)}$; $ {q_{{\rm{mf}}}^{n + 1} = {{\left( {\frac{1}{{\mu B}}} \right)}^{n + 1}}\left( {p_{\rm{m}}^{n + 1} - p_{\rm{f}}^{n + 1}} \right);q_{{\rm{fm}}}^{n + 1} = {{\left( {\frac{1}{{\mu B}}} \right)}^{n + 1}}\left( {p_{\rm{f}}^{n + 1} - p_{\rm{m}}^{n + 1}} \right)}$。

式（8）—式（10）中：T_ij为传导率；grad（）为压力梯度函数；div（）为梯度算子函数；上标n代表不同时间步。

将本文简化模型（不考虑拟启动压力梯度）与商业软件进行对比，绘制不同非线性系数下的井底压力动态特征曲线。图 1中，区域1的宽度为140 m，渗透率为1 mD，区域2和区域3的宽度均为180 m，渗透率为5 mD。生产井裂缝导流能力为2 000 mD·m，裂缝半长为100 m；注入井裂缝导流能力为2 000 mD·m，裂缝半长为150 m。由图 2可看出，本文方法计算得到的井底压力与Saphir数值解计算得到的井底压力拟合度较高，井底压力双对数曲线总体上分为6个流动阶段：井储阶段（阶段1）、表皮影响阶段（阶段2）、双线性流阶段（阶段3）、线性流阶段（阶段4）、井间干扰阶段（阶段5）及边界控制流阶段（阶段6）。在网格剖分方面商业数值模拟软件计算速度优于Matlab油藏模拟器软件（MRST），实际计算模拟方面MRST计算速度优于商业数值模拟软件。随着拟启动压力系数c的增大，流体渗流阻力增大，压力向外传播的速度减慢，压力及导数曲线上翘幅度变大，生产井和注入井产生的压力波优先沿离散裂缝传播（图 1c—1e）。

图 1

Fig. 1

下载原图 图 1 PEBI网格剖分及不同非线性系数压力分布 Fig. 1 PEBI mesh generation and pressure distribution of different nonlinearity coefficient

图 2

Fig. 2

下载原图 图 2 井底压力动态特征曲线 Fig. 2 Dynamic characteristic curve of wellbore pressure

诱导缝是指与压裂缝直接沟通的具有一定导流能力的裂缝，离散裂缝是指不直接和压裂缝沟通的具有一定导流能力的裂缝。储层中2口压裂直井分别位于不同的区域，其中区域1的渗透率高于区域2和区域3，2口井通过诱导缝和离散裂缝沟通。基于该物理模型和网格剖分结果计算测试井和邻井的井底压力，利用基本参数（表 1）和物理模型（图 3）研究邻井参数、裂缝参数和非均质储层参数对测试井双对数曲线的影响。

表 1(Table 1)

表 1 模拟基本参数 Table 1 Basic Simulation Parameters 参数 取值 参数 取值 原始地层压力/MPa 40 压裂裂缝导流能力/（mD·m） 2 500 储层温度/℃ 30 压裂裂缝半长/m 100 区域1渗透率/mD 2.0 诱导缝导流能力/（mD·m） 256 区域2渗透率/mD 0.2 诱导缝平均半长/m 90 区域3渗透率/mD 0.4 离散裂缝导流能力/（mD·m） 156 区域1宽度/m 240 离散裂缝平均半长/m 200 区域2宽度/m 130 井储系数/（m3·MPa-1） 0.3 区域3宽度/m 130 表皮系数 0 储层长度/m 1 000 非线性系数（/MPa·m-1） 0.3

下载CSV 表 1 模拟基本参数 Table 1 Basic Simulation Parameters

图 3

Fig. 3

下载原图 图 3 非均质储层压裂复杂裂缝井多井干扰裂缝几何模型 Fig. 3 Geometric model of fractured complex fractures in heterogeneous reservoirs with multiple well interference

基于图 3a物理模型示意图，井1为测试井，井2为邻井，多井干扰情况下绘制测试井井底压降及压降导数双对数曲线，分析邻井产量及邻井井底压力对测试井井底压降及压降导数曲线的影响。由图 4和图 5可看出，邻井生产方式对测试井中后期流动阶段井底压降及导数曲线的影响较大，邻井不同的生产方式主要影响测试井压降及压降导数中后期流动阶段的曲线形态。邻井生产方式包括定压生产和定产量生产，当测试井定产气量生产、邻井定压生产时，由于测试井与邻井通过诱导缝和离散裂缝沟通，邻井产生的压力波沿诱导缝和压裂缝传播到测试井，压降导数曲线出现不同程度的先上翘后变缓，邻井井底压降越高，压降导数曲线中后期流动阶段越低（图 4a），压力场分布如图 4b—4d所示。当测试井和邻井均定产气量生产时，邻井产量越大，晚期阶段压降及导数曲线上翘开始时间越早（图 5a），压力场分布如图 5b—5d所示。

图 4

Fig. 4

下载原图 图 4 邻井井底压力对测试井双对数曲线和压力场的影响 Fig. 4 The influence of adjacent well wellbore pressure on the wellbore pressure and pressure distribution

图 5

Fig. 5

下载原图 图 5 测试井与邻井产量比对测试井双对数曲线和压力场的影响 Fig. 5 The influence of adjacent well and test well rate ratio on the wellbore pressure and pressure distribution

诱导缝是指与压裂缝沟通且裂缝导流能力低于压裂缝导流能力的裂缝。由图 5可看出，诱导缝导流能力可影响双线性流阶段井底压降及导数曲线形态，诱导缝导流能力越大，诱导缝内流体向压裂缝“补充”的能力则越强。双线性流阶段压降导数曲线出现明显的“下凹”，诱导缝导流能力越大，压力波在诱导缝的传播速度越快，压降导数曲线“下凹”程度越明显（图 6）。

图 6

Fig. 6

下载原图 图 6 诱导缝导流能力对测试井双对数曲线和压力场的影响 Fig. 6 The influence of the induced fractures conductivity on the wellbore pressure curve and pressure distribution

假定测试井定产量生产、邻井定压力生产，分析测试井井底压降及导数曲线形态，由图 7可看出，诱导缝和离散裂缝条数可影响井底压降及导数中后期曲线形态。诱导缝和离散裂缝条数的增加与减少对观测井井底压降曲线的影响较小。离散裂缝条数分布如图 7所示，比较情形1和情形2可知，情形2未考虑离散裂缝，但2口井通过诱导缝沟通，离散裂缝的缺失对井底压降曲线影响较小。情形2和情形3相比较，情形3未考虑离散裂缝且缺少1条沟通2口井的诱导缝，曲线分析结果表明，诱导缝的存在对晚期阶段压力曲线的影响较明显。当诱导缝直接沟通2口井压裂裂缝时，观测井中后期阶段压力导数曲线上翘开始时间较早，反之亦然（图 7）。

图 7

Fig. 7

下载原图 图 7 离散裂缝条数对测试井双对数曲线和压力场分布的影响 Fig. 7 The influence of the induced fractures number on the wellbore pressure curve and pressure distribution

针对多井压裂形成的非均质区域，根据分区渗透率差异影响，分3种情况讨论。情形1：K₁ = 2 mD，K₂ = 0.5 mD，K₃ = 0.5 mD；情形2：K₁ = 2 mD，K₂ = 1 mD，K₃ = 1 mD；情形3：K₁ = 2 mD，K₂ = 2 mD，K₃ = 2 mD。测试井定产气量生产，邻井定井底压降生产时，区域渗透率差对中后期井底压降导数曲线有明显的影响（图 8）。当测试井所在区域渗透率不变时，尽管测试井与邻井通过诱导缝沟通，但当邻井所在区域渗透率越大时，邻井产生的压力波可更快的传播到储层和测试井，导致中后期阶段压力导数曲线上翘开始的时间越早，晚期阶段压力导数曲线下掉开始的时间越早（图 8a—8d）。

图 8

Fig. 8

下载原图 图 8 多区渗透率对测试井双对数曲线和压力场分布的影响 Fig. 8 The influence of the region permeability on the wellbore pressure curve and pressure distribution

假设区域1宽度分别为200 m，240 m和280 m，区域2和区域3的宽度相等，3个区域的渗透率分别为K₁ = 2 mD，K₂ = 1 mD，K₃ = 1 mD。区域1的宽度对井底压降动态特征曲线无明显影响，主要是由于诱导缝直接沟通了2口井。区域1宽度越大，邻井产生的压力波会很快传播到区域1并通过诱导缝传播到测试井压裂裂缝，因此，定压生产井干扰阶段压力导数曲线越高（图 9a—9d）。

图 9

Fig. 9

下载原图 图 9 区域1宽度对对测试井双对数曲线和压力场分布的影响 Fig. 9 The influence of the region1 width on the wellbore pressure curve and pressure distribution

（1）针对非均质致密气非线性渗流特征，以基于嵌入式离散裂缝模型为基础，建立了考虑拟启动压力梯度影响的压裂井复杂裂缝非均质储层多井干扰数值试井数学模型。

（2）将本文简化模型与商业软件模拟结果进行了对比，并根据压力及压降导数曲线验证了模型的准确性。

（3）当邻井定压生产时，测试井压降导数曲线后期表现为先上翘后变缓或下掉特征；当邻井定产量生产时，测试井井底压降曲线上翘反映开始的时间越早，拟启动压力梯度系数越大，曲线上翘幅度越大。