经过多段压裂改造后的致密气藏水平井的产能计算方法

摘要

本发明公开了一种经过多段压裂改造后的致密气藏水平井的产能计算方法，它包括以下步骤：建立物理模型；得到n条裂缝中的任一条引起的无量纲压力PD与该条裂缝无量纲裂缝流量qD(α)之间的关系式；根据气体在裂缝中的流动关系以及裂缝与地层之间的边界耦合关系，得到n条裂缝中的任一条在与井筒接触处产生的压力扰动Pwfn与标况下该条裂缝流量的Qscn的关系式；根据气体在井筒中的流动关系以及裂缝与井筒之间的边界耦合关系，得到n条裂缝中的任一条在与井筒接触处产生的压力扰动Pwfi、与该条裂缝相邻的裂缝在与井筒接触处产生的压力扰动Pwfi-1和这两条裂缝在井筒内之间的流量Qsci的关系式；利用数值迭代法得到该水平井的产能。

说明书

技术领域

本发明涉及了一种产能计算方法，尤其涉及了一种经过多段压裂改造后的致密气藏水平 井的产能计算方法。本发明涉及了一种产能计算系统，尤其涉及了一种经过多段压裂改造后 的致密气藏水平井的产能计算系统。

背景技术

致密气藏（Tight Gas）是指渗透率小于0.1毫达西（mD）的砂岩地层天然气藏。致密 气藏作为一种重要的天然气资源，已经逐渐成为天然气产量的主要增长点。现有技术中，往 往采用水平井对致密气藏进行开采。水平井是指井斜角达到或接近90°且井身沿着水平方向 钻进一定长度的井。水平井经过多段压裂改造后，形成多条形态不同的横向裂缝（以下简称 “裂缝”）。裂缝沿垂直于井筒的方向起裂。裂缝大大增加了气井与地层的接触面积，同时改 善了井底周围储层的渗流条件，增加了油层泄油面积。致密气藏的气体从地层孔隙中流入裂 缝，在裂缝中流向井筒，进而沿井筒流向井口。

为了对经过多段压裂改造后的致密气藏水平井进行产能评估，本领域技术人员进行了不 断深入的研究。在裂缝有限导流研究方面，裂缝模型从单裂缝发展到多裂缝再到体积压裂， 求解方法从解析法发展到半解析解再到数值法，模拟精度不断提高。压裂水平井产能评价方 面，范子菲利用水平井筒的流动特征建立了储层与井筒的流动耦合模型，李晓平应用体积平 衡原理修正了耦合模型。但是他们所提出的水平井模型是以裸眼、割缝衬管或割缝筛管的方 式完井。相较于其他完井方式，分段压裂后致密气藏水平井的裂缝与井筒相比流动面积要大 得多，所以需要充分考虑人工裂缝流动对气井产能的影响。

发明内容

本发明的发明目的在于，提供了一种经过多段压裂改造后的致密气藏水平井的产能计算 方法，该方法能够提高经过多段压裂改造后的致密气藏水平井产能预测的准确性。

本发明公开了一种经过多段压裂改造后的致密气藏水平井的产能计算方法，其特征在 于：它包括以下步骤

步骤A），建立物理模型，所述物理模型具有以下定义:A1)地层均质等厚，地层的俯视面 为矩形，该矩形具有四条闭合且等压的边界，所述矩形的宽为x_e，该值通过对地层进行试井 解释得到，所述矩形的长为y_e，该值通过对地层进行试井解释得到；A2)具有n条裂缝，所 有n条裂缝完全贯穿地层，其中n=1,2,3……，第1条裂缝位于该水平井井筒的根端，第i 条裂缝逐渐向该水平井井筒的趾端排布,其中i=1,2,3……,n；

所述物理模型定义如下无量纲量：

$P_D=\frac{78.55\mathrm{kh}(P_i^2-P^2)}{\mu {\mathrm{ZTQ}}_{\mathrm{sc}}},q_D=\frac{{2x}_f{q}_{\mathrm{sc}}}{Q_{\mathrm{sc}}},j_D=\frac{j}{x_f}(j=x,y),C_{\mathrm{fD}}=\frac{k_f{w}_f}{{\mathrm{kx}}_f}$

式中：

P_D代表无量纲压力；P_i代表原始地层压力；P代表地层压力；T代表地层温度；k代表地 层渗透率；h代表地层厚度；μ代表气体粘度；Z代表气体偏差系数，其通过室内实验得到； Q_sc代表标况下裂缝流量；q_D代表无量纲裂缝流量；q_sc为标况下单位裂缝长度流量；j_D为无量 纲长度；x_f为裂缝半长；C_fD为无量纲裂缝导流能力；k_f为裂缝渗透率；w_f为裂缝宽度；

步骤B），基于物理模型，根据气体在地层中的渗流规律，得到n条裂缝中的任一条引起 的无量纲压力P_D与该条裂缝无量纲裂缝流量q_D(α)之间的关系式，

