对均质有界储集层中1口有限导流垂直压裂气井, 物理模型假设条件为：①储集层具有应力敏感效应, 考虑渗透率应力敏感, 忽略孔隙度应力敏感; ②致密气为干气或者低含水饱和度致密气, 水处于束缚状态, 不存在启动压力梯度效应, 气体在储集层中的流动为单相渗流, 且满足达西定律; ③水力压裂后形成1条有限导流裂缝, 裂缝具有空间变导流和时间变导流特性, 气体在裂缝中的流动为一维流动, 考虑压裂液对储集层的伤害, 用裂缝表皮系数来衡量其对储集层伤害的大小; ④裂缝体积与井控体积相比非常小, 裂缝渗透率远大于储集层渗透率, 裂缝压力降低造成裂缝内气体体积膨胀对整个流动的影响非常小, 可忽略裂缝控制方程中拟压力对时间的导数项; ⑤考虑气体压缩系数及黏度随压力变化而变化, 气体压缩系数及偏差因子采用Dranchuk-Purvis-Robinson方法^[17]计算, 气体黏度采用Lee方法^[18]计算。⑥不考虑井筒储集效应、温度变化等其他因素对流动的影响。

据文献[19], 储集层控制方程为：

$\frac{\partial }{\partial x}\left[ {{K}_{\text{r}}}\left( p \right)\frac{\partial \psi }{\partial x} \right]+\frac{\partial }{\partial y}\left[ {{K}_{\text{r}}}\left( p \right)\frac{\partial \psi }{\partial y} \right]-$${{\alpha }_{3}}T{{q}_{f}}={{\alpha }_{1}}\phi {{C}_{t}}\left( p \right)\mu \left( p \right)\frac{\partial \psi }{\partial t}$ （1）

其中 $\psi \left( p \right)=2\int_{0}^{p}{\frac{p}{\mu Z}}\text{d}p$

裂缝控制方程为：

$\frac{\partial }{\partial l}\left[ {{K}_{\text{f}}}\left( l, t \right)\frac{\partial \psi }{\partial l} \right]+{{\alpha }_{3}}T{{q}_{f}}=0$ （2）

初始条件为：

$\psi \left( x, y, 0 \right)={{\psi }_{\text{i}}}$ （3）

内边界条件为：

已知井口产量时,

${{\left. \frac{\partial \psi }{\partial l} \right|}_{{{\Gamma }_{\text{in}}}}}=\frac{{{\alpha }_{2}}{{q}_{\text{sc}}}\left( t \right)T}{{{w}_{f}}{{K}_{f}}\left( 0, t \right)h}$ （4）

已知井底压力时,

${{\left. \psi \right|}_{{{\Gamma }_{\text{in}}}}}={{\psi }_{\text{w}}}\left( t \right)$ （5）

封闭外边界为：

${{\left. \frac{\partial \psi }{\partial N} \right|}_{{{\Gamma }_{\text{out}}}}}=0$ （6）

物质平衡方程为：

$\frac{{\bar{p}}}{Z}=\frac{{{p}_{\text{i}}}}{{{Z}_{\text{i}}}}\left( 1-\frac{{{G}_{\text{p}}}}{G} \right)$ （7）

1.3.1 空间变导流

裂缝导流能力随裂缝长度呈线性、指数或对数关系变化, 计算公式分别为^[14]：

${{F}_{\text{cl}}}={{F}_{\text{c0}}}\left( 1-{{a}_{s}}{{l}_{\text{D}}} \right)$ （8）

${{F}_{\text{cl}}}={{F}_{\text{c0}}}{{\text{e}}^{-{{b}_{s}}{{l}_{\text{D}}}}}$ （9）

${{F}_{\text{cl}}}={{F}_{\text{c0}}}\left[ 1-{{c}_{\operatorname{s}}}\ln \left( 1+{{l}_{\text{D}}} \right) \right]$ （10）

1.3.2 时间变导流

裂缝导流能力随生产时间呈对数或指数关系变化, 计算公式分别为^{[10, 11]}:

${{F}_{\text{ct}}}={{F}_{\text{ci}}}\left[ 1-\chi \ln \left( 1+t \right) \right]$ （11）

${{F}_{\text{ct}}}={{F}_{\text{ci}}}{{\text{e}}^{-\gamma t}}+{{F}_{\text{cr}}}$ （12）

对目前常用的指数变化关系（12）式进行修改：

${{F}_{\text{ct}}}={{F}_{\text{ci}}}\left( \eta {{\text{e}}^{-\frac{t}{8\ 760C}}}+1-\eta \right)$ （13）

系数η 控制导流能力衰减幅度, 系数C控制导流能力的衰减速度, 修改后系数的物理意义更明确, 有助于生产数据解释及分析。

