引言

页岩气是一种潜在的、资源量巨大的非常规天然气资源，具有开采技术要求高、开采寿命长、 稳产周期长的特点。近些年来，严峻的能源形势和能源价格的快速增长，使页岩气资源在全球受到了广泛的重视^[1-5]。 针对页岩气储层特有的解吸吸附特征，国内外许多学者通过修正物质平衡时间^[6-8]、建立半解析数学模型^[9]、整合Blasingame产能递减及修正双孔隙模型^[10]等方法在页岩气渗流机理方面取得了一系列研究成果。

针对页岩气渗流模式的压裂水平井^[11-15]的产能计算相当复杂，综合应用各级渗流模式需要建立复杂的数值模型，并需要提供使用渗流模式的判断依据，应用过程中有一定难度。为了便于简单分析，有些学者提出了一种三线性渗流模型^[11-12]，但模型分析较为复杂，更重要的是未考虑页岩气储层解吸-吸附机理，为了更合理地计算预测储层的生产能力，本文重新构建了页岩气藏三线性渗流模型，建立了考虑解吸吸附的数学模型，进行了产量公式和压力公式的推导，并应用此模型进行了某页岩储层压裂水平井的产能计算和分析，得到了与实际较为符合的结论。

1 三线性渗流模型

页岩气藏水平井压裂会形成带有一条主裂缝的裂缝网络，页岩气则通过基岩解吸、扩散和渗流到裂缝网络中，然后通过主裂缝流入井内，形成了完整的渗流过程。然而基岩、裂缝网络和主裂缝的渗流有着比较明显的区别。基岩是页岩气的存储空间，其渗流多级模式突出，气井开采初期，首先产出的是其中的游离气，之后由于气藏压力降低，基质表面的吸附气开始解吸形成新的游离气从而补充气源，解吸持续时间较长，所以气井具有开采初期产量较高、随后迅速递减且开采寿命长的特点；裂缝网络增加渗流通道，形成基岩与主裂缝的连通桥梁，同时也有部分游离气参与渗流；主裂缝是主要渗流通道，是井与储层之间的接触桥梁，其渗透率很高。以上3种介质中的渗流很复杂，按其关系进行简化。即分成如图 1中三部分线性渗流模式。

(1) 基岩线性渗流：基岩中页岩气通过裂缝面渗入裂缝网络简化为通过基岩整体以线性渗流模式流入裂缝网络系统，其几何规模为平面上的压裂网络区域的3/4；

(2) 裂缝网络线性渗流：把页岩气在裂缝网络体系中的渗流简化为由两主裂缝中间分别向各自主裂缝的线性渗流，渗流通道是裂缝网络等效渗透率，其值可等价为天然裂缝渗透率，一般为基质渗透率的100倍；

(3) 主裂缝线性渗流：将主裂缝内的渗流简化为接受两边裂缝网络窜流的主裂缝空间的线性渗流，其渗透率可达到几百$\sim$几千毫达西。

2 考虑解吸-吸附数学模型

对上述渗流模型的3个部分建立渗流数学模型。

2.1 基岩线性渗流

考虑基岩中的解吸-吸附作用，建立质量守恒方程为

$-\frac{\partial }{\partial x}\left( {{\rho }_{\text{gm}}}{{v}_{\text{m}}} \right)=\frac{\partial }{\partial t}({{\phi }_{\text{m}}}{{\rho }_{\text{gm}}}+{{\rho }_{\text{gsc}}}{{V}_{\text{E}}})$

(1)

其中：

$v_{\rm{m}} = - \dfrac{{K_{\rm{m}} }}{\mu }\dfrac{{\partial p_{\rm{m}} }}{{\partial x}}$

$\rho _{{\rm{gm}}} = \dfrac{{p_{\rm{m}} M}}{{z{\rm{R}}T}}$

$V_{\rm{E}} = \dfrac{{V_{\rm{L}} p_{\rm{m}} }}{{p_{\rm{L}} + p_{\rm{m}} }}$

为 Langmuir吸附气^[13]。

式中：

$\rho _{\rm{g}}$—地层状态气体密度，g/m³；

v—渗流速度，m/s；

t—生产时间，s；

$\phi$—孔隙度，%；

$\rho _{\rm{gsc}}$—标准状态气体密度，g/m³；

$V_{\rm{E}}$—页岩吸附气量，m³/t；

K—渗透率，mD ；

$\mu$—气体黏度，mPa$\cdot$s；

p—压力，MPa；

M—气体摩尔质量，g/mol；

z—气体偏差因子，无因次；

R—摩尔常数，J/(mol$\cdot$K)；

T—气藏温度，K；

$V_{\rm{L}}$—Langmuir体积，m³/t；

$p_{\rm{L}}$—Langmuir压力，MPa。

下标

m-基质系统。

式(1)左边整理为

$ - \dfrac{\partial }{{\partial x}}\left( {\rho _{{\rm{gm}}} v_{\rm{m}} } \right)= \dfrac{\partial }{{\partial x}}\left( {\dfrac{{K_{\rm{m}} }}{\mu }\dfrac{{p_{\rm{m}} M}}{{z{\rm{R}}T}}\dfrac{{\partial p_{\rm{m}} }}{{\partial x}}} \right)= \dfrac{M}{{{\rm{R}}T}}\dfrac{\partial }{{\partial x}}\left( {\dfrac{{K_{\rm{m}} }}{\mu }\dfrac{{p_{\rm{m}} }}{z}\dfrac{{\partial p_{\rm{m}} }}{{\partial x}}} \right)$