$\left(\begin{array}{c}{P}_D(x_D,y_D;x_{\mathrm{wD}},y_{\mathrm{wD}})\\ =2\underset{x_{\mathrm{wD}}-1}{\overset{x_{\mathrm{wD}}+1}{\int }}\left\{\underset{m=1}{\overset{\infty }{\Sigma }}\frac{q_D\left(\alpha \right)}{\mathrm{m\pi }}\sin \frac{{\mathrm{m\pi x}}_D}{x_{\mathrm{eD}}}\sin \frac{\mathrm{m\pi \alpha }}{x_{\mathrm{eD}}}\frac{\cosh \frac{\mathrm{m\pi }(y_{\mathrm{eD}}-|y_D-y_{\mathrm{wD}}\left|\right)}{x_{\mathrm{eD}}}-\cosh \frac{\mathrm{m\pi }(y_{\mathrm{eD}}-|y_D+y_{\mathrm{wD}}\left|\right)}{x_{\mathrm{eD}}}}{\sinh \frac{{\mathrm{m\pi y}}_{\mathrm{eD}}}{x_{\mathrm{eD}}}}\right\}\mathrm{d\alpha }\end{array}\right)$

式中：x_w，y_w为裂缝中心坐标，其通过物理模型的定义得到；

步骤C），基于物理模型，根据气体在裂缝中的流动关系以及裂缝与地层之间的边界耦合 关系，得到n条裂缝中的任一条在与井筒接触处产生的压力扰动P_wfn与标况下该条裂缝流量 的Q_scn的关系式，

$\left(\begin{array}{c}\frac{78.55\mathrm{kh}(P_i^2-P_{\mathrm{wfn}}^2)}{{\mathrm{\mu ZTQ}}_{\mathrm{scn}}}=\frac{1}{C_{\mathrm{fD}}}\frac{h}{x_f}[\ln \frac{h}{{2r}_w}-\frac{\pi }{2}]+f\left(C_{\mathrm{fD}}\right)\\ +2\left\{\underset{m=1}{\overset{\infty }{\Sigma }}\frac{x_e^2}{m^3{\pi }^2{x}_f^2}{\sin }^2\frac{{\mathrm{m\pi x}}_f}{x_e}{\sin }^2\frac{{\mathrm{m\pi x}}_w}{x_e}\frac{\cosh \frac{{\mathrm{m\pi y}}_e}{x_e}-\cosh \frac{\mathrm{m\pi }(y_e-{2y}_w)}{x_e}}{\sinh \frac{{\mathrm{m\pi y}}_e}{x_e}}\right\}\end{array}\right),$

式中：P_wfn为n条裂缝中的任一条在与井筒接触处产生的压力扰动；r_w为水平井筒半径；

步骤D），基于物理模型，根据气体在井筒中的流动关系以及裂缝与井筒之间的边界耦合 关系，得到n条裂缝中的任一条在与井筒接触处产生的压力扰动P_wfi、与该条裂缝相邻的裂 缝在与井筒接触处产生的压力扰动P_wfi-1和这两条裂缝在井筒内之间的流量Q_sci的关系式，

$p_{\mathrm{wfi}}^2-P_{\mathrm{wfi}-1}^2=9\times {10}^{-12}\frac{{\mathrm{ZT\gamma }}_g}{r_w^5}{f}_i{d}_i{\left(\underset{j=i}{\overset{n}{\Sigma }}{Q}_{\mathrm{scj}}\right)}^2$

其中，P_wf0=P_wf；

式中：Z为气体偏差系数；γ_g为气体相对密度；f为摩阻系数；d为裂缝间距；

步骤E），估算第n条裂缝的流量Q_scn，测量得到第1条裂缝在与井筒接触处产生的压力 扰动P_wf，利用数值迭代法得到该水平井的产能。

优选地，在步骤C)中，气体在裂缝中的流动关系要考虑到裂缝的渗流作用、无限导流裂 缝与有限导流裂缝之间的差异以及井筒的径向聚流效应。

优选地，在步骤E)中，还包括以下步骤：

E1）根据实际生产情况估算第n条裂缝流量的最大值Q_scn(max)和最小值Q_scn(min)，取最大值 与最小值的算术平均Q_scn(mid)=0.5×[Q_scn(max)+Q_scn(min)]；

E2）根据步骤C）中的公式计算得到对应的P_wfn(max)、P_wfn(min)和P_wfn(mid)，根据步骤D）中公 式得到对应的P_wfn-1(max)、P_wfn-1(min)和P_wfn-1(mid)，根据步骤C）中的公式计算得到对应的Q_scn-1(max)、 Q_scn-1(min)、Q_scn-1(mid)；

E3）重复步骤E2）得到P_wf0(max)、P_wf0(min)、P_wf0(mid)；

E4）将P_wf0(mid)与P_wf进行差值比较，如果P_wf0(mid)与P_wf的差值符合精度要求则判断计算值 P_wf0(mid)正确，如果P_wf0(mid)与P_wf的差值不符合精度要求，则如果(P_wf0(max)-P_wf)×(P_wf0(min)-P_wf)<0， 那么Q_scn(min)=Q_scn(mid)，否则Q_scn(max)=Q_scn(mid)；

E5）重复步骤E1）至E4），直至P_wf0(mid)与P_wf的差值符合精度要求。

优选地，在步骤D）中，所述摩擦系数f由以下公式计算得到：

$f_i=[1.14-21g(\frac{e}{1000D}+\frac{21.25}{R_{\mathrm{ei}}^{0.9}}){]}^{-2},$

其中，e为井筒内壁粗糙度。

优选地，该水平井的产能为其中，j=1,2,3……，n。

优选地，所述P_i通过对未开采的地层测量得到，所述P通过对开采后的地层进行测量得 到，所述T通过对开采后的地层的温度进行测量得到，所述k通过室内实验或试井解释得到， 所述h通过测井解释得到，所述μ通过室内实验得到，所述Z通过室内实验得到，所述x_f通过试井解释得到，所述k_f通过试井解释得到，所述w_f通过试井解释得到。