1.3.3 双重变导流

同时考虑空间变导流和时间变导流, 即裂缝导流能力不仅随裂缝位置而变, 同时还随着时间而变。以空间变导流和时间变导流均采用指数关系变化为例, 推导出双重变导流计算公式：

${{F}_{\text{c}}}\left( l, t \right)={{F}_{\text{ci0}}}{{\text{e}}^{-{{b}_{s}}{{l}_{\text{D}}}}}\left( \eta {{\text{e}}^{-\frac{t}{8\ 760C}}}+1-\eta \right)$ （14）

以此类推可以获得不同函数组合情况下的双重变导流公式。

致密气储集层可能会存在应力敏感、启动压力梯度及滑脱效应等非线性渗流特征。应力敏感效应主要表现为储集层渗透率不再为常数, 而是随有效应力的增大而减小。大量实验^{[20, 21, 22]}表明致密气储集层具有较强的应力敏感效应, 对于分析长时间生产数据的产量递减分析来说必须予以考虑。气藏中存在启动压力梯度主要是由于气水作用造成^{[23, 24]}, 这与油藏中的启动压力梯度有本质区别, 因此气藏中存在一定程度的水是产生启动压力梯度的必要条件。此外, 致密气储集层可能还存在滑脱效应, 滑脱效应会增加视渗透率, 储集层越致密, 储集层压力越低, 滑脱效应越强, 目前考虑滑脱效应的渗流模型主要应用于埋深较浅的煤层气中^{[25, 26]}, 而对于储集层压力普遍较高的致密砂岩气藏, 滑脱效应的影响还有待考察。

1.4.1 滑脱效应

滑脱效应的大小主要由滑脱因子与储集层平均压力决定, 可由Klinkenberg公式表征^[27]：

${{K}_{\text{g}}}={{K}_{\infty }}\left( 1+\frac{b}{{{p}_{\text{m}}}} \right)$ （15）

为了考察致密砂岩气藏中滑脱效应对渗流影响的大小, 采用苏里格气田岩心做了96次滑脱效应实验, 结果见图1, 经双对数线性回归得滑脱因子与绝对渗透率之间的关系式为：

$b=0.009\ 11K_{\infty }^{^{-0.606\ \ 1}}$ （16）

图1

	Figure Option ViewDownloadNew Window
	图1 滑脱效应实验结果

假设致密气储集层绝对渗透率为0.1× 10^-3 μ m², 原始储集层压力30 MPa, 生产后期的储集层平均压力为10 MPa, 根据（15）式与（16）式可计算得到滑脱效应对早期、后期储集层视渗透率的影响分别为0.12%和0.37%。即使储集层绝对渗透率为0.01× 10^-3 μ m², 对早期、后期储集层视渗透率的影响也仅有0.49%和1.48%。计算该储集层条件下滑脱效应对Blasingame曲线的影响（见图2）, 对比可知, 考虑和不考虑滑脱效应情况下, 两者曲线基本重合, 说明其对产量递减曲线的影响非常小, 可以在致密气井产量递减分析过程中将其忽略。

图2

	Figure Option ViewDownloadNew Window
	图2 滑脱效应对Blasingame曲线的影响

1.4.2 应力敏感效应

储集层应力敏感的评价方法有很多种^[28], 目前多采用指数和幂律关系式。鉴于实验测试数据幂律关系式回归的相关度更高, 本文采用渗透率随有效应力呈幂律关系变化的公式^[29]：

${{K}_{\text{r}}}\left( p \right)={{K}_{\text{ri}}}{{\left( \frac{{{\sigma }_{\text{e}}}}{{{\sigma }_{\text{ei}}}} \right)}^{-{{S}_{\text{p}}}}}={{K}_{\text{ri}}}{{\left( \frac{\sigma -p}{\sigma -{{p}_{\text{i}}}} \right)}^{-{{S}_{\text{p}}}}}$ （17）

采用混合有限元方法对模型进行求解^{[30, 31]}：

$\iint_{\Omega }{{{F}_{\text{eq}}}}\text{d}\Omega =\iint_{{{\Omega }_{\text{m}}}}{{{F}_{\text{eq}}}}\text{d}{{\Omega }_{\text{m}}}+{{w}_{\text{f}}}\int_{{{\overline{\Omega }}_{\text{f}}}}{{{F}_{\text{eq}}}}\text{d}{{\overline{\Omega }}_{\text{f}}}$ （18）

将整个计算区域划分为两个部分, 一个是二维流动的储集层区域, 一个是一维流动的裂缝区域。

利用Galerkin加权余量法^[32]分别离散储集层和裂缝的控制方程, 得到储集层区域二维有限元方程为：