(2)

式(1)的右边两项整理为

第一项

$\dfrac{\partial }{{\partial t}}\left( {\phi _{\rm{m}} \dfrac{{p_{\rm{m}} M}}{{z{\rm{R}}T}}} \right) = \\ \phi _{\rm{m}} \dfrac{M}{{{\rm{R}}T}}\dfrac{\partial }{{\partial t}}\left( {\dfrac{{p_{\rm{m}} }}{z}} \right)= \\ \phi _{\rm{m}} \dfrac{M}{{{\rm{R}}T}}\left[{p_{\rm{m}} \dfrac{\partial }{{\partial t}}\left( {\dfrac{1}{z}} \right) + \dfrac{1}{z}\dfrac{{\partial p_{\rm{m}} }}{{\partial t}}} \right] = \\ \phi _{\rm{m}} \dfrac{M}{{{\rm{R}}T}}\left( { - \dfrac{{p_{\rm{m}} }}{{z^2 }}\dfrac{{\partial z}}{{\partial p_{\rm{m}} }}\dfrac{{\partial p_{\rm{m}} }}{{\partial t}} + \dfrac{1}{z}\dfrac{{\partial p_{\rm{m}} }}{{\partial t}}} \right) = \\ \phi _{\rm{m}} \dfrac{M}{{{\rm{R}}T}}\dfrac{{p_{\rm{m}} }}{z}\left( {\dfrac{1}{{p_{\rm{m}} }} - \dfrac{1}{z}\dfrac{{\partial z}}{{\partial p_{\rm{m}} }}} \right)\dfrac{{\partial p_{\rm{m}} }}{{\partial t}} = \\ \phi _{\rm{m}} \dfrac{M}{{{\rm{R}}T}}\dfrac{{p_{\rm{m}} }}{z}C_{{\rm{gm}}} \dfrac{{\partial p_{\rm{m}} }}{{\partial t}}$

(3)

其中： $ C_{{\rm{gm}}} = \dfrac{1}{{p_{\rm{m}} }} - \dfrac{1}{z}\dfrac{{\partial z}}{{\partial p_{\rm{m}} }} $

第二项

$\dfrac{\partial }{{\partial t}}\left( {\rho _{{\rm{gsc}}} \dfrac{{V_{\rm{L}} p_{\rm{m}} }}{{p_{\rm{L}} \!+\! p_{\rm{m}} }}} \right)\!=\!\rho _{{\rm{gsc}}} V_{\rm{L}} \dfrac{{(p_{\rm{m}}\!+\!p_{\rm{L}} )\dfrac{{\partial p_{\rm{m}} }}{{\partial t}}\!-\!p_{\rm{m}} \dfrac{{\partial p_{\rm{m}} }}{{\partial t}}}}{{(p_{\rm{m}}\!+\!p_{\rm{L}} )^2 }} \!=\!\rho _{{\rm{gsc}}} V_{\rm{L}} \dfrac{{p_{\rm{L}} }}{{(p_{\rm{m}} + p_{\rm{L}} )^2 }}\dfrac{{\partial p_{\rm{m}} }}{{\partial t}}$

(4)

式(3)+式(4)得

${\phi _{\rm{m}}}\dfrac{M}{{{\rm{R}}T}}\dfrac{{{p_{\rm{m}}}}}{z}{C_{{\rm{gm}}}}\dfrac{{\partial {p_{\rm{m}}}}}{{\partial t}} + {\rho _{{\rm{gsc}}}}{V_{\rm{L}}}\dfrac{{{p_{\rm{L}}}}}{{{{\left( {{p_{\rm{m}}} + {p_{\rm{L}}}} \right)}^2}}}\dfrac{{\partial {p_{\rm{m}}}}}{{\partial t}} = \\ {\phi _{\rm{m}}}\dfrac{M}{{{\rm{R}}T}}\dfrac{{{p_{\rm{m}}}}}{z}\left[{{C_{{\rm{gm}}}} + \dfrac{{{\rho _{{\rm{gsc}}}}{V_{\rm{L}}}{p_{\rm{L}}}}}{{{\phi _{\rm{m}}}\dfrac{M}{{{\rm{R}}T}}\dfrac{{{p_{\rm{m}}}}}{z}{{\left( {{p_{\rm{m}}}\!+\!{p_{\rm{L}}}} \right)}^2}}}} \right]\dfrac{{\partial {p_{\rm{m}}}}}{{\partial t}}\!$

(5)

定义考虑解吸-吸附的基质综合压缩系数^[14]为

$C_{\rm{tm}} = C_{\rm{gm}} + \dfrac{{\rho _{\rm{gsc}} V_{\rm{L}} p_{\rm{L}} }}{{\phi _{\rm{m}} \overline {\rho _{\rm{g}} } (p_{\rm{m}} + p_{\rm{L}} )^2 }}$

(6)

式中：