优选地，所述Z通过室内实验得到；所述γ_g通过室内实验得到，所述d通过压裂设计 数据得到。

本发明还公开了一种采用上述计算方法的产能计算系统，它包括

建模单元，其用于建立物理模型；

第一计算单元，其用于根据气体在地层中的渗流规律得到n条裂缝中的任一条引起的无 量纲压力P_D与该条裂缝无量纲裂缝流量q_D(α)之间的关系式；

第二计算单元，其用于根据气体在裂缝中的流动关系以及裂缝与地层之间的边界耦合关 系得到n条裂缝中的任一条在与井筒接触处产生的压力扰动P_wfn与标况下该条裂缝流量的Q_scn的关系式；

第三计算单元，其用于根据气体在裂缝中的流动关系以及裂缝与地层之间的边界耦合关 系得到n条裂缝中的任一条在与井筒接触处产生的压力扰动P_wfi、与该条裂缝相邻的裂缝在 与井筒接触处产生的压力扰动P_wfi-1和这两条裂缝在井筒内之间的流量Q_sci的关系式；

第四计算单元，其用于接收估算第n条裂缝的流量Q_scn，其用于接收由测量得到第1条 裂缝在与井筒接触处产生的压力扰动P_wf，其用于采用数值迭代法得到该水平井的产能。

本发明人引入致密气藏压裂水平井产能单压裂段评价思路，以单压裂段为单元，应用质 量守恒原理将储层内渗流、裂缝内变质量渗流和井筒变流量管流耦合，研究水平井不同压裂 段的流量变化，推导建立适合于致密气藏水平井压裂后产能评价的理论公式，给出相应的迭 代算法并结合实例进行分析解释，从而形成一套致密气藏分段压裂水平井产能评价的新方 法。本发明基于渗流力学理论，研究气体在地层、裂缝、井筒三个独立系统内的流动规律， 借助质量守恒原理将储层内渗流、裂缝内变质量渗流和井筒变流量管流耦合，并通过迭代算 法确定多段压裂水平井产能。

附图说明

图1A显示了水平井的物理模型。

图1B显示了图1A中单段裂缝中的气体流动示意图。

图2显示了气体在裂缝中渗流的示意图。

图3显示了不同影响参数下导流能力C_fD与井底压力P_wD的变化关系。

图4显示了有限导流裂缝井底压力与导流能力变化关系。

图5显示了裂缝中的井筒聚流作用。

图6显示了水平井筒中水平管流剖面图。

图7显示了水平井的产量与井底压力变化关系图。

图8显示了不同井底流压下各裂缝沿水平井筒的压力分布。

图9显示了该井的井眼轨迹。

具体实施方式

下面结合附图对本发明的较佳实施例进行详细阐述，以使本发明的优点和特征能更易于 被本领域的技术人员理解，从而对本发明的保护范围作出更为清楚明确的界定。

本发明中产能计算方法的第一实施例，经过多段压裂改造后的致密气藏水平井的产能计 算方法，它包括以下步骤：建立物理模型。基于物理模型，根据气体在地层中的渗流规律， 得到n条裂缝中的任一条引起的无量纲压力P_D与该条裂缝无量纲裂缝流量q_D(α)之间的关系 式。基于物理模型，根据气体在裂缝中的流动关系以及裂缝与地层之间的边界耦合关系，得 到n条裂缝中的任一条在与井筒接触处产生的压力扰动P_wfn与标况下该条裂缝流量的Q_scn的 关系式。基于物理模型，根据气体在井筒中的流动关系以及裂缝与井筒之间的边界耦合关系， 得到n条裂缝中的任一条在与井筒接触处产生的压力扰动P_wfi、与该条裂缝相邻的裂缝在与 井筒接触处产生的压力扰动P_wfi-1和这两条裂缝在井筒内之间的流量Q_sci的关系式。估算第n 条裂缝的流量Q_scn，测量得到第1条裂缝在与井筒接触处产生的压力扰动P_wf，利用数值迭代 法得到该水平井的产能。

具体的，在步骤A）中，建立物理模型，所述物理模型具有以下定义:A1)地层均质等厚。 地层的俯视面为矩形。该矩形具有四条闭合且等压的边界。所述矩形的宽为x_e，该值通过对 地层进行试井解释得到，所述矩形的长为y_e，该值通过对地层进行试井解释得到。A2)具有 n条裂缝，所有n条裂缝完全贯穿地层，其中n=1,2,3……，第1条裂缝位于该水平井的根 端，第i条裂缝逐渐向该水平井的趾端排布,其中i=1,2,3……,n。

所述物理模型定义如下无量纲量：

$P_D=\frac{78.55\mathrm{kh}(P_i^2-P^2)}{\mu {\mathrm{ZTQ}}_{\mathrm{sc}}},q_D=\frac{{2x}_f{q}_{\mathrm{sc}}}{Q_{\mathrm{sc}}},j_D=\frac{j}{x_f}(j=x,y),C_{\mathrm{fD}}=\frac{k_f{w}_f}{{\mathrm{kx}}_f}$