$\begin{align} & A{{K}_{\text{r}}}\left( p \right)\left[ {{b}_{i}}{{b}_{i}}+{{c}_{i}}{{c}_{i}}+\frac{{{\alpha }_{1}}\phi \mu \left( p \right){{C}_{\text{t}}}\left( p \right)}{6\Delta t{{K}_{\text{r}}}\left( p \right)} \right]\psi _{i}^{n+1}+ \\ & A{{K}_{\text{r}}}\left( p \right)\left[ {{b}_{i}}{{b}_{j}}+{{c}_{i}}{{c}_{j}}+\frac{{{\alpha }_{1}}\phi \mu \left( p \right){{C}_{\text{t}}}\left( p \right)}{12\Delta t{{K}_{\text{r}}}\left( p \right)} \right]\psi _{j}^{n+1}+ \\ & A{{K}_{\text{r}}}\left( p \right)\left[ {{b}_{i}}{{b}_{k}}+{{c}_{i}}{{c}_{k}}+\frac{{{\alpha }_{1}}\phi \mu \left( p \right){{C}_{\text{t}}}\left( p \right)}{12\Delta t{{K}_{\text{r}}}\left( p \right)} \right]\psi _{k}^{n+1}- \\ & \frac{{{K}_{\text{r}}}\left( p \right)L}{3}\frac{\partial \psi _{i}^{n+1}}{\partial N}-\frac{{{K}_{\text{r}}}\left( p \right)L}{6}\frac{\partial \psi _{j(k)}^{n+1}}{\partial N}= \\ & \frac{{{\alpha }_{1}}A\phi \mu \left( p \right){{C}_{\text{t}}}\left( p \right)}{6\Delta t}\psi _{i}^{n}+\frac{{{\alpha }_{1}}A\phi \mu \left( p \right){{C}_{\text{t}}}\left( p \right)}{12\Delta t}\psi _{j}^{n}+ \\ \end{align}$ $\frac{{{\alpha }_{1}}A\phi \mu \left( p \right){{C}_{\text{t}}}\left( p \right)}{12\Delta t}\psi _{k}^{n}$ （19）

裂缝区域一维有限元方程为：

$\frac{{{F}_{\text{c}}}\left( l, t \right)}{L}\psi _{i}^{n+1}-\frac{{{F}_{\text{c}}}\left( l, t \right)}{L}\psi _{j}^{n+1}+{{F}_{\text{c}}}\left( l, t \right)\frac{\partial {{\psi }_{i}}}{\partial l}=0$ （20）

由（19）式和（20）式建立储集层区域的有限元刚度矩阵和裂缝区域的刚度矩阵, 再将两者组合成系统刚度矩阵。根据窜流关系式, 通过储集层和裂缝单元叠加消除裂缝与储集层交界处的边界项, 即（19）式等号左端最后两项以及（20）式等号左端最后一项, 具体组合方法可参照文献[30, 31]。

裂缝导流能力、储集层渗透率、气体黏度以及压缩系数均会随着时间而变, 实际计算中每个时间步根据各参数的具体计算方法来计算。此外, 模型采用两种方法计算储集层平均压力：①根据生产数据迭代获取井控储量时采用（7）式计算储集层平均压力; ②全历史曲线拟合时采用网格单元压力加权平均的方法计算储集层平均压力。最后利用线性方程组求解器求解组合成的线性方程组, 从而获得模型的解。

假定储集层渗透率为0.1× 10^-3 μ m², 原始储集层压力为30 MPa, 储集层温度为100 ℃, 有效厚度为8 m, 孔隙度为10%, 圆形封闭半径为500 m, 天然气相对密度为0.7, 裂缝半长为100 m, 单独考虑裂缝空间变导流, 缝口导流能力为200× 10^-3 μ m²· m。为了便于对比分析, 可设置变导流系数a_s, b_s, c_s的值使得缝端导流能力均为20× 10^-3 μ m²· m。

通过计算得不同裂缝空间变导流条件下的Blasingame典型曲线（见图3）。对比计算结果可知, 裂缝空间变导流能力对产量递减典型曲线影响主要在早期, 变导流会降低早期典型曲线值, 类似于裂缝表皮系数的影响, 常规不考虑变导流的模型会造成表皮系数解释结果偏大。不同变导流条件下晚期典型曲线几乎重合, 尤其β 曲线完全重合可以说明空间变导流不影响压力的扩散速度, 拟稳态阶段的出现时间不受其影响。另外, 对比不同空间变导流关系曲线可以发现, 在相同的缝端和缝口导流能力情况下, 指数变化对典型曲线的降低幅度最大, 线性关系最小, 对数关系居中。

图3

	Figure Option ViewDownloadNew Window
	图3 空间变导流对Blasingame曲线的影响