$C_{\rm{t}} $—压缩系数，MPa$^{-1}$；

$\overline {\rho _{\rm{g}} } $—地层气体平均密度，g/m³。

由此，基质渗流基本微分方程式(1)整理为

$\dfrac{\partial }{{\partial x}}\left( {\dfrac{{{K_{\rm{m}}}}}{\mu }\dfrac{{{p_{\rm{m}}}}}{z}\dfrac{{\partial {p_{\rm{m}}}}}{{\partial x}}} \right) = {\phi _{\rm{m}}}{C_{{\rm{tm}}}}\dfrac{{{p_{\rm{m}}}}}{z}\dfrac{{\partial {p_{\rm{m}}}}}{{\partial t}}$

(7)

定义拟压力函数及导压系数为

${m_{\rm{m}}} = 2\int_0^{{p_{\rm{m}}}} {\dfrac{p}{{\mu z}}{\rm{d}}p}$

(8)

${\eta _{\rm{m}}} = 2.637 \times {10^{ - 4}}\dfrac{{{K_{\rm{m}}}}}{{{\phi _{\rm{m}}}{C_{\rm{tm}}}\mu }}$

(9)

式中：

${m_{\rm{m}}}$—基岩拟压力函数，Pa/s；

${\eta}$—导压系数，m²/s。

式(8)、式(9)代入式(7)，得到基质渗流拟压力控制方程

$\dfrac{{{\partial ^2}{m_{\rm{m}}}(x,t)}}{{\partial {x^2}}} = \dfrac{{{\phi _{\rm{m}}}{C_{{\rm{tm}}}}\mu }}{{{K_{\rm{m}}}}}\dfrac{{\partial {m_{\rm{m}}}}}{{\partial t}} = \dfrac{1}{{{\eta _{\rm{m}}}}}\dfrac{{\partial {m_{\rm{m}}}}}{{\partial t}}$

(10)

定义无因次量

${m_{{\rm{mD}}}} = \dfrac{{{K_{\rm{m}}}{y_{\rm{e}}}}}{{{\rm{1422}}{q_{\rm{F}}}T}}\left[{{m_{\rm{m}}}({p_{\rm{i}}}) - {m_{\rm{m}}}(p)} \right]$

${t_{\rm{D}}} = 1.69 \times {10^{ - 7}}\dfrac{{{K_{\rm{f}}}}}{{{\phi _{\rm{f}}}\mu {C_{{\rm{tf}}}}{x_{\rm{F}}}^2}}t$

${x_{\rm{D}}} = \dfrac{x}{{{x_{\rm{F}}}}}$

${y_{\rm{D}}} = \dfrac{y}{{{x_{\rm{F}}}}}$

${\eta _{{\rm{mD}}}} = \dfrac{{{\eta _{\rm{m}}}}}{{{\eta _{\rm{f}}}}}$

${\eta _{\rm{f}}} = 1.82 \times {10^{ - 6}}\dfrac{{{K_{\rm{f}}}}}{{{\phi _{\rm{f}}}\mu {C_{{\rm{tf}}}}}}$

式中： $y_{\rm{e}}$—两条主裂缝间距的一半，m；

$q_{\rm{F}}$—每条主裂缝流入水平井的流量，m³/d；

$p_{\rm{i}}$—原始地层压力，MPa；

%$C_{{\rm{tf}}}$\dash

$x_{\rm{F}}$—主裂缝半长，m。

下标：

D—无因次，下同；

f—裂缝网络系统；

F—主裂缝系统。

则式(10)无因次化为

$\dfrac{{{\partial ^2}{m_{{\rm{mD}}}}({x_{\rm{D}}},{t_{\rm{D}}})}}{{\partial x_{\rm{D}}^2}} = \dfrac{1}{{{\eta _{{\rm{mD}}}}}}\dfrac{{\partial {m_{{\rm{mD}}}}}}{{\partial t{}_{\rm{D}}}}$

(11)

由拉普拉斯变换^[15-17]公式

$L[{m_{\rm{mD}}}({t_{\rm{D}}})] = {\overline m _{\rm{mD}}}({t_{\rm{D}}}) = \int_0^\infty {{{\rm{e}}^{ - s{t_{\rm{D}}}}}} {m_{\rm{mD}}}({t_{\rm{D}}}){\rm{d}}{t_{\rm{D}}}$

(12)

式中：

$\overline m _{\rm{D}}$—拉普拉斯空间下的无因次拟压力；

s—拉普拉斯变换算子，无因次。

在式(11)两边同乘以e$^{-st_{\rm{D}}}$并且相对于时间从0到无限大进行积分实现，即

$\int_0^\infty {{{\rm{e}}^{ - s{t_{\rm{D}}}}}} \dfrac{{{\partial ^2}{m_{{\rm{mD}}}}}}{{\partial x_{\rm{D}}^2}}{\rm{d}}{t_{\rm{D}}} = \dfrac{1}{{{\eta _{{\rm{mD}}}}}}\int_0^\infty {{{\rm{e}}^{ - s{t_{\rm{D}}}}}} \dfrac{{\partial {m_{{\rm{mD}}}}}}{{\partial {t_{\rm{D}}}}}{\rm{d}}{t_{\rm{D}}}$