在上述式中：P_D为无量纲压力。P_i为原始地层压力，其通过对未开采的地层测量得到， 单位为兆帕（MPa）。P为地层压力，其通过对开采后的地层进行测量得到，单位为兆帕（MPa）。 T为地层温度，其通过对开采后的地层的温度进行测量得到，单位为开尔文（K）；k为地层 渗透率，其通过室内实验或试井解释得到，其单位为达西（D）。h为地层厚度，其通过测井 解释得到，其单位为米(m).μ为气体粘度，其通过室内实验得到，其单位为毫帕·秒 （mPa·s）。Z为气体偏差系数，其通过室内实验得到。Q_sc为标况下裂缝流量，其单位为10⁴* 立方米/日(10⁴m³/d).q_D为无量纲裂缝流量。q_sc为标况下单位裂缝长度流量，其单位为10⁴*立 方米/日/米(10⁴m³/d/m)。j_D为无量纲长度。x_f为裂缝半长，其通过试井解释得到，其单位为 米（m）。C_fD为无量纲裂缝导流能力。k_f为裂缝渗透率，其通过试井解释得到，其单位为达西 （D）。w_f为裂缝宽度，其通过试井解释得到，其单位为米（m）。

试井解释就是以渗流力学理论为基础，通过对油井测试信息的研究，确定反映测试井和 地层特性的各种物理参数、生产能力等的方法。

测井解释确定测井信息与地质信息之间应用的关系，采用正确的方法把测井信息加工成 地质信息。

图1A显示了水平井的物理模型。如图1A所示，水平井的井筒沿着水平方向延伸。裂缝 沿纵向垂直于水平井井筒。图1B显示了图1A中单条裂缝中的气体流动示意图。在图1B中， 气体从地层中沿双曲线形态流线流向裂缝。气体从地层流入到裂缝后，通过渗流作用运动至 水平井筒。由于裂缝完全贯穿地层，可以认为裂缝内的流动为一维线性流动，而且这种流动 是一种沿裂缝长度变化的变质量流动。由于井筒内径远远大于地层和裂缝中的流动通道尺 寸，因此井筒内的气体流动按单相管流进行计算。

在该物理模型中，地层内的气体遵循达西定律（Darcy’s Law），裂缝内为变质量的Darcy 流动，水平井筒内为变质量的管流。三种流动方式相互干扰，通过相交的边界耦合。

在步骤B）中，基于物理模型，根据气体在地层中的渗流规律，得到n条裂缝中的任一 条引起的无量纲压力P_D与该条裂缝无量纲裂缝流量q_D(α)之间的关系式。

如图1B所示，在裂缝周边的地层中气体的流线形态近似于双曲线。在地层呈双曲线流 线以外的区域，气体的流线形态呈辐射状，表现为拟径向流特征。

根据拉普拉斯变换公式得到式（1）：

$\frac{{\partial }^2{P}_D}{{\partial x}_D^2}+\frac{{\partial }^2{P}_D}{{\partial y}_D^2}+q_D\left(x_D\right)\delta (y_D-y_{\mathrm{wD}})=0---\left(1\right)$

基于物理模型中定义的地层的四条边界等压得到式（2）和（3）：

P_D(x_D,0)=P_D(x_D,x_eD)                       （2）

P_D(0,y_D)=P_D(y_eD,y_D)                      （3）

式中：x_w，y_w为裂缝中心坐标，其通过物理模型的定义得到，其单位为米（m）。δ为狄 里克函数。

对式（1）中的各个变量沿x_D和y_D方向做傅里叶（Fourier）有限正弦积分，分别记为：

${\hat{P}}_D={\int }_0^{x_{\mathrm{eD}}}{P}_D\sin \left({\beta }_m{x}_D\right){\mathrm{dx}}_D---\left(4\right)$

${\stackrel{-}{P}}_D={\int }_0^{y_{\mathrm{eD}}}{P}_D\sin \left({\gamma }_n{y}_D\right){\mathrm{dy}}_D---\left(5\right)$

${\tilde{q}}_D={\int }_0^{x_{\mathrm{eD}}}{q}_D\left(x_D\right)\sin \left({\gamma }_n{x}_D\right){\mathrm{dx}}_D---\left(6\right)$

同时利用公式（2）（3）处理式（1），可以得到双重傅里叶（Fourier）积分变换下的压 力与裂缝流量的关系式：

${-\pi }^2(\frac{m^2}{x_{\mathrm{eD}}^2}+\frac{n^2}{y_{\mathrm{eD}}^2}){\stackrel{\widehat{‾}}{P}}_D+{\tilde{q}}_D\sin {\gamma }_n{y}_{\mathrm{wD}}=0---\left(7\right)$

利用反变换公式进行两次反演可以得到压力函数：

$P_D=\underset{n=1}{\overset{\infty }{\Sigma }}\frac{\sin \left({\gamma }_n{y}_D\right)}{N\left({\beta }_n\right)}\left(\underset{m=1}{\overset{\infty }{\Sigma }}\frac{\sin \left({\beta }_m{x}_D\right)}{N\left({\beta }_m\right)}{\stackrel{\widehat{‾}}{P}}_D\right)---\left(8\right)$

其中特征值：

β_m=mπ/x_eD；γ_n=nπ/y_eD                  （9）

范数倒数满足：

$N\left({\beta }_m\right)={\int }_0^{x_{\mathrm{eD}}}{\sin }^2\left({\beta }_m{x}_D\right){\mathrm{dx}}_D=\frac{x_{\mathrm{eD}}}{2};N\left({\gamma }_n\right)={\int }_0^{y_{\mathrm{eD}}}{\sin }^2\left({\gamma }_n{y}_D\right){\mathrm{dy}}_D=\frac{y_{\mathrm{eD}}}{2}---\left(10\right)$

将式（7）、式（9）和式（10）代入式（8）可以得到压力公式为：