采用相同的储集层参数, 单独考虑裂缝时间变导流的影响。假定裂缝初始导流能力为200× 10^-3 μ m²· m, 生产1年后导流能力值衰减至20× 10^-3 μ m²· m, 随后基本稳定。分别根据对数和指数变化关系计算时间变导流能力对Blasingame典型曲线的影响（见图4）。分析计算结果可知, 相比于标准定导流典型曲线, 裂缝时间变导流能力对产量递减曲线的影响比较明显, 其主要影响在中期, 表现为产量曲线及产量积分曲线下掉, β 曲线上升, 形成S型。与空间变导流的影响类似, 时间变导流对典型曲线晚期的影响非常小, 从晚期β 曲线重合可知时间变导流同样不会影响储集层的压力扩散速度。此外, 对比对数关系与指数关系条件下的典型曲线可知, 对数关系下的典型曲线光滑度不如指数关系, 另外对数关系对典型曲线的影响幅度要小于指数关系, 这是因为对数关系整个过程的变化率比较平稳, 而指数关系的变化率由大变小。

图4

	Figure Option ViewDownloadNew Window
	图4 时间变导流对Blasingame曲线的影响

以指数变导流关系为例来分析变导流系数对典型曲线的影响。图5为设定变导流时间系数为0.1 a, 不同变导流衰减系数条件下典型曲线对比。图6为设定变导流衰减系数为0.9, 不同变导流时间系数条件下的典型曲线对比。由计算结果可知, 裂缝变导流衰减系数越大, 产量曲线及产量积分曲线下掉的幅度越大, S型曲线越明显。裂缝变导流时间系数越大, 即裂缝导流能力衰减速度越小, 产量曲线及产量积分曲线下掉的开始时间越晚, 其对典型曲线的影响也越小, 当变导流时间系数很大时, 其影响很难辨别, 基本可以忽略。由此可知时间变导流模型中裂缝衰减系数控制典型曲线下掉的幅度大小, 而时间系数控制着典型曲线下掉的时间早晚。这一特点对于时间变导流裂缝井的产量递减曲线拟合具有重要指导意义。

图5

	Figure Option ViewDownloadNew Window
	图5 裂缝变导流衰减系数对Blasingame曲线的影响（C=0.1 a）

图6

	Figure Option ViewDownloadNew Window
	图6 裂缝变导流时间系数对Blasingame曲线的影响（η =0.9）

图7为不同裂缝变导流条件下的Blasingame曲线对比, 其中空间变导流和时间变导流均采用指数变化关系, 变导流系数b_s值为2, η 值为0.9, C值为0.2 a。对比曲线可知, 典型曲线早期定导流曲线与时间变导流曲线重合, 而空间变导流曲线与双重变导流曲线重合, 这说明典型曲线早期主要受空间变导流的影响; 典型曲线中期空间变导流曲线与双重变导流相差较大, 而时间变导流曲线与双重变导流曲线差别较小, 这说明典型曲线中期受时间变导流和空间变导流双重影响, 但时间变导流起主导作用; 典型曲线晚期基本不受变导流的影响, 所有曲线基本重合。综合对空间变导流和时间变导流的单因素影响分析可知, 同时考虑空间变导流和时间变导流条件下的Blasingame曲线是两者效应的叠加。

图7

	Figure Option ViewDownloadNew Window
	图7 双重变导流对Blasingame曲线的影响

致密气储集层实际生产井中可能会存在裂缝双重变导流和应力敏感共同作用的情况, 为此在模型中同时考虑裂缝双重变导流和地层应力敏感效应。假定上覆岩石压力为60 MPa, 空间变导流和时间变导流均采用指数变化关系, 变导流系数b_s值为2, η 值为0.9, C值为0.2 a, 分别计算双重变导流作用下不同应力敏感系数的Blasingame典型曲线（见图8）。对比可知, 不同应力敏感系数条件下的曲线在早期阶段基本重合, 但一段时间后曲线差异逐渐加大。与时间变导流的影响类似, 应力敏感效应会造成中后期阶段产量递减典型曲线值减小, 应力敏感系数越大曲线值减小的幅度越大。与时间变导流效应不同的是, 应力敏感会使得地层渗透率减小, 由此减缓压力扩散的速度, 使得边界响应的时间推迟, 图中β 曲线后期能够清晰地反映出该特征, 应力敏感效应越强拟稳态开始的时间越晚, 这种特点可用于区分应力敏感和时间变导流。

图8

	Figure Option ViewDownloadNew Window
	图8 应力敏感与双重变导流共同作用下的Blasingame曲线

另外, 晚期进入拟稳态后, 变导流能力效应和应力敏感效应影响极小, 即：变导流能力效应和应力敏感效应不会对井控动态储量的拟合结果产生影响。