(13)

由于${m_{\rm{mD}}}$是${x_{\rm{D}}}$和时间${t_{\rm{D}}}$的函数，对时间积分将自动地消除时间函数，只剩下${m_{\rm{mD}}}$对${x_{\rm{D}}}$的函数，所以对式(13)的左边简化为

$\int_0^\infty {{{\rm{e}}^{ - s{t_{\rm{D}}}}}} \dfrac{{{\partial ^2}{m_{{\rm{mD}}}}}}{{\partial x_{\rm{D}}^2}}{\rm{d}}{t_{\rm{D}}} = \dfrac{{{\partial ^2}\int_0^\infty {{{\rm{e}}^{ - s{t_{\rm{D}}}}}{m_{{\rm{mD}}}}{\rm{d}}{t_{\rm{D}}}} }}{{\partial x_{\rm{D}}^2}} = \dfrac{{{\partial ^2}{{\overline m }_{{\rm{mD}}}}({x_{\rm{D}}},s)}}{{\partial x_{\rm{D}}^2}}$

(14)

假设页岩气藏原始压力处处相等，那么${m_{{\rm{mD}}}}({t_{\rm{D}}} = 0)=0$， 式(14)中的积分变为

$\int_0^\infty {{{\rm{e}}^{ - s{t_{\rm{D}}}}}} \dfrac{{\partial {m_{{\rm{mD}}}}}}{{\partial {t_{\rm{D}}}}}{\rm{d}}{t_{\rm{D}}} = {{\rm{e}}^{ - s{t_{\rm{D}}}}}{m_{{\rm{mD}}}}({t_{\rm{D}}})\left| {_0^\infty } \right. + s\int_0^\infty {{{\rm{e}}^{ - s{t_{\rm{D}}}}}{m_{{\rm{mD}}}}({t_{\rm{D}}}){\rm{d}}} {t_{\rm{D}}}= \\ s\int_0^\infty {{{\rm{e}}^{ - s{t_{\rm{D}}}}}{m_{{\rm{mD}}}}({t_{\rm{D}}}){\rm{d}}} {t_{\rm{D}}} = s{\overline m _{{\rm{mD}}}}({x_{\rm{D}}},s)$

(15)

综上所述，基质渗流方程的拉普拉斯变换形式

$\dfrac{{{\partial ^2}{{\overline m }_{{\rm{mD}}}}({x_{\rm{D}}},s)}}{{\partial x_{\rm{D}}^2}} - \dfrac{s}{{{\eta _{{\rm{mD}}}}}}{\overline m _{{\rm{mD}}}}({x_{\rm{D}}},s) = 0$

(16)

则基岩线性渗流的数学模型为

$\dfrac{{{\partial ^2}{{\overline m }_{{\rm{mD}}}}({x_{\rm{D}}},s)}}{{\partial x_{\rm{D}}^2}} - \dfrac{s}{{{\eta _{{\rm{mD}}}}}}{\overline m _{{\rm{mD}}}}({x_{\rm{D}}},s) = 0$

(17)

又边界条件为

$\dfrac{{\partial {{\overline m }_{{\rm{mD}}}}}}{{\partial {x_{\rm{D}}}}}\left| {_{{x_{\rm{D}}} = {x_{{\rm{eD}}}}}} \right. = 0 %\sweqno$

${\overline m _{{\rm{mD}}}}\left| {_{{x_{\rm{D}}} = 1}} \right. = {\overline m _{{\rm{fD}}}}\left| {_{{x_{\rm{D}}} = 1}} \right. %\sweqno$

其中：

${x_{{\rm{eD}}}} = \dfrac{{{x_{\rm{e}}}}}{{{x_{\rm{F}}}}}$

得到渗流方程在拉普拉斯空间下的解

$\overline m _{{\rm{mD}}} = {\overline {m} _{{\rm{fD}}}} _{{x_{{\rm{D}} } = 1} } \dfrac{{\cosh [\sqrt {s/\eta _{{\rm{mD}}} } (x_{{\rm{eD}}} - x_{\rm{D}} )]}}{{\cosh [\sqrt {s/\eta _{{\rm{mD}}} } (x_{{\rm{eD}}} - 1)]}}$

(18)

式中：

$x_{\rm{e}}$—体积压裂区域半长，m。

2.2 裂缝网络线性渗流

裂缝网络中的气体主要沿y方向渗流

$ - \dfrac{\partial }{{\partial y}}\left( {\rho _{{\rm{gf}}} v_{\rm{f}} } \right) + 2q_{{\rm{mf}}} = \dfrac{\partial }{{\partial t}}(\phi _{\rm{f}} \rho _{{\rm{gf}}} )$

(19)

式中：

$v_{\rm{f}} = - \dfrac{{K_{\rm{f}} }}{\mu }\dfrac{{\partial p_{\rm{f}} }}{{\partial y}}$

$\rho _{{\rm{gf}}} = \dfrac{{p_{\rm{f}} M}}{{{\rm{R}}Tz}}$

$q_{{\rm{mf}}} = \rho _{{\rm{gm}}} q_{\rm{g}} = \dfrac{{p_{\rm{m}} }}{z}\dfrac{M}{{{\rm{R}}T}}\dfrac{{K_{\rm{m}} }}{\mu }\dfrac{{\partial p_{\rm{m}} }}{{\partial x}}\left| {_{x = x_{\rm{F}} } } \right.$