$P_D=\underset{m=1}{\overset{\infty }{\Sigma }}\frac{2{\tilde{q}}_D}{x_{\mathrm{eD}}{y}_{\mathrm{eD}}}\sin \frac{{\mathrm{m\pi x}}_D}{x_{\mathrm{eD}}}\left[\underset{n=1}{\overset{\infty }{\Sigma }}\frac{\cos \mathrm{n\pi }(y_D-y_{\mathrm{wD}})/y_{\mathrm{eD}}-\cos \mathrm{n\pi }(y_D+y_{\mathrm{wD}})/y_{\mathrm{eD}}}{{\pi }^2(m^2/x_{\mathrm{eD}}^2+n^2/y_{\mathrm{eD}}^2)}\right]---\left(11\right)$

只有沿裂缝才有流量分布，所以式（6）可以改写为关于裂缝长度的积分关系式：

${\tilde{q}}_D={\int }_{x_{\mathrm{wD}}-1}^{x_{\mathrm{wD}}+1}{q}_D\left(\alpha \right)\sin \left(\frac{\mathrm{m\pi \alpha }}{x_{\mathrm{eD}}}\right)\mathrm{d\alpha }---\left(12\right)$

同时注意变换关系：

$\underset{k=1}{\overset{\infty }{\Sigma }}\frac{\cos \mathrm{k\pi x}}{k^2+a^2}=\frac{\pi }{2a}\frac{\cosh \left[\mathrm{a\pi }(1-x)\right]}{\sinh \left(\mathrm{a\pi }\right)}-\frac{1}{{2a}^2};[0\le x\le 2\pi ]---\left(13\right)$

利用式（12）、式（13）改写式（11），得到n条裂缝中的任一条引起的无量纲压力P_D与该条裂缝无量纲裂缝流量q_D(α)之间的关系式：

$\left(\begin{array}{c}{P}_D(x_D,y_D;x_{\mathrm{wD}},y_{\mathrm{wD}})\\ =2\underset{x_{\mathrm{wD}}-1}{\overset{x_{\mathrm{wD}}+1}{\int }}\left\{\underset{m=1}{\overset{\infty }{\Sigma }}\frac{q_D\left(\alpha \right)}{\mathrm{m\pi }}\sin \frac{{\mathrm{m\pi x}}_D}{x_{\mathrm{eD}}}\sin \frac{\mathrm{m\pi \alpha }}{x_{\mathrm{eD}}}\frac{\cosh \frac{\mathrm{m\pi }(y_{\mathrm{eD}}-|y_D-y_{\mathrm{wD}}\left|\right)}{x_{\mathrm{eD}}}-\cosh \frac{\mathrm{m\pi }(y_{\mathrm{eD}}-|y_D+y_{\mathrm{wD}}\left|\right)}{x_{\mathrm{eD}}}}{\sinh \frac{{\mathrm{m\pi y}}_{\mathrm{eD}}}{x_{\mathrm{eD}}}}\right\}\mathrm{d\alpha }\end{array}\right)---\left(14\right)$

式中：x_w，y_w为裂缝中心坐标，其通过物理模型的定义得到，其单位为米（m）。

在步骤C）中，基于物理模型，根据气体在裂缝中的流动关系以及裂缝与地层之间的边 界耦合关系，得到n条裂缝中的任一条在与井筒接触处产生的压力扰动P_wfn与标况下该条裂 缝流量的Q_scn的关系式：

$\left(\begin{array}{c}\frac{78.55\mathrm{kh}(P_i^2-P_{\mathrm{wfn}}^2)}{{\mathrm{\mu ZTQ}}_{\mathrm{scn}}}=\frac{1}{C_{\mathrm{fD}}}\frac{h}{x_f}[\ln \frac{h}{{2r}_w}-\frac{\pi }{2}]+f\left(C_{\mathrm{fD}}\right)\\ +2\left\{\underset{m=1}{\overset{\infty }{\Sigma }}\frac{x_e^2}{m^3{\pi }^2{x}_f^2}{\sin }^2\frac{{\mathrm{m\pi x}}_f}{x_e}{\sin }^2\frac{{\mathrm{m\pi x}}_w}{x_e}\frac{\cosh \frac{{\mathrm{m\pi y}}_e}{x_e}-\cosh \frac{\mathrm{m\pi }(y_e-{2y}_w)}{x_e}}{\sinh \frac{{\mathrm{m\pi y}}_e}{x_e}}\right\}\end{array}\right),$

式中：P_wfn为n条裂缝中的任一条在与井筒接触处产生的压力扰动，其单位为兆帕（MPa）。 r_w为井筒的半径，其通过技术手册得到，其单位为米（m）。

图2显示了气体在裂缝中渗流的示意图。如图2所示，气体从地层流入到裂缝后，通过 渗流作用运动至水平井筒。由于裂缝完全贯穿地层，可以认为裂缝内的流动为一维线性流动， 而且这种流动是一种沿裂缝长度变化的变质量流动。

由于裂缝体积较小，弹性较小，可以忽略裂缝弹性的影响，将流体在裂缝中的无量纲渗 流方程简化为稳态形式，得到式（15）和式（16）。

$\frac{d^2{P}_{\mathrm{fD}}}{{\mathrm{dx}}_D^2}+\frac{2}{C_{\mathrm{fD}}}{q}_D\left(x_D\right)=0,[-1\le x_D\le 1]---\left(15\right)$

$\frac{{\mathrm{dP}}_{\mathrm{fD}}\left(x_{\mathrm{wD}}\right)}{{\mathrm{dx}}_D}=-\frac{\pi }{C_{\mathrm{fD}}}---\left(16\right)$

其中，式中P_fD为由导流作用引起的压力扰动。

对式（15）关于x_D进行两次积分，可有：

$P_{\mathrm{wD}}-P_{\mathrm{fD}}\left(x_D\right)=\frac{\pi }{C_{\mathrm{fD}}}\left[\right|x_D-x_{\mathrm{wD}}|-{\int }_{x_{\mathrm{wD}}}^{x_D}\mathrm{dv}{\int }_{x_{\mathrm{wD}}}^v{q}_D\left(u\right)\mathrm{du}]---\left(17\right)$

因为压力是关于位置的函数，因此，在裂缝与地层的交界处的压力相同。裂缝与地层的 耦合条件为：

P_fD(x_D)=P_D(x_D,y_wD;x_wD,y_wD)，[-1≤x_D≤1]                   （18）

将式（14）代入式（18），形成Fredholm型积分方程，该方程无法解析求解，这里采用 数值解法：将裂缝等分为n份，等分段的流量、压力均匀，这样就会形成n+1阶变量为各段 流量q_Dj（j=1,2,3….,n）及井底压力P_wD的线性方程组，式（19）：

$\left(\begin{array}{c}{P}_{\mathrm{wD}}+2\underset{i=1}{\overset{n}{\Sigma }}{q}_{\mathrm{Di}}\underset{m=1}{\overset{\infty }{\Sigma }}\left(\begin{array}{c}\frac{x_{\mathrm{eD}}}{m^2{\pi }^2}\sin \mathrm{m\pi }\frac{x_{\mathrm{wD}}+(j-0.5)\mathrm{\Delta x}}{x_{\mathrm{eD}}}[\cos \mathrm{m\pi }\frac{x_{\mathrm{wD}}+\mathrm{i\Delta }{x}_D}{x_{\mathrm{eD}}}-\cos \mathrm{m\pi }\frac{x_{\mathrm{wD}}+(i-1)\Delta x_D}{x_{\mathrm{eD}}}\\ \times \frac{\cosh [\mathrm{m\pi }{y}_{\mathrm{eD}}/x_{\mathrm{eD}}]-\cosh [\mathrm{m\pi }(y_{\mathrm{eD}}-|2y_{\mathrm{wD}}\left|\right)/x_{\mathrm{eD}}]}{\sinh (\mathrm{m\pi }{y}_{\mathrm{eD}}/x_{\mathrm{eD}})}\end{array}\right)\\ =\frac{\pi }{C_{\mathrm{fD}}}\{(x_{\mathrm{wD}}+(j-0.5)\mathrm{\Delta x})(1-\underset{i=1}{\overset{j-1}{\Sigma }}{q}_{\mathrm{Di}}{\mathrm{\Delta x}}_D-\frac{{\mathrm{\Delta x}}_D}{2}{q}_{\mathrm{Dj}})+\underset{i=1}{\overset{j-1}{\Sigma }}{q}_{\mathrm{Di}}\Delta x_D[x_{\mathrm{wD}}+(i-0.5)\Delta x_D]+q_{\mathrm{Dj}}\Delta x_D\frac{x_{\mathrm{wD}}+(j-0.75)\Delta x_D}{2}\}\end{array}\right)---\left(19\right)$

流量约束方程

$\underset{i=1}{\overset{n}{\Sigma }}{q}_{\mathrm{Di}}=1---\left(20\right)$

利用牛顿迭代法求解线性方程组（19），并计算不同影响参数下导流能力C_fD与井底压力 P_wD的变化关系。得到的结果如图3所示，根据图3得知有限导流裂缝井底压力P_wD随C_fD增加 而减小。当C_fD>300（设定为300）时P_wD趋近于常数，即为无限导流裂缝对应的井底压力P_infwD。 而且这种变化趋势只与C_fD有关，不受其他参数的影响。图4显示了有限导流裂缝井底压力与 导流能力变化关系。如图4所示，通过数据回归可以得到无限导流裂缝与有限导流裂缝之间 的差值函数f(C_fD)：

$f\left(C_{\mathrm{fD}}\right)=\frac{1.65-0.328\ln C_{\mathrm{fD}}+0.116{\left(\ln C_{\mathrm{fD}}\right)}^2}{1.0+0.18\ln C_{\mathrm{fD}}+0.064{\left(\ln C_{\mathrm{fD}}\right)}^2+0.005{\left(\ln C_{\mathrm{fD}}\right)}^3}---\left(21\right)$

式（21）也是裂缝有限导流能力的影响函数，其中无限导流裂缝的井底压力为：

$P_{\mathrm{infwD}}=2\left\{\underset{n=1}{\overset{\infty }{\Sigma }}\frac{x_{\mathrm{eD}}}{n^2\pi }\sin \mathrm{n\pi }\frac{x_D}{x_{\mathrm{eD}}}\sin \mathrm{n\pi }\frac{1}{x_{\mathrm{eD}}}\sin \mathrm{n\pi }\frac{x_{\mathrm{wD}}}{x_{\mathrm{eD}}}\frac{\cosh \frac{{\mathrm{n\pi y}}_{\mathrm{eD}}}{x_{\mathrm{eD}}}-\cosh \frac{\mathrm{n\pi }(y_{\mathrm{eD}}-{2y}_{\mathrm{wD}})}{x_{\mathrm{eD}}}}{\sinh \frac{{\mathrm{n\pi y}}_{\mathrm{eD}}}{x_{\mathrm{eD}}}}\right\}---\left(22\right)$

图5显示了裂缝中的井筒聚流作用。如图5所示，同时在裂缝的近井筒附近会形成径向 流动，与垂直裂缝井相比会产生附加的压力降，即井筒的径向聚流效应。此处引入表皮因子 来加以考虑。

$\mathrm{skin}=\frac{\mathrm{kh}}{k_f{w}_f}[\ln \frac{h}{{2r}_w}-\frac{\pi }{2}]---\left(23\right)$

所以，考虑了井筒聚流作用的有限导流裂缝的井底压力P_finwD可以得到：

P_finwD=P_infwD+f(C_fD)+skin                         （24）

有量纲展开后的井底压力为：

$\left(\begin{array}{c}\frac{78.55\mathrm{kh}(P_i^2-P_{\mathrm{wfn}}^2)}{{\mathrm{\mu ZTQ}}_{\mathrm{scn}}}=\frac{1}{C_{\mathrm{fD}}}\frac{h}{x_f}[\ln \frac{h}{{2r}_w}-\frac{\pi }{2}]+f\left(C_{\mathrm{fD}}\right)\\ +2\left\{\underset{m=1}{\overset{\infty }{\Sigma }}\frac{x_e^2}{m^3{\pi }^2{x}_f^2}{\sin }^2\frac{{\mathrm{m\pi x}}_f}{x_e}{\sin }^2\frac{{\mathrm{m\pi x}}_w}{x_e}\frac{\cosh \frac{{\mathrm{m\pi y}}_e}{x_e}-\cosh \frac{\mathrm{m\pi }(y_e-{2y}_w)}{x_e}}{\sinh \frac{{\mathrm{m\pi y}}_e}{x_e}}\right\}\end{array}\right)---\left(25\right)$

式中：P_wfn为n条裂缝中的任一条在与井筒接触处产生的压力扰动，MPa。

式（25）得到n条裂缝中的任一条在与井筒接触处产生的压力扰动P_wfn与标况下该条裂 缝流量的Q_scn的关系式。

在步骤D）中，基于物理模型，根据气体在井筒中的流动关系以及裂缝与井筒之间的边 界耦合关系，得到n条裂缝中的任一条在与井筒接触处产生的压力扰动P_wfi、与该条裂缝相 邻的裂缝在与井筒接触处产生的压力扰动P_wfi-1和这两条裂缝在井筒内之间的流量Q_sci的关系 式。

图6显示了水平井筒中水平管流剖面图。如图6所示，由于井筒内径远远大于地层和裂缝 中的流动通道尺寸，因此井筒内的气体流动按单相管流进行计算。水平井筒内的流量是变化 的，因此，本文采用分段计算。忽略由于流速增大引起的动能压降，根据（李士伦,等.天 然气工程[M].北京:石油工业出版社,2008.Li S L,et al.Natural Gas Engineering[M]. Beijing:Petroleum Industry Press,2008.）中的内容，总压力梯度为：

$\frac{\mathrm{dP}}{\mathrm{dy}}=f\frac{{\mathrm{\rho v}}^2}{{2r}_w}---\left(26\right)$

采用平均参数法分离变量积分，得到n条裂缝中的任一条在与井筒接触处产生的压力扰 动P_wfi、与该条裂缝相邻的裂缝在与井筒接触处产生的压力扰动P_wfi-1和这两条裂缝在井筒内 之间的流量Q_sci的关系式：

$p_{\mathrm{wfi}}^2-P_{\mathrm{wfi}-1}^2=9\times {10}^{-12}\frac{{\mathrm{ZT\gamma }}_g}{r_w^5}{f}_i{d}_i{\left(\underset{j=i}{\overset{n}{\Sigma }}{Q}_{\mathrm{scj}}\right)}^2,$（其中，P_wf0=P_wf）            （27）

其中，e为井筒内壁粗糙度，其通过技术手册得到，其单位为毫米（mm）。Z为气体偏差 系数，其通过室内实验得到。γ_g为气体相对密度，其通过室内实验得到。f为摩阻系数。d 为裂缝间距，其通过压裂设计数据得到，其单位为米（m）。摩阻系数f由Jain提出的显式 公式计算，式中的雷诺数均考虑紊流状态，由（27）计算。

$f_i=[1.14-21g(\frac{e}{1000D}+\frac{21.25}{R_{\mathrm{ei}}^{0.9}}){]}^{-2}---\left(28\right)$

$R_{\mathrm{ei}}=\frac{{177.1\gamma }_g\underset{j=i}{\overset{n}{\Sigma }}{Q}_{\mathrm{scj}}}{{2\mu }_g{r}_w}---\left(29\right)$

在步骤E）中，估算第n条裂缝的流量Q_scn，测量得到第1条裂缝在与井筒接触处产生 的压力扰动P_wf，利用数值迭代法得到该水平井的产能。

E1）根据实际生产情况估算第n条裂缝流量的最大值Q_scn(max)和最小值Q_scn(min)，取最大值 与最小值的算术平均Q_scn(mid)=0.5×[Q_scn(max)+Q_scn(min)]；

E2）根据步骤C）中的公式计算得到对应的P_wfn(max)、P_wfn(min)和P_wfn(mid)，根据步骤D）中公 式得到对应的P_wfn(max)、P_wfn(min)和P_wfn(mid)，根据步骤C）中的公式计算得到对应的Q_scn-1(max)、Q_scn-1(min)、 Q_scn-1(mid)；

E3）重复步骤E2）得到P_wf0(max)、P_wf0(min)、P_wf0(mid)；

E4）将P_wf0(mid)与P_wf进行差值比较，如果P_wf0(mid)与P_wf的差值符合精度要求则判断计算值 P_wf0(mid)正确，如果P_wf0(mid)与P_wf的差值不符合精度要求，则如果(P_wf0(max)-P_wf)×(P_wf0(min)-P_wf)<0， 那么Q_scn(min)=Q_scn(mid)，否则Q_scn(max)=Q_scn(mid)；

E5）重复步骤E1）至E4），直至P_wf0(mid)与P_wf的差值符合精度要求。

在上述步骤中，P_wf0(mid)符合精度要求是指P_wf0(mid)与实际P_wf相同或者与P_wf的差值在允许的 范围内。P_wf0(mid)不符合精度要求是指P_wf0(mid)与P_wf的差值超出了允许的范围。

当P_wf0(mid)符合精度要求时，可得到该水平井的产能。

该水平井的产能为其中，j=1,2,3……，n。

当P_wf0(mid)符合精度要求时，还得到了不同裂缝参数和井底压力条件下的沿水平井筒压降 曲线、各条裂缝的流量和流量等产能指标。

与传统技术相比，本发明提出的经过多段压裂改造后的致密气藏水平井的产能计算方 法，考虑了地层-裂缝-井筒全过程的耦合流动产能理论预测方法，从理论上大大提高了致密 气藏多段压裂水平井产能预测的准确度。

本发明的又一实施例中，以苏里格气田某井为例进行计算分析。根据地质资料显示，该 井的控制面积约为1600m×600m（x_e*y_e），利用表1和表2气井的实际数据和本文推导的产能 公式，计算水平井产能。

表1苏里格气田某井的基本参数

表2苏里格气田某井经过人工压裂后的多段裂缝参数

当井底压力P_wf=0.1MPa时，考虑井筒摩阻各裂缝流量分别为Q_sc1=22.08×10⁴m³/d， Q_sc2=16.85×10⁴m³/，Q_sc3=2.18×10⁴m³/d，气井的产能为Q_AOF为41.12×10⁴m³/d现场利用压力恢 复测试数据和Topaze试井分析软件评价该井的产能为40.72×10⁴m³/d，相对误差为0.98%， 验证了模型及算法的正确性。

计算不同井底压力（0.1MPa，1MPa，5MPa，10MPa，20MPa）下的气井流量Q_sc并绘制图7， 各裂缝流量（表3）及沿井筒压力分布并绘制图8。图7显示了水平井的流量与井底压力变 化关系图。图8显示了不同井底流压下各裂缝沿水平井筒的压力分布。如图7所示，井筒摩 阻对气井流量Q_sc的影响随井底流压的减小而增加，当P_wf>25MPa时井筒摩阻影响可忽略。如 图8所示，各条裂缝沿井筒的压力分布均一，差异较小。表3反映了当考虑摩阻时，裂缝自 井筒的趾端到跟端流量不断增加。

图9显示了该井的井眼轨迹。从图9该井的井眼轨迹图上看，靠近趾端的储层几乎没有 实施加砂压裂，且井筒趾端钻遇了一段无效储层，因此井筒趾端部分几乎没有流量贡献。

表3不同井底流压下的各裂缝流量分布

本发明还公开了一种采用上述计算方法的产能计算系统，它包括建模单元、第一计算单 元、第二计算单元、第三计算单元和第四计算单元。

建模单元用于建立物理模型。该物理模型具有以下定义:A1)地层均质等厚。地层的俯视 面为矩形。该矩形具有四条闭合且等压的边界。所述矩形的宽为x_e，所述矩形的长为y_e。A2) 具有n条裂缝，所有n条裂缝完全贯穿地层，其中n=1,2,3……，第1条裂缝位于该水平井 的根端，第i条裂缝逐渐向该水平井的趾端排布,其中i=1,2,3……,n。