式中：

$q_{\rm{g}}$—基岩区流入裂缝网络区流量，m³/d。

代入式(19)，得到

$\dfrac{\partial }{{\partial y}}\left( {\dfrac{{K_{\rm{f}} }}{\mu }\dfrac{{p_{\rm{f}} }}{z}\dfrac{M}{{{\rm{R}}T}}\dfrac{{\partial p_{\rm{f}} }}{{\partial y}}} \right) + 2\dfrac{{p_{\rm{m}} M}}{{z{\rm{R}}T}}\dfrac{{K_{\rm{m}} }}{\mu }\dfrac{{\partial p_{\rm{m}} }}{{\partial x}}\left| {_{x = x_{\rm{F}} } } \right. =\phi _{\rm{f}} \dfrac{\partial }{{\partial t}}\left( {\dfrac{{p_{\rm{f}} M}}{{z{\rm{R}}T}}} \right)$

(20)

又

$\phi _{\rm{f}} \dfrac{\partial }{{\partial t}}\left( {\dfrac{{p_{\rm{f}} }}{z}} \right) = \phi _{\rm{f}} \dfrac{{p_{\rm{f}} }}{z}C_{{\rm{gf}}} \dfrac{{\partial p_{\rm{f}} }}{{\partial t}}$ $C_{{\rm{gf}}} = \dfrac{1}{{p_{\rm{f}} }} - \dfrac{1}{z}\dfrac{{\partial z}}{{\partial p_{\rm{f}} }}$

式中： $C_{{\rm{gf}}}$—裂缝网络区气体压综合系数，1/MPa。

定义拟压力函数

$m_{\rm{f}} = 2\int_0^{p_{_{\rm{f}} } } {\dfrac{p}{{\mu z}}} {\rm{d}}p$

(21)

式中：

$m_{{\rm{f}}}$—裂缝网络拟压力函数，1/s。

式(21)代入式(20)，得

$\dfrac{{\partial ^2 m_{\rm{f}} (y,t)}}{{\partial y^2 }} + 2\dfrac{{K_{\rm{m}} }}{{K_{\rm{f}} }}\dfrac{{\partial m_{\rm{m}} }}{{\partial x}}\left| {_{x = x_{\rm{F}} } } \right. = \dfrac{1}{{\eta _{\rm{f}} }}\dfrac{{\partial m_{\rm{f}} }}{{\partial t}}$

(22)

定义无因次参数

$m_{{\rm{fD}}} = \dfrac{{K_{\rm{f}} x_{\rm{F}} }}{{{\rm{1422}}q_{\rm{F}} T}}\left[{m_{\rm{f}} (p_{\rm{i}} ) - m_{\rm{f}} (p)} \right] %\sweqno$

$\eta _{{\rm{fD}}} = \dfrac{{\eta _{\rm{f}} }}{{\eta _{\rm{f}} }} = 1 %\sweqno$

则式(22)变为

$\dfrac{{\partial ^2 m_{{\rm{fD}}} (y_{\rm{D}} ,t_{\rm{D}} )}}{{\partial y_{\rm{D}}^2 }} + 2\dfrac{{K_{\rm{m}} }}{{K_{\rm{f}} }}\dfrac{{\partial m_{{\rm{mD}}} }}{{\partial x_{\rm{D}} }}\left| {_{x_{\rm{D}} = 1} } \right. = \dfrac{{\partial m_{{\rm{fD}}} }}{{\partial t_{\rm{D}} }}$

(23)

同基质渗流一样，进行拉普拉斯变换，得到拉普拉斯空间的公式为

$\dfrac{{\partial^2 \overline m _{y{\rm{D}}} (y_{\rm{D}} ,s)}}{{\partial y_{\rm{D}}^2 }}\!+\! 2\dfrac{{K_{\rm{m}} }}{{K_{\rm{f}} }}\dfrac{{\partial \overline m _{{\rm{mD}}} }}{{\partial x_{\rm{D}} }}|_{x_{\rm{D}} = 1}\!-\!s\overline m _{{\rm{fD}}} (y_{\rm{D}} ,s)\!=\!0$

(24)

令

$y_{{\rm{eD}}} = \dfrac{{y_{\rm{e}} }}{{x_{\rm{F}} }}$

$C_{{\rm{fD}}} = \dfrac{{K_{\rm{f}} x_{\rm{F}} }}{{K_{\rm{m}} y_{\rm{e}} }}$

则式(24)表示为

$\dfrac{{\partial ^2\overline m_{y{\rm{D}}}(y_{\rm{D}},s)}}{{\partial y_{\rm{D}}^2}}+ \dfrac{2}{{y_{{\rm{eD}}}C_{{\rm{fD}}}}}\dfrac{{\partial \overline m_{{\rm{mD}}}}}{{\partial x_{\rm{D}}}}|_{x_{\rm{D}}=1}-s\overline m_{{\rm{fD}}}(y_{\rm{D}},s)=0$

(25)

由基岩渗流解析解公式(18)得到