第一计算单元，其用于根据气体在地层中的渗流规律得到n条裂缝中的任一条引起的无 量纲压力P_D与该条裂缝无量纲裂缝流量q_D(α)之间的关系式，具体式可参见式（14）。

第二计算单元，其用于根据气体在裂缝中的流动关系以及裂缝与地层之间的边界耦合关 系得到n条裂缝中的任一条在与井筒接触处产生的压力扰动P_wfn与标况下该条裂缝流量的Q_scn的关系式，具体式可参见式（25）。

第三计算单元，其用于根据气体在裂缝中的流动关系以及裂缝与地层之间的边界耦合关 系得到n条裂缝中的任一条在与井筒接触处产生的压力扰动P_wfi、与该条裂缝相邻的裂缝在 与井筒接触处产生的压力扰动P_wfi-1和这两条裂缝在井筒内之间的流量Q_sci的关系式，具体式 可参见式（27）。

第四计算单元，其用于接收估算第n条裂缝的流量Q_scn，其用于接收由测量得到第1条 裂缝在与井筒接触处产生的压力扰动P_wf，其用于采用数值迭代法得到该水平井的产能，具 体步骤如步骤E）中所示。

本说明书中的各个实施例均采用递进的方式描述，每个实施例重点说明的都是与其他实 施例的不同之处，各个实施例之间相同相似的部分互相参见即可。

上述实施例只为说明本发明的技术构思及特点，其目的在于让熟悉此项技术的人士能够 了解本发明的内容并据以实施，并不能以此限制本发明的保护范围。凡根据本发明精神实质 所作的等效变化或修饰，都应涵盖在本发明的保护范围之内。