$\dfrac{{\partial \overline m _{{\rm{mD}}} }}{{\partial x_{\rm{D}} }}\left| {_{x_{\rm{D}} = 1} = - \beta _0 \overline m _{{\rm{fD}}} \left| {_{x_{\rm{D}} = 1} } \right.} \right.$

(26)

式中：

${{\beta }_{\text{o}}}=\sqrt{s\text{ }/{{\eta }_{\text{mD}}}}\tanh \left[ \sqrt{\frac{s}{{{\eta }_{\text{mD}}}}}({{x}_{\text{eD}}}-1) \right]$

式(25)变为

$\dfrac{{\partial ^2 \overline m _{{\rm{fD}}} }}{{\partial y_{\rm{D}} ^2 }} - \alpha _{\rm{o}} \overline m _{{\rm{fD}}} = 0$

(27)

式中：

$\alpha _{\rm{o}} = \dfrac{{2\beta _{\rm{o}} }}{{C_{{\rm{fD}}} y_{{\rm{eD}}} }} + s$

边界条件为

$\left\{ \begin{array}{l} \dfrac{{\partial \overline m _{{\rm{fD}}} }}{{\partial y_{\rm{D}} }}\left| {_{y_{\rm{D}} = y_{{\rm{eD}}} } } \right. = 0 \\ \overline m _{{\rm{fD}}} \left| {_{y_{\rm{D}} = w_{\rm{D}} /2} = } \right.\overline m _{{\rm{FD}}} \left| {_{y_{\rm{D}} = w_{\rm{D}} /2} } \right. \\ \end{array} \right.$

(28)

式中：

$w_{\rm{D}}$—主裂缝宽度的无因次值，定义为

$w_{\rm{D}} = \dfrac{w}{{x_{\rm{F}} }}$

(29)

由此得到裂缝网络渗流方程的解为

$\overline m _{{\rm{fD}}} = (\overline m _{{\rm{FD}}} \left| {_{y_{\rm{D}} = w_{\rm{D}} /2} } \right.)\dfrac{{\cosh [\sqrt {\alpha _{\rm{o}} } (y_{{\rm{eD}}} - y_{\rm{D}} )]}}{{\cosh [\sqrt {\alpha _{\rm{o}} } (y_{{\rm{eD}}} - w_{\rm{D}} /2)]}}$

(30)

2.3 主裂缝线性渗流

同理应用拉普拉斯空间变换可得到主裂缝渗流方程

$\dfrac{{\partial ^2 \overline m _{{\rm{FD}}} }}{{\partial x_{\rm{D}}^2 }} + \dfrac{4}{{C_{{\rm{FD}}} }}\dfrac{{\partial \overline m _{{\rm{FD}}} }}{{\partial y_{\rm{D}} }}\left| {_{y_{\rm{D}} = w_{\rm{D}} /2} } \right. - \dfrac{s}{{\eta _{{\rm{FD}}} }}\overline m _{{\rm{FD}}} = 0$

(31)

式中：

$\eta _{{\rm{FD}}} = \dfrac{{\eta _{\rm{F}} }}{{\eta _{\rm{f}} }}$

$\eta _{\rm{F}} = 1.82 \times 10^{ - 6} \dfrac{{K_{\rm{F}} }}{{\phi _{\rm{F}} \mu C_{{\rm{tF}}} }}$

$C_{{\rm{FD}}} = \dfrac{{K_{\rm{F}} w}}{{K_{\rm{f}} x_{\rm{F}} }}$

$m_{{\rm{FD}}} = \dfrac{{K_{\rm{F}} x_{\rm{F}} }}{{{\rm{1422}}q_{\rm{F}} T}}\left[{m_{\rm{F}} (p_{\rm{i}} ) - m_{\rm{F}} (p)} \right]$

由裂缝网络渗流方程解析解公式(30)，得到

$\dfrac{{\partial \overline m _{{\rm{fD}}} }}{{\partial y_{\rm{D}} }}\left| {_{y_{\rm{D}} = w_{\rm{D}} /2} } \right. = - \beta _{\rm{F}} \overline m _{{\rm{FD}}} \left| {_{y_{\rm{D}} = w_{\rm{D}} /2} } \right.$

(32)

式中：

$\beta _{\rm{F}} = \sqrt {\alpha _{\rm{o}} } \tanh [\sqrt {\alpha _{\rm{o}} } (y_{{\rm{eD}}} - w_{\rm{D}} /2)]$

所以式(31)变化为

$\dfrac{{\partial ^2 \overline m _{{\rm{FD}}} }}{{\partial x_{\rm{D}}^2 }} - \alpha _{\rm{F}} \overline m _{{\rm{FD}}} = 0$

(33)

式中：

$\alpha _{\rm{F}} = \dfrac{{4\beta _{\rm{F}} }}{{C_{{\rm{FD}}} }} + \dfrac{s}{{\eta _{{\rm{FD}}} }}$

边界条件

$\dfrac{{\partial \overline m _{{\rm{FD}}} }}{{\partial x_{\rm{D}} }}\left| {_{x_{\rm{D}} = 1} } \right. = 0$ $\dfrac{{\partial \overline m _{{\rm{FD}}} }}{{\partial x_{\rm{D}} }}\left| {_{x_{\rm{D}} = 0} } \right. = - \dfrac{\pi }{{C_{{\rm{FD}}} s}}$

由此，得到主裂缝渗流方程的解

$\overline m _{{\rm{FD}}} = \dfrac{\pi }{{C_{{\rm{FD}}} s^{1/4} \sqrt {\alpha _{\rm{F}} } }}\dfrac{{\cosh [\sqrt {\alpha _{\rm{F}} } (1 - x_{\rm{D}} )]}}{{\sinh \left( {\sqrt {\alpha _{\rm{F}} } } \right)}}$

(34)

$x_{\rm{D}}$=0时，得到井底压力公式

$\overline m _{{\rm{wD}}} = \dfrac{\pi }{{C_{{\rm{FD}}} s^{1/4} \sqrt {\alpha _{\rm{F}} } }}\dfrac{{\cosh (\sqrt {\alpha _{\rm{F}} } )}}{{\sinh (\sqrt {\alpha _{\rm{F}} } )}} = \\ \dfrac{\pi }{{C_{{\rm{FD}}} s^{1/4} \sqrt {\alpha _{\rm{F}} } \tanh (\sqrt {\alpha _{\rm{F}} } )}} $

(35)

3 渗流模型的求解

对于双曲正切函数有

$\tanh (x) = x - \dfrac{{x^3 }}{3} + \dfrac{{2x^5 }}{{15}} + \ldots \ldots$

(36)

约去无穷小量，$\tanh (x) \approx x$ 从而$\tanh (\sqrt {\alpha _{\rm{F}} } ) \approx \sqrt {\alpha _{\rm{F}} }$，则

$\overline m _{{\rm{wD}}} \approx \dfrac{\pi }{{C_{{\rm{FD}}} s^{1/4} \alpha _{\rm{F}} }}$

(37)

将$\alpha _{\rm{F}}$的表达式代入式(37)，整理化简得到

$\overline m _{{\rm{wD}}} = \dfrac{{\pi \eta _{{\rm{FD}}} }}{{s^{1/4} (4\beta _{\rm{F}} \eta _{{\rm{FD}}} + sC_{{\rm{FD}}} )}}$

(38)

再代入对应参数并化简，得到

$\overline m _{{\rm{wD}}} = \dfrac{{\pi y_{{\rm{eD}}} \eta _{\rm{m}} }}{{4C_{{\rm{FD}}} x_{\rm{e}} K_{\rm{m}}^2 (y_{{\rm{eD}}} - w_{\rm{D}} /2)[2\eta _{\rm{f}} (x_{{\rm{eD}}} - 1) + C_{{\rm{RD}}} y_{{\rm{eD}}} \eta _{\rm{m}} + C_{{\rm{FD}}} ^2 y_{{\rm{eD}}} \eta _{\rm{m}}]s^{1.25} }}$

(39)

其中

$x_{{\rm{eD}}} = \dfrac{{x_{\rm{e}} }}{{x_{\rm{F}} }}$

由Laplace反变换：$\dfrac{{\Gamma (m + 1)}}{{s^{m + 1} }} = t^m (m > - 1)$，得到

$m_{{\rm{wD}}} = \dfrac{{\pi y_{{\rm{eD}}} \eta _{\rm{m}} t_{\rm{D}}^{1/4} }}{{4C_{{\rm{FD}}} x_{\rm{e}} K_{\rm{m}}^2 (y_{{\rm{eD}}} - w_{\rm{D}} /2)[2\eta _{\mathop{\rm f}\nolimits} (x_{{\rm{eD}}} - 1) + C_{{\rm{RD}}} y_{{\rm{eD}}} \eta _{\rm{m}} + C_{{\rm{FD}}} ^2 y_{{\rm{eD}}} \eta _{\rm{m}}]\Gamma (1.25)}}$

(40)

由于$\Gamma (1.25) \approx 0.908$，对式(40)进行单位换算，并化简

$m_{{\rm{wD}}} = \dfrac{{2.04 \times 10^{ - 10} \pi y_{{\rm{eD}}} \eta _{\rm{m}} t_{\rm{D}}^{1/4} }}{{C_{{\rm{FD}}} x_{\rm{e}} K_{\rm{m}}^2 (y_{{\rm{eD}}} - w_{\rm{D}} /2)[2\eta _{\rm{f}} (x_{{\rm{eD}}} - 1) + C_{{\rm{RD}}} y_{{\rm{eD}}} \eta _{\rm{m}} + C_{{\rm{FD}}} ^2 y_{{\rm{eD}}} \eta _{\rm{m}}]}}$

(41)

又$m_{{\rm{wD}}} = \dfrac{{K_{\rm{f}} h}}{{1422q_{\rm{F}} T}}[m(p_{\rm{i}} ) - m(p_{{\rm{wf}}} )]$ 代入拟压力函数定义式，得到

$q_{\rm{F}}\!=\!\dfrac{{2K_{\rm{f}} h(p_{\rm{i}}^2\!-\!p_{{\rm{wf}}}^2)}}{{1422\mu zm_{{\rm{wD}}} T}} \!= \\ \!\dfrac{{6.89 \times 10^6 (p_{\rm{i}}^2\!-\!p_{{\rm{wf}}}^2 )K_{\rm{f}} hC_{{\rm{FD}}} x_{\rm{e}} K_{\rm{m}}^2 (y_{{\rm{eD}}}\!-\!w_{\rm{D}} /2)[2\eta _{\rm{f}} (x_{{\rm{eD}}}\!-\! 1)\!+\!C_{{\rm{fD}}} y_{{\rm{eD}}} \eta_{\rm{m}} + C_{{\rm{FD}}} ^2 y_{{\rm{eD}}} \eta _{\rm{m}}]}}{{\pi y_{{\rm{eD}}} \eta_{\rm{m}} T\mu z}}t_{\rm{D}}^{ - 1/4} $

(42)

则水平井总的产量为

$q = q_{\rm{H}} n_{\rm{F}}$

也可得到井底压力与生产时间的关系式

$p_{{\rm{wf}}} = \left\{ {p_{\rm{i}}^2 - \dfrac{{6.89 \times 10^{ - 8} \pi y_{{\rm{eD}}} \eta _m^{} T\mu zq_{\rm{g}} t_{\rm{D}}^{1/4} }}{{C_{{\rm{FD}}} x_{\rm{e}} K_{\rm{m}}^2 K_{\rm{f}} h(y_{{\rm{eD}}} - w_{\rm{D}} /2)[\eta _{\rm{f}} (x_{{\rm{eD}}} - 1) + C_{{\rm{fD}}} y_{{\rm{eD}}} \eta _{\rm{m}} + C_{{\rm{FD}}} ^2 y_{{\rm{eD}}} \eta _{\rm{m}}]}}} \right\}^{1/2}$

(43)

4 模型的应用

某页岩气储层基本参数如表 1所示。由所建立三线性渗流模型进行产能分析，得到以下3个方面的影响。

4.1 解吸/吸附机理对页岩气藏生产的影响

使用表 1中页岩气藏相关参数对页岩气储层进行模拟，井底压力为5 MPa时分析解吸吸附机理对页岩气藏开采的影响。

Langmuir体积常数是指单位体积的页岩气储层吸收气体的最大体积，Langmuir压力常数是指吸附气体体积达最大吸附量的一半时对应的压力，它们表征了页岩气藏的吸附能力。由图 2可知：气井开始生产短时间内产量会迅速降低，之后下降速率变慢；改变Langmuir体积可以明显地改变气井的稳产年限，例如$V_{\rm{L}}$=3 m³/t时，气井以10 000 m³/d可生产5 a，当$V_{\rm{L}}$=5 m³/t时，气井以10 000 m³/d可以稳产14 a。

图 2、图 3表明：Langmuir压力的变化也会影响页岩气井的稳产时间，但效果相对不明显；页岩气藏压裂水平井初期产量在9 000$\sim$20 000 m³/d，但产量递减较快，10 000 m³/d稳产最大达到5 a，之后以3 000$\sim$6 000 m³/d稳定生产可延续几十年，体现了页岩气藏低产量、生产时间长的特点。

4.2 人工裂缝半长对气井生产的影响

利用上述数学模型及表 1中的相关数据模拟了不同人工裂缝半长条件下气井产量的变化(如图 4)，由图可以看出：当裂缝半长为100 m时，气井初期产量在11 000 m³/d；当裂缝半长增加到400 m左右时，气井初期产量可达到20 000 m³/d，同时，裂缝半长较大者后期产量也明显较高，因此，理论上认为水平井压裂时，裂缝体积缝网越大越有利于页岩气的解吸和渗流。

4.3 人工裂缝渗透率对气井生产的影响

人工裂缝渗透率增大，减小了气体的流动阻力，因此产量随之增加，如图 5所示。当人工裂缝渗透率为2 500 mD时，13 a后气井产量降到10 000 m³/d以下；当渗透率为1 000 mD，气井产量还达不到10 000 m³/d，这说明人工裂缝渗透率增加，稳产时间随之增大。然而，压裂高导流能力的裂缝成本相对较高，从油气田开采经济最优化的原则来说，并不是裂缝渗透率越高越好。

5 结 论

(1) 页岩气压裂水平井渗流很复杂，多尺度的渗流机理计算产能势必更复杂，本文在新设计的页岩气藏渗流模型基础上建立的数学模型及其解析解提供了相对简单的计算方法，并对某页岩气藏经过计算分析，得到的结果符合生产规律，因此，该方法可应用于页岩气藏压裂水平井产能的计算分析。

(2) 解吸/吸附现象的存在是页岩气藏与常规气藏最重要的区别，本文在建立的模型中加入了该元素，通过计算结果形象地描述了解吸吸附对页岩气藏开采的影响： Langmuir压力常数增加，页岩气井日产量增加，增幅减小；Langmuir体积常数越大， 表示页岩气的吸附量越大，页岩气越易解吸，产量越大。

(3) 计算结果表明人工裂缝渗透率及裂缝半长的增加都有助于提高页岩气井的日产量，但是满足这两个条件时付出的代价和成本也会很高，从油气藏经济最大化原则方面讲，一味增加人工裂缝渗透率和半长是不科学的，因此在开采页岩气藏时，要寻求最优的人工裂缝渗透率和半长等气藏参数。