EI、Scopus 收录
中文核心期刊

三维压裂缝网不稳定压力半解析求解方法

王志凯, 程林松, 曹仁义, 王进, 贾品, 王选茹

王志凯, 程林松, 曹仁义, 王进, 贾品, 王选茹. 三维压裂缝网不稳定压力半解析求解方法. 力学学报, 2021, 53(8): 2246-2256. DOI: 10.6052/0459-1879-21-183
引用本文: 王志凯, 程林松, 曹仁义, 王进, 贾品, 王选茹. 三维压裂缝网不稳定压力半解析求解方法. 力学学报, 2021, 53(8): 2246-2256. DOI: 10.6052/0459-1879-21-183
Wang Zhikai, Cheng Linsong, Cao Renyi, Wang Jin, Jia Pin, Wang Xuanru. Semi-analytical method for unstable pressure of 3D fracture network. Chinese Journal of Theoretical and Applied Mechanics, 2021, 53(8): 2246-2256. DOI: 10.6052/0459-1879-21-183
Citation: Wang Zhikai, Cheng Linsong, Cao Renyi, Wang Jin, Jia Pin, Wang Xuanru. Semi-analytical method for unstable pressure of 3D fracture network. Chinese Journal of Theoretical and Applied Mechanics, 2021, 53(8): 2246-2256. DOI: 10.6052/0459-1879-21-183
王志凯, 程林松, 曹仁义, 王进, 贾品, 王选茹. 三维压裂缝网不稳定压力半解析求解方法. 力学学报, 2021, 53(8): 2246-2256. CSTR: 32045.14.0459-1879-21-183
引用本文: 王志凯, 程林松, 曹仁义, 王进, 贾品, 王选茹. 三维压裂缝网不稳定压力半解析求解方法. 力学学报, 2021, 53(8): 2246-2256. CSTR: 32045.14.0459-1879-21-183
Wang Zhikai, Cheng Linsong, Cao Renyi, Wang Jin, Jia Pin, Wang Xuanru. Semi-analytical method for unstable pressure of 3D fracture network. Chinese Journal of Theoretical and Applied Mechanics, 2021, 53(8): 2246-2256. CSTR: 32045.14.0459-1879-21-183
Citation: Wang Zhikai, Cheng Linsong, Cao Renyi, Wang Jin, Jia Pin, Wang Xuanru. Semi-analytical method for unstable pressure of 3D fracture network. Chinese Journal of Theoretical and Applied Mechanics, 2021, 53(8): 2246-2256. CSTR: 32045.14.0459-1879-21-183

三维压裂缝网不稳定压力半解析求解方法

基金项目: 国家自然科学基金资助项目(U1762210, 51774297), 中国石油科技项目重大项目(ZLZX2020-02-04), 中国石油天然气股份有限公司长庆油田分公司横向项目(ZY19-XA411-FW1128)资助
详细信息
    作者简介:

    程林松, 教授, 主要研究方向: 油气渗流理论与应用. E-mail: lscheng@cup.edu.cn

    曹仁义, 副教授, 主要研究方向: 油气渗流理论与应用. E-mail: caorenyi@cup.edu.cn

  • 中图分类号: TE33

SEMI-ANALYTICAL METHOD FOR UNSTABLE PRESSURE OF 3D FRACTURE NETWORK

  • 摘要: 受地应力及压裂工艺影响, 大斜度井水力压裂缝网展布复杂, 缝网中存在不同倾斜方向、不同展布形态及不同贯穿程度的压裂缝. 本文通过将裂缝面离散为若干矩形微元实现裂缝形态有效表征, 将渗流过程划分为基质向裂缝流动及裂缝向井筒流动两阶段, 采用有限差分方法构建离散裂缝面内不稳定渗流数值解, 结合封闭边界面源函数及叠加原理构建基质内不稳定渗流解析解, 耦合裂缝内流动数值解与基质内流动解析解, 求解了三维压裂缝网不稳定压力. 基于积分中值定理提出了点源、特殊线源代替面源求解基质内渗流的求解方法, 分析了该方法的可行性及适用条件, 在保证模型精度的同时提升了计算效率. 研究表明, 在基质内采用点源函数面积分求解面源的方法可准确求解三维压裂缝网井底压力动态但计算效率极低, 基于积分中值定理的点源、特殊线源近似面源求解方法可大大提升计算效率, 且在裂缝微元划分较为精细(微元无因次边长小于0.15)时可取得较高精度, 基于该模型分析了裂缝导流能力、裂缝倾角、裂缝高度及裂缝段间距对压裂大斜度井典型试井曲线的影响.
    Abstract: Due to the influence of in-situ stress and fracturing technology, the distribution of hydraulic fracture network in highly deviated wells is complex with different inclined directions, different distribution forms and different penetration degrees. In this paper, the fracture surface is discretized into several rectangular micro elements to realize the effective characterization of fracture morphology. The seepage process is divided into two stages: matrix flow to fracture and fracture flow to wellbore. The numerical solution of unsteady seepage in discrete fracture surface is constructed by using finite difference method, and the analytical solution of unsteady seepage in matrix is constructed by combining closed boundary source function and superposition principle. The solution of the unstable pressure of the 3D fracture network is obtained by coupling the numerical solution of flow in fracture and the analytical solution of flow in matrix. Based on the integral mean value theorem, a solution method of point source and special line source instead of surface source is proposed to solve the seepage in matrix. The feasibility and applicable conditions of this method are analyzed, which can ensure the accuracy of the model and improve the calculation efficiency. The research shows that the point source function area fraction method can accurately solve the bottom hole pressure dynamic of 3D pressure fracture network in the matrix, but the calculation is quite inefficient. The point source and special line source approximate surface source method based on the integral mean value theorem can greatly improve the calculation efficiency, and can achieve higher accuracy when the fracture micro element division is more precise (the dimensionless side length of micro element is less than 0.15). Based on this model, the flow regimes and related sensitivity analysis are analyzed. The results show that the fracture conductivity, fracture dip angle, fracture height and fracture interval have great influence on the typical well test curves of highly deviated wells. Fracture dip angle and fracture interval mainly affect the pressure and pressure derivative curve from linear flow to early radial flow, especially when the fracture height is small, the pressure derivative has backflow phenomenon.
  • 随着常规油气资源储量不断降低, 全球非常规油气资源勘探开发进入活跃期, 致密油、页岩油等非常规油气资源探明储量不断上升[1-3]. 特别是近5年, 中国非常规油气进入了规模发现与开发上产阶段, 致密、页岩油气展现出了巨大的资源潜力[4-7]. 与常规储层不同, 致密储层多为源储一体, 储集层致密, 未经压裂储层难以形成工业油流. 开发中需采用体积压裂技术将储层“打碎”, 形成“人工油气藏”[8,9]. 通常认为压裂缝沿最大主应力方向延伸[10-12], 但受原始地应力及压裂工艺影响, 压裂阶段储层应力状态非常复杂, 现场微地震数据及耦合应力场压裂模拟结果表明, 体积压裂后储层缝网展布复杂, 存在不同的倾角、形态及穿透比[13-16].

    对于上述储层, 其渗流规律较为复杂, 单一线性流难以表征其流动状态[17-19], 樊冬艳等[20]、Al Rbeawi等[21-22]基于点源函数叠加原理研究了相同倾斜方向下多段压裂无限导流倾斜缝典型试井曲线特征. 贾品等[23-25]、Teng和Li[26]、万义钊等[27]通过将裂缝离散为微元的方法实现不同裂缝形态的表征, 同时Teng和Li[28-29]、王苏冉[30]通过引入空间倾斜面源实现单条倾斜缝压力、产量瞬态分析. 但上述研究所考虑的裂缝只能沿单一方向倾斜, 在此条件下, 可通过调整坐标系的方法降低模型复杂度实现模型有效求解. 但唐慧莹等[14]研究表明, 即使在常规多段压裂过程中, 压裂缝也会因缝间应力干扰发生裂缝面弯曲, 出现裂缝倾斜方向不同的现象. 因此为更准确表征水力压裂缝空间展布, 建立可表征任意方向倾斜、任意形态展布的多段压裂井压力瞬态分析模型[31-34]是非常必要的.

    本文通过离散裂缝网络表征不同裂缝形态, 结合点源函数面积分与裂缝微元局部坐标系, 构建了不受空间坐标系约束的封闭盒式储层任意平面缝源函数, 并验证了其准确性. 针对三重积分面源积分函数存在的计算效率较低的问题, 提出了应用点源或特殊线源代替面源的方法实现了倾斜裂缝快速求解, 探讨了此类求解方法的准确性及可行性, 基于此模型划分了单段、多段倾斜压裂缝流态, 分析了裂缝导流能力、裂缝倾角、裂缝高度以及裂缝段间距等因素的影响.

    图1所示, 均质封闭储层 (xe, ye, h) 中存在一口体积压裂大斜度井, 沿井筒延伸方向存在M条倾斜缝, 裂缝半缝长为$ {l}_{\text{f/2}} $, m; 沿裂缝面延伸高度为$ {h}_{\text{f}} $, m; 裂缝宽度为w, m; 各裂缝倾斜方向不同, 每条裂缝内采用局部坐标$ ({e_{{\rm{\varepsilon I}}}},{e_{{\rm{\xi I}}}})$进行定位.

    图  1  多段压裂井三维缝物理模型
    Figure  1.  Physical model of multi-fractured well with 3D fractures

    储层及裂缝孔隙度为$ {\phi }_{\mathrm{m}},\; {\phi }_{\mathrm{f}} $; 渗透率为km, kf, mD; 综合压缩系数为Ctm, Ctf, MPa−1. 流体黏度为μ, mPa·s; 体积系数为B; 油井定产qw生产, m3/d; 井储系数为C, m3/MPa; 不考虑表皮影响.

    模型基本假设如下:

    (1) 将渗流阶段划分为基质向裂缝渗流及裂缝向井筒流动两部分, 不考虑基质向井筒的直接渗流;

    (2) 不考虑井筒沿程压力降, 忽略渗流过程中的重力影响;

    (3) 裂缝有限导流, 单条裂缝内孔渗数据保持一致;

    (4) 储层内流体为油单相, 无毛管力影响.

    水力裂缝形态受储层物性及压裂工艺影响[35], 常见裂缝形态有矩形缝、椭圆缝、菱形缝等. 本文将裂缝离散为矩形微元, 实现不同形态裂缝有效表征. 如图2所示, 井筒穿过椭圆缝中心, 结合椭圆裂缝形态将其离散为17个矩形微元, 其中中心点微元(编号9)为井筒所在微元.

    图  2  椭圆缝离散表征
    Figure  2.  Discrete characterization of elliptical fracture

    模型推导与计算过程中, 裂缝微元相关参数定义如下: 第I条裂缝微元总量为NfI, 其中井筒所在微元编号NfwI; 裂缝系统微元总量为Nf; 第i个微元的长、宽为$\Delta {\varepsilon }_{{i}}$, $\Delta {\xi }_{{i}}$, m; 定义每个裂缝微元的渗透率为Ki, 相应的可确定裂缝微元有限导流能力为$ {K}_{i}{w}_{i} $, 通过给定不同的Ki值即可实现裂缝不同有限导流能力求解.

    基于裂缝内二维不稳定渗流,建立裂缝内流动综合控制方程.

    $$\begin{split} &\frac{\partial }{{\partial \varepsilon }}\left( {\beta \frac{{\Delta {\xi _{{i}}}w{k_{\rm{f}}}}}{{\mu B}}\frac{{\partial p}}{{\partial \varepsilon }}} \right)\Delta {\varepsilon _{{i}}} +\frac{\partial }{{\partial \xi }}\left( {\beta \frac{{\Delta {\varepsilon_{{i}}}w{k_{\rm{f}}}}}{{\mu B}}\frac{{\partial p}}{{\partial \xi }}} \right)\Delta {\xi _{{i}}} + \\&\qquad{q_{{\rm{s}}{{\rm{c}}_{{i}}}}} = \left( {\frac{{\Delta {\varepsilon _{{i}}}\Delta {\xi _{{i}}}w{\phi _{\rm{f}}}{C_{{\rm{tf}}}}}}{B}} \right)\frac{{\partial p}}{{\partial t}} \end{split} $$ (1)

    考虑井筒与椭圆缝相交于微元9, 以该微元为研究对象, 采用有限差分方法将公式(1)离散处理

    $$\begin{split} &\frac{{\beta \Delta {\xi _9}w{k_{\rm{f}}}}}{{\mu B{d_{(8,9)}}}}\left( {p_{{\rm{f9}}}^n - p_{{\rm{f8}}}^n} \right) + \frac{{\beta \Delta {\xi _9}w{k_{\rm{f}}}}}{{\mu B{d_{(9,10)}}}}\left( {p_{{\rm{f9}}}^n - p_{{\rm{f10}}}^n} \right) +\\ &\qquad\frac{{\beta \Delta {\varepsilon _9}w{k_{\rm{f}}}}}{{\mu B\Delta {d_{(3,9)}}}}\left( {p_{{\rm{f9}}}^n - p_{{\rm{f3}}}^n} \right) + \frac{{\beta \Delta {\varepsilon _9}w{k_{\rm{f}}}}}{{\mu B{d_{(9,15)}}}}\left( {p_{{\rm{f9}}}^n - p_{{\rm{f15}}}^n} \right) + \\ &\qquad(q_{{\rm{f9}}}^n - q_{{\rm{fw1}}}^n) = \frac{{\Delta {\varepsilon _9}\Delta {\xi _9}w{\phi _{\rm{f}}}{C_{{\rm{tf}}}}}}{{B\Delta t}}\left( {p_{{\rm{f9}}}^{n - 1} - p_{{\rm{f9}}}^n} \right)\end{split} $$ (2)

    式中, $\;\beta$为单位转换系数, 取值0.0864; $ {q}_{\text{fi}}^{{n}} $为第n时间步基质向第i个微元的供给, m3/d; $ {q}_{\text{f-w}\text{I}}^{n} $为第n时间步第i条裂缝向井筒的 供给, m3/d; $ {p}_{\text{fi}}^{n} $为第n时间步第i个微元的压力, MPa; $ {\Delta d}_{(i,j)} $为微元i, j中心点连线的距离, m.

    为简化模型推导过程, 引入无量纲参数

    ${x_{\rm{D}}} = \dfrac{x}{{{l_{{\rm{f/2}}}}}}$, ${y_{\rm{D}}} = \dfrac{y}{{{l_{{\rm{f/2}}}}}}$, ${{\textit{z}}_{\rm{D}}} = \dfrac{{\textit{z}}}{{{l_{{\rm{f/2}}}}}}$, ${w_{\rm{D}}} = \dfrac{w}{{{l_{{\rm{f/2}}}}}}$, ${h_{\rm{D}}} = \dfrac{h}{{{l_{{\rm{f/2}}}}}}$, ${x_{{\rm{eD}}}} = \dfrac{{{x_e}}}{{{l_{{\rm{f/2}}}}}}$, ${y_{{\rm{eD}}}} = \dfrac{{{y_e}}}{{{l_{{\rm{f/2}}}}}}$, $\Delta {l_{\rm{D}}} = \dfrac{{\Delta \varepsilon }}{{{l_{{\rm{f/2}}}}}}$, $\Delta {d_{{{{\rm{D}}(i,j)}}}} = \dfrac{{\Delta {d_{{{(i,j)}}}}}}{{{l_{{\rm{f/2}}}}}}$, $\Delta {\xi _{\rm{D}}} = \dfrac{{\Delta \xi }}{{{l_{{\rm{f/2}}}}}}$, ${q_{\rm{D}}} = \dfrac{q}{{{q_{\rm{w}}}}}$, ${C_{{\rm{fD}}}} = \dfrac{{{k_{\rm{f}}}w}}{{{k_m}{l_{{\rm{f/2}}}}}}$, ${C_{\rm{s}}} = \dfrac{{{\phi _{\rm{f}}}{C_{{\rm{tf}}}}}}{{{\phi _{\rm{m}}}{C_{{\rm{tm}}}}}}$, ${p_{\rm{D}}} = $$ \dfrac{{2{\text{π}} \beta {k_m}h\Delta p}}{{{q_{\rm{w}}}B\mu }}$, ${p_{{\rm{fD}}}} = \dfrac{{2{\text{π}} \beta {k_m}h\Delta {p_{\rm{f}}}}}{{{q_{\rm{w}}}B\mu }}$, ${t_{\rm{D}}} = \dfrac{{\beta {k_m}t}}{{{\phi _{\rm{m}}}\mu {C_{{\rm{tm}}}}l_{{\rm{f/2}}}^2}}$, ${C_{\rm{D}}} = $$ \dfrac{C}{{2{\text{π}} h{\phi _{\rm{m}}}{C_{{\rm{tm}}}}r_w^2}}$, $\gamma = \dfrac{{0.0434Br_w^2}}{{l_{{\rm{f/2}}}^2}}$, ${h_{{\rm{fD}}}} = \dfrac{{{h_{\rm{f}}}}}{{l_{{\rm{f/2}}}^{}}}$.

    将无量纲参数代入式(2)

    $$\begin{split} & \frac{{{C_{{\rm{fD}}}}\Delta {\xi _{{\rm{D9}}}}}}{{2{\text{π}} {h_{\rm{D}}}{d_{{\rm{D(8,9)}}}}}}\left( {p_{{\rm{fD}}9}^n - p_{{\rm{fD}}8}^n} \right) + \frac{{{C_{{\rm{fD}}}}\Delta {\xi _{{\rm{D9}}}}}}{{2{\text{π}} {h_{\rm{D}}}{d_{{\rm{D(9,10)}}}}}}\left( {p_{{\rm{fD}}9}^n - p_{{\rm{fD}}10}^n} \right) + \\ &\qquad\frac{{{C_{{\rm{fD}}}}\Delta {\varepsilon _{{\rm{D9}}}}}}{{2{\text{π}} {h_{\rm{D}}}{d_{{\rm{D(3,9)}}}}}}\left( {p_{{\rm{fD}}9}^n \!-\! p_{{\rm{fD}}3}^n} \right) \!+\! \frac{{{C_{{\rm{fD}}}}\Delta {\varepsilon _{{\rm{D9}}}}}}{{2{\text{π}} {h_{\rm{D}}}{d_{{\rm{D(9,15)}}}}}}(p_{{\rm{fD}}9}^n \!-\! p_{{\rm{fD}}15}^n) + \\ &\qquad (q_{{\rm{fD9}}}^n -q_{{\rm{fwD1}}}^n) = \frac{{\Delta {\varepsilon _{{\rm{D9}}}}\Delta {\xi _{{\rm{D9}}}}{w_{\rm{D}}}{C_{\rm{s}}}}}{{2{\text{π}} \Delta{t_{\rm{D}}}{h_{\rm{D}}}}}(p_{{\rm{fD9}}}^{n - 1} - p_{{\rm{fD9}}}^n) \end{split} $$ (3)

    化简可得

    $$\begin{split} & \sum\limits_{{{j = 1}}}^{{N_{\rm{f}}}} {{T_{{{{\rm{D}}(i,j)}}}}} p_{{{{\rm{fD}}i}}}^{{n}} - \left( {\sum\limits_{{{j}} = {\rm{1}}}^{{N_{\rm{f}}}} {{T_{{{{\rm{D}}(i,j)}}}}} + {\alpha _{{\rm{D9}}}}} \right)p_{{\rm{fD9}}}^{{n}} -\\ &\qquad (q_{{\rm{fD9}}}^n - q_{{\rm{fwD1}}}^n) = - {\alpha _{{\rm{D9}}}}p_{{\rm{fD9}}}^{{{n - 1}}} \\ \end{split} $$ (4)

    式中, TD(i,j)i, j两裂缝微元间的无因次传导率, $ {\alpha }_{{{\rm{D}}i}} $为相邻时间步间的影响.

    $${T_{{{{\rm{D}}(i,j)}}}} = \frac{{{C_{{\rm{fD}}}}\Delta {n_{{{{\rm{D}}(i,j)}}}}}}{{2{\text{π}} {h_{\rm{D}}}{d_{{\rm{D}}({{i,j}})}}}}\qquad$$ (5)
    $${\alpha _{{{{\rm{D}}i}}}} = \frac{{\Delta {\varepsilon _{{{{\rm{D}}i}}}}\Delta {\xi _{{{{\rm{D}}i}}}}{w_{\rm{D}}}{C_{\rm{s}}}}}{{2{\text{π}} \Delta {t_{\rm{D}}}{h_{\rm{D}}}}}\;\;$$ (6)
    $$\Delta {n_{{{{\rm{D}}(i,j)}}}} = \left\{ {\begin{array}{*{20}{l}} {\Delta {\xi _{{{{\rm{D}}i}}}}},&{i ,j\;\;{\text{在}} \varepsilon{\text{方向上相邻}} }\\ {\Delta {\varepsilon _{{{{\rm{D}}i}}}}},&{i ,j\;\;{\text{在}} \varepsilon{\text{方向上相邻}}}\\ 0,&{i ,j\;\;{\text{不相邻}}} \end{array}} \right.$$ (7)

    将式(4)应用到离散裂缝系统所有微元上, 可得包含Nf个裂缝微元的流动方程有限差分方程组

    $$\left[ {{\boldsymbol{T}}, - {{\boldsymbol{I}}_{{N_{\rm{f}}}}},{\boldsymbol{0}},{\boldsymbol{a}}} \right] \cdot \left[ {\begin{array}{*{20}{c}} {{{\boldsymbol{p}}_{{\rm{fD}}}}} \\ {{{\boldsymbol{q}}_{{\rm{fD}}}}} \\ {{{\boldsymbol{p}}_{{\rm{fwD}}}}} \\ {{{\boldsymbol{q}}_{{\rm{fwD}}}}} \end{array}} \right] = {{\boldsymbol{R}}_1}$$ (8)

    式中, 矩阵$ \boldsymbol{T} $大小为Nf × Nf, T(i, j) = TD(i,j), 表征微元间不稳定渗流影响; ${\boldsymbol{I}}_{{{N}}_{\rm{f}}}$Nf × Nf的单位矩阵, 表征基质源汇项影响; 矩阵0大小为Nf × M的零矩阵; 矩阵a大小为Nf × M, 表征井筒源汇项影响; 矩阵R1大小为Nf × 1, 由第n−1时间步参数决定

    $${\boldsymbol{a}} \!= \!{\left[ {\begin{array}{*{20}{l}} {\left. {\begin{array}{*{20}{c}} 0 \\ \vdots \\ 0 \\ 1 \end{array}} \right\}{{{N}}_{{\rm{FW1}}}}}&\!\!\!{}&\!\!\!{\left. {\begin{array}{*{20}{c}} 0 \\ \vdots \\ 0 \\ 1 \end{array}} \right\}{{{N}}_{{{{\rm{FW}}i}}}}}&\!\!\!{}&\!\!\!{\left. {\begin{array}{*{20}{c}} 0 \\ \vdots \\ 0 \\ 1 \end{array}} \right\}{N_{{{{\rm{FW}}M}}}}} \\ 0&\!\!\! \cdots &\!\!\!0&\!\!\! \cdots &\!\!\!0 \\ \vdots &\!\!\!{}&\!\!\! \vdots &\!\!\!{}&\!\!\! \vdots \\ 0&\!\!\!{}&\!\!\!0&\!\!\!{}&\!\!\!0 \!\!\!\!\! \end{array}} \right]_{{N_{\rm{f}}} \times {{M}}}}$$ (9)
    $${{\boldsymbol{p}}_{{\rm{fD}}}} \!=\! {\left[ {\begin{array}{*{20}{c}} {{\boldsymbol{p}}_{{\rm{fD1}}}^n} \\ {{\boldsymbol{p}}_{{\rm{fD2}}}^n} \\ \vdots \\ {{\boldsymbol{p}}_{{\rm{fD}}{{\rm{N}}_{\rm{f}}}}^n} \!\!\!\! \end{array}} \right]_{{N_{\rm{f}}} \times 1}},{{\boldsymbol{q}}_{{\rm{fD}}}} \!=\! {\left[ {\begin{array}{*{20}{c}} {{\boldsymbol{q}}_{{\rm{fD1}}}^n} \\ {{\boldsymbol{q}}_{{\rm{fD2}}}^n} \\ \vdots \\ {{\boldsymbol{q}}_{{\rm{fD}}{{\rm{N}}_{\rm{f}}}}^n} \!\!\!\! \end{array}} \right]_{{N_{\rm{f}}} \times 1}},{{\boldsymbol{p}}_{{\rm{fw}}}} \!=\! {\left[ {\begin{array}{*{20}{c}} {{\boldsymbol{p}}_{{\rm{fwD1}}}^n} \\ {{\boldsymbol{p}}_{{\rm{fwD2}}}^n} \\ \vdots \\ {{\boldsymbol{p}}_{{\rm{fwDM}}}^n} \!\!\!\! \end{array}} \right]_{{{M}} \times 1}}$$ (10)
    $${{\boldsymbol{q}}_{{\rm{fw}}}} = {\left[ {\begin{array}{*{20}{c}} {{\boldsymbol{q}}_{{\rm{fwD1}}}^n} \\ {{\boldsymbol{q}}_{{\rm{fwD2}}}^n} \\ \vdots \\ {{\boldsymbol{q}}_{{\rm{fwDM}}}^n} \end{array}} \right]_{{\boldsymbol{M }}\times 1}},{{\boldsymbol{R}}_1} = {\left[ {\begin{array}{*{20}{c}} { - {{\boldsymbol{\alpha }}_{{\rm{D}}1}}{\boldsymbol{p}}_{{\rm{fD1}}}^{n - 1}} \\ { - {{\boldsymbol{\alpha }}_{{\rm{D}}2}}{\boldsymbol{p}}_{{\rm{fD2}}}^{n - 1}} \\ \vdots \\ { - {{\boldsymbol{\alpha}} _{{\rm{D}}{N_{\rm{f}}}}}{\boldsymbol{p}}_{{\rm{fD}}{{\rm{N}}_{\rm{f}}}}^{n - 1}} \end{array}} \right]_{{{{N}}_{\rm{f}}} \times 1}}$$ (11)

    考虑上述离散方法所得第I个微元中心点坐标为$ ({x}_{\text{DI}},{y}_{\text{DI}},{{\textit{z}}}_{\text{DI}}) $, 其局部坐标系由矩形微元相邻边所在方向的单位向量$({{\boldsymbol{e}}_{{\rm{\varepsilon I}}}},{{\boldsymbol{e}}_{{\rm{\xi I}}}})$构成, 基于中心点坐标及局部坐标系可确定微元内任意点坐标

    $$\left. {\begin{array}{*{20}{c}} {{{\boldsymbol{e}}_{{\rm{\varepsilon I}}}} = ({x_{{\rm{\varepsilon I}}}},{y_{{\rm{\varepsilon I}}}},{{\textit{z}}_{{\rm{\varepsilon I}}}})} \\ {{{\boldsymbol{e}}_{\xi {\rm{I}}}} = ({x_{\xi {\rm{I}}}},{y_{\xi {\rm{I}}}},{{\textit{z}}_{\xi {\rm{I}}}})} \end{array}} \right\}$$ (12)
    $$\left. {\begin{array}{*{20}{c}} {{x_{\rm{D}}}(\varepsilon,\xi ) = {x_{{\rm{DI}}}} + \varepsilon \cdot {x_{\varepsilon {\rm{I}}}} + \xi \cdot {x_{\xi {\rm{I}}}}} \\ {{y_{\rm{D}}}(\varepsilon,\xi ) = {y_{{\rm{DI}}}} + \varepsilon \cdot {y_{\varepsilon {\rm{I}}}} + \xi \cdot {y_{\xi {\rm{I}}}}} \\ {{{\textit{z}}_{\rm{D}}}(\varepsilon,\xi ) = {{\textit{z}}_{{\rm{DI}}}} + \varepsilon \cdot {{\textit{z}}_{\varepsilon {\rm{I}}}} + \xi \cdot {{\textit{z}}_{\xi {\rm{I}}}}} \end{array}} \right\}$$ (13)

    结合Newman乘积与封闭边界无限大平面源经典解[36]可得封闭储层(x, y, z)处单一点源以定产q生产时在点(x0, y0, z0)处产生压力降

    $$\begin{split} {p_{\rm{e}}} - &{p_{\rm{0}}} = \frac{{Bq}}{{{\phi _{\rm{m}}}{c_{{\rm{tm}}}}}} \cdot \\ & \frac{1}{{{x_e}}}\left[ {1 + 2\sum\limits_{m = 1}^\infty {\exp } \left( { - \frac{{{m^2}{{\text{π}} ^2}\beta \eta t}}{{x_e^2}}} \right)\cos\left( \frac{{m{\text{π}} {x_0}}}{{{x_e}}}\right)\cos \left( \frac{{m{\text{π}} x}}{{{x_e}}}\right)} \right] \cdot \\ & \frac{1}{{{y_e}}}\left[ {1 + 2\sum\limits_{m = 1}^\infty {\exp } \left( { - \frac{{{m^2}{{\text{π}} ^2}\beta \eta t}}{{y_e^2}}} \right)\cos\left( \frac{{m{\text{π}} {y_0}}}{{{y_e}}}\right)\cos \left(\frac{{m{\text{π}} y}}{{{y_e}}}\right)} \right] \cdot \\ & \frac{1}{h}\left[ {1 + 2\sum\limits_{m = 1}^\infty {\exp } \left( { - \frac{{{m^2}{{\text{π}} ^2}\beta \eta t}}{{{h^2}}}} \right)\cos \left(\frac{{m{\text{π}} {{\textit{z}}_0}}}{h}\right)\cos\left( \frac{{m{\text{π}} {\textit{z}}}}{h}\right)} \right] \\ \end{split} $$ (14)

    将无量纲参数代入式(14)

    $$\begin{split} {p_{{\rm{D0}}}} = &\frac{{2{\text{π}} {q_{\rm{D}}}}}{{{x_{{\rm{eD}}}}{y_{{\rm{eD}}}}}} \cdot \\ & \left[ {1 + 2\sum\limits_{m = 1}^\infty {\exp } \left( { - \frac{{{m^2}{{\text{π}} ^2}{t_{\rm{D}}}}}{{x_{{\rm{eD}}}^2}}} \right)\cos \left(\frac{{m{\text{π}} {x_{{\rm{D0}}}}}}{{{x_{{\rm{eD}}}}}} \right)\cos \left( \frac{{m{\text{π}} {x_{\rm{D}}}}}{{{x_{{\rm{eD}}}}}}\right)} \right] \cdot \\ & \left[ {1 + 2\sum\limits_{m = 1}^\infty {\exp } \left( { - \frac{{{m^2}{{\text{π}} ^2}{t_{\rm{D}}}}}{{y_{{\rm{eD}}}^2}}} \right)\cos \left(\frac{{m{\text{π}} {y_{{\rm{D0}}}}}}{{{y_{{\rm{eD}}}}}}\right)\cos \left( \frac{{m{\text{π}} {y_{\rm{D}}}}}{{{y_{{\rm{eD}}}}}}\right)} \right] \cdot \\ & \left[ {1 + 2\sum\limits_{m = 1}^\infty {\exp } \left( { - \frac{{{m^2}{{\text{π}} ^2}{t_{\rm{D}}}}}{{h_{\rm{D}}^2}}} \right)\cos \left(\frac{{m{\text{π}} {{\textit{z}}_{{\rm{D0}}}}}}{{{h_{\rm{D}}}}}\right)\cos \left( \frac{{m{\text{π}} {{\textit{z}}_{\rm{D}}}}}{{{h_{\rm{D}}}}}\right)} \right] \\ \end{split} $$ (15)

    由叠加原理可知, 空间面源可由空间点源沿微元面积分获得, 因此空间面源函数表达式

    $$\begin{split} &{p_{{\rm{D0}}}} = \dfrac{{2{\text{π}} {q_{\rm{D}}}}}{{{x_{{\rm{eD}}}}{y_{{\rm{eD}}}}\Delta {\varepsilon _{\rm{D}}}\Delta {\xi _{\rm{D}}}}} \cdot \displaystyle\int_{\varepsilon = - \frac{{\Delta {\varepsilon _{\rm{D}}}}}{2}}^{\varepsilon = \frac{{\Delta {\varepsilon _{\rm{D}}}}}{2}} {\displaystyle\int_{\xi = - \frac{{\Delta {\xi _{\rm{D}}}}}{2}}^{\xi = \frac{{\Delta {\xi _{\rm{D}}}}}{2}} {} } \\ & \left[ {1 + 2\sum\limits_{m = 1}^\infty {\exp } \left( { - \dfrac{{{m^2}{{\text{π}} ^2}{t_{\rm{D}}}}}{{x_{{\rm{eD}}}^2}}} \right)\cos \left(\dfrac{{m{\text{π}} {x_{{\rm{D0}}}}}}{{{x_{{\rm{eD}}}}}}\right)\cos \left(\dfrac{{m{\text{π}} {x_{\rm{D}}}(\varepsilon,\xi )}}{{{x_{{\rm{eD}}}}}}\right)} \right] \cdot \\ & \left[ {1 + 2\sum\limits_{m = 1}^\infty {\exp } \left( { - \dfrac{{{m^2}{{\text{π}} ^2}{t_{\rm{D}}}}}{{y_{{\rm{eD}}}^2}}} \right)\cos \left(\dfrac{{m{\text{π}} {y_{{\rm{D0}}}}}}{{{y_{{\rm{eD}}}}}}\right)\cos \left(\dfrac{{m{\text{π}} {y_{\rm{D}}}(\varepsilon,\xi )}}{{{y_{{\rm{eD}}}}}}\right)} \right] \cdot \\ & \left[ {1 \!\!+\!\! 2\sum\limits_{m = 1}^\infty {\exp } \left( { - \dfrac{{{m^2}{{\text{π}} ^2}{t_{\rm{D}}}}}{{h_{\rm{D}}^2}}} \right)\cos \left(\dfrac{{m{\text{π}} {{\textit{z}}_{{\rm{D0}}}}}}{{{h_{\rm{D}}}}}\right)\cos \left(\dfrac{{m{\text{π}} {{\textit{z}}_{\rm{D}}}(\varepsilon,\xi )}}{{{h_{\rm{D}}}}}\right)} \right]{\rm{d}}\varepsilon {\rm{d}}\xi \\ \end{split} $$ (16)

    基于空间面源压力分布, 结合压力叠加原理及杜哈美原理得第i个微元第n时间步在第j个微元处产生的压力降

    $$\begin{split} & p_{{{{\rm{D}}j}}}^n = \frac{{2{\text{π}} }}{{\Delta {\varepsilon _{\rm{D}}}\Delta {\xi _{\rm{D}}}{x_{{\rm{eD}}}}{y_{{\rm{eD}}}}}}\sum\limits_{i = 1}^{{N_{\rm{f}}}} {\sum\limits_{k = 1}^n {q_{{{{\rm{D}}i}}}^k\int_{t_{\rm{D}}^{{{k - 1}}}}^{t_{\rm{D}}^{{k}}} {\int_{\varepsilon = - \frac{{\Delta {\varepsilon _{\rm{D}}}}}{2}}^{\varepsilon = \frac{{\Delta {\varepsilon _{\rm{D}}}}}{2}} {\int_{\xi = - \frac{{\Delta {\xi _{\rm{D}}}}}{2}}^{\xi = \frac{{\Delta {\xi _{\rm{D}}}}}{2}} {} } } } } \cdot \\ & \left\{ \!\!{1 \!\!+\!\! 2\sum\limits_{m = 1}^\infty {\!\!\exp \!\!} \left[ { - \frac{{{m^2}{{\text{π}} ^2}({t_{\rm{D}}} - \tau )}}{{x_{{\rm{eD}}}^2}}} \right]\!\cos\!\left( \frac{{m{\text{π}} {x_{{\rm{Dj}}}}}}{{{x_{{\rm{eD}}}}}}\right)\!\cos\! \left( \frac{{m{\text{π}} {x_{{\rm{Di}}}}(\varepsilon,\xi )}}{{{x_{{\rm{eD}}}}}}\right)} \right\} \cdot \\ & \left\{ \!\!{1 \!\!+\!\! 2\sum\limits_{m = 1}^\infty {\!\!\exp \!\!} \left[ { - \frac{{{m^2}{{\text{π}} ^2}({t_{\rm{D}}} - \tau )}}{{y_{{\rm{eD}}}^2}}} \right]\!\cos\! \left(\frac{{m{\text{π}} {y_{{\rm{Dj}}}}}}{{{y_{{\rm{eD}}}}}}\right)\!\cos\! \left(\frac{{m{\text{π}} {y_{{\rm{Di}}}}(\varepsilon,\xi )}}{{{y_{{\rm{eD}}}}}}\right)} \right\} \cdot \\ & \left\{ \!\!{1 \!\!+\!\! 2\sum\limits_{m = 1}^\infty {\!\!\exp \!\!} \left[ { - \frac{{{m^2}{{\text{π}} ^2}({t_{\rm{D}}} - \tau )}}{{h_{\rm{D}}^2}}} \right]\!\cos\! \left(\frac{{m{\text{π}} {z_{{\rm{Dj}}}}}}{{{h_{\rm{D}}}}}\right)\!\cos\! \left(\frac{{m{\text{π}} {z_{{\rm{Di}}}}(\varepsilon,\xi )}}{{{h_{\rm{D}}}}}\right)} \right\}\\ & {\rm{d}}\varepsilon {\rm{d}}\xi {\rm{d}}\tau \end{split} $$ (17)

    为简化方程形式, 定义式(16)中的积分项为参数$ {G}_{i,j}^{n,k} $

    $$\begin{split} & G_{i,j}^{n,k} = \frac{{2{\text{π}} }}{{\Delta {\varepsilon _{\rm{D}}}\Delta {\xi _{\rm{D}}}{x_{{\rm{eD}}}}{y_{{\rm{eD}}}}}}\int_{t_{\rm{D}}^{{{k - 1}}}}^{t_{\rm{D}}^{{k}}} {\int_{\varepsilon = - \frac{{\Delta {\varepsilon _{\rm{D}}}}}{2}}^{\varepsilon = \frac{{\Delta {\varepsilon _{\rm{D}}}}}{2}} {\int_{\xi = - \frac{{\Delta {\xi _{\rm{D}}}}}{2}}^{\xi = \frac{{\Delta {\xi _{\rm{D}}}}}{2}} {} } } \\ & \left\{\!\!{1 \!\!+\!\! 2\sum\limits_{m = 1}^\infty {\!\!\exp \!\!} \left[ { - \frac{{{m^2}{{\text{π}} ^2}({t_{\rm{D}}} - \tau )}}{{x_{{\rm{eD}}}^2}}} \right]\!\cos\! \left( \frac{{m{\text{π}} {x_{{\rm{Dj}}}}}}{{{x_{{\rm{eD}}}}}}\right)\!\cos\! \left( \frac{{m{\text{π}} {x_{{\rm{Di}}}}(\varepsilon,\xi )}}{{{x_{{\rm{eD}}}}}}\right)} \right\} \cdot \\ & \left\{ \!\!{1 \!\!+\!\! 2\sum\limits_{m = 1}^\infty {\!\!\exp \!\!} \left[ { - \frac{{{m^2}{{\text{π}} ^2}({t_{\rm{D}}} - \tau )}}{{y_{{\rm{eD}}}^2}}} \right]\!\cos\! \left( \frac{{m{\text{π}} {y_{{\rm{Dj}}}}}}{{{y_{{\rm{eD}}}}}}\right)\!\cos\! \left( \frac{{m{\text{π}} {y_{{\rm{Di}}}}(\varepsilon,\xi )}}{{{y_{{\rm{eD}}}}}}\right)} \right\} \cdot \\ & \left\{ \!\!{1 \!\!+\!\! 2\sum\limits_{m = 1}^\infty {\!\!\exp \!\!} \left[ { - \frac{{{m^2}{{\text{π}} ^2}({t_{\rm{D}}} - \tau )}}{{h_{\rm{D}}^2}}} \right]\!\cos\! \left( \frac{{m{\text{π}} {z_{{\rm{Dj}}}}}}{{{h_{\rm{D}}}}}\right)\!\cos\! \left( \frac{{m{\text{π}} {z_{{\rm{Di}}}}(\varepsilon,\xi )}}{{{h_{\rm{D}}}}}\right)} \right\}\\ &{\rm{d}}\varepsilon {\rm{d}}\xi {\rm{d}}\tau \end{split} $$ (18)

    从而

    $$p_{{{{\rm{D}}j}}}^n = \sum\limits_{i = 1}^{{N_{\rm{f}}}} {\sum\limits_{k = 1}^n {q_{{{{\rm{D}}i}}}^kG_{i,j}^{n,k}} } = q_{{{{\rm{D}}i}}}^k\sum\limits_{i = 1}^{{N_{\rm{f}}}} {G_{i,j}^{n,k}} + \sum\limits_{i = 1}^{{N_{\rm{f}}}} {\sum\limits_{k = 1}^{n - 1} {q_{{{{\rm{D}}i}}}^kG_{i,j}^{n,k}} } $$ (19)

    基于式(19)建立基质向裂缝流动矩阵

    $${\boldsymbol{B}} \cdot \left[ {\begin{array}{*{20}{c}} {{{\boldsymbol{p}}_{{\rm{fD}}}}} \\ {{{\boldsymbol{q}}_{{\rm{fD}}}}} \\ {{{\boldsymbol{p}}_{{\rm{fwD}}}}} \\ {{{\boldsymbol{q}}_{{\rm{fwD}}}}} \end{array}} \right] = {{\boldsymbol{R}}_2}$$ (20)

    式中, 矩阵B大小为Nf × 2(Nf + M), 表征基质向各微元流动的影响; 矩阵R2大小为Nf × 1, 由n − 1时间步中压力及流量参数构成

    $${\boldsymbol{B}} \!\!=\!\! {\left[ {\begin{array}{*{20}{c}} {\underbrace {1{\rm{ }}0 \cdots 0}_{{N_{\rm{f}}}}}&{\underbrace { - {\boldsymbol{G}}_{1,1}^{n,n} \cdots - {\boldsymbol{G}}_{{N_{\rm{f}}},1}^{n,n}}_{{N_{\rm{f}}}}}&{\underbrace {0 \cdots 0}_{2M}} \\ {}& \vdots &{} \\ {\underbrace {\overbrace {0 \cdots 0}^{g - 1}{\rm{ }}1{\rm{ }}0 \cdots 0}_{{N_{\rm{f}}}}}&{\underbrace { - {\boldsymbol{G}}_{1,g}^{n,n} \cdots - {\boldsymbol{G}}_{{N_{\rm{f}}},g}^{n,n}}_{{N_{\rm{f}}}}}&{\underbrace {0 \cdots 0}_{2M}} \\ {}& \vdots &{} \\ {\underbrace {0 \cdots 0{\rm{ }}1}_{{N_{\rm{f}}}}}&{\underbrace { - {\boldsymbol{G}}_{1,{N_{\rm{f}}}}^{n,n} \cdots - {\boldsymbol{G}}_{{N_{\rm{f}}},{N_{\rm{f}}}}^{n,n}}_{{N_{\rm{f}}}}}&{\underbrace {0 \cdots 0}_{2M}} \end{array}} \right]_{{N_{\rm{f}}} \!\times\! 2({N_{\rm{f}}} \!\!+\!\! M)}}$$ (21)
    $${{\boldsymbol{R}}_2} = {\left[ {\begin{array}{*{20}{c}} {\displaystyle\sum\limits_{k = 1}^{n - 1} {\displaystyle\sum\limits_{i = 1}^{{N_{\rm{f}}}} {{\boldsymbol{q}}_{{\rm{fi}}}^k{\boldsymbol{G}}_{i,1}^{n,k}} } } \\ {\displaystyle\sum\limits_{k = 1}^{n - 1} {\displaystyle\sum\limits_{i = 1}^{{N_{\rm{f}}}} {{\boldsymbol{q}}_{{\rm{fi}}}^k{\boldsymbol{G}}_{i,2}^{n,k}} } } \\ \vdots \\ {\displaystyle\sum\limits_{k = 1}^{n - 1} {\displaystyle\sum\limits_{i = 1}^{{N_{\rm{f}}}} {{\boldsymbol{q}}_{{\rm{fi}}}^k{\boldsymbol{G}}_{i,{N_{\rm{f}}}}^{n,k}} } } \end{array}} \right]_{{N_{\rm{f}}} \times 1}}$$ (22)

    由式(19)可知, 随着Nfn的增大, 式(18)的调用次数迅速升高, 而式(18)中的三重积分会导致参数$ {G}_{i,j}^{n,k} $求解过程中计算效率较低, 因此考虑以下两种求解方法.

    (1) 线源代替面源

    当多条裂缝法向量共面时, 可建立坐标轴(x, y, z)使其中一条坐标轴与裂缝微元局部坐标系$ \left( {{{\boldsymbol{e}}_{{\rm{\varepsilon I}}}},{{\boldsymbol{e}}_{{\rm{\xi I}}}}} \right) $的一条向量平行, 实现模型简化, 如图3所示.

    图  3  裂缝微元局部坐标系示意图
    Figure  3.  Schematic diagram of local coordinate system of fracture elements

    以局部坐标轴eε与坐标轴x平行为例, 此时xξ,$ {y}_{\varepsilon },\;{{\textit{z}}}_{\varepsilon } $均为0, 式(12)中$ {y}_{{\rm{D}}}(\varepsilon,\xi ) $, $ {{\textit{z}}}_{{\rm{D}}}(\varepsilon,\xi ) $将转化成ξ的函数

    $$\left. {\begin{array}{*{20}{c}} {{x_{\rm{D}}}(\varepsilon ) = {x_{{\rm{DI}}}} + \varepsilon {e_x}} \\ {{y_{\rm{D}}}(\xi ) = {y_{{\rm{DI}}}} + \xi {y_{\xi {\rm{I}}}}} \\ {{{\textit{z}}_{\rm{D}}}(\xi ) = {{\textit{z}}_{{\rm{DI}}}} + \xi {{\textit{z}}_{\xi {\rm{I}}}}} \end{array}} \right\}$$ (23)

    相应的参数$ {G}_{i,j}^{n,k} $可简化为双重积分形式

    $$\begin{split} & {\rm{G}}_{i,j}^{n,k} = \frac{{2{\text{π}} }}{{\Delta {\xi _{\rm{D}}}{x_{{\rm{eD}}}}{y_{{\rm{eD}}}}}}\int_{t_{\rm{D}}^{{{k - 1}}}}^{t_{\rm{D}}^{{k}}} {\int_{\xi = - \frac{{\Delta {\xi _{\rm{D}}}}}{2}}^{\xi = \frac{{\Delta {\xi _{\rm{D}}}}}{2}} {\left\{ {1 + \frac{{4{x_{{\rm{eD}}}}}}{{{\text{π}} \Delta {x_{\rm{D}}}}}} \right.} } \cdot \sum\limits_{m = 1}^\infty {\frac{1}{m}} \cdot\\ & \left. {\!\!\exp \!\! \left[ { \!\!-\! \frac{{{m^2}{{\text{π}} ^2}(t_{\rm{D}}^n - {\tau _{\rm{D}}})}}{{x_{{\rm{eD}}}^2}}} \right]\sin \left(\frac{{m{\text{π}} \Delta {x_{\rm{D}}}}}{{2{x_{{\rm{eD}}}}}}\right)\!\cos\! \left( \frac{{m{\text{π}} {x_{{\rm{D}}i}}}}{{{x_{{\rm{eD}}}}}}\right)\!\cos\! \left(\frac{{m{\text{π}} {x_{{\rm{DI}}}}}}{{{x_{\rm{e}}}}}\right)} \right\} \cdot \\ & \left\{ {1 + 2\sum\limits_{m = 1}^\infty {\!\!\exp \!\! } \left[ { - \frac{{{m^2}{{\text{π}} ^2}({t_{\rm{D}}} - \tau )}}{{y_{{\rm{eD}}}^2}}} \right]\!\cos\! \left(\frac{{m{\text{π}} {y_{{\rm{Dj}}}}}}{{{y_{{\rm{eD}}}}}}\right)\!\cos\! \left(\frac{{m{\text{π}} {y_{{\rm{Di}}}}(\xi )}}{{{y_{{\rm{eD}}}}}}\right)} \right\} \cdot \\ & \left\{ {1 \!\!+\!\! 2\sum\limits_{m = 1}^\infty {\!\!\exp \!\! } \left[{ - \frac{{{m^2}{{\text{π}} ^2}({t_{\rm{D}}} - \tau )}}{{h_{\rm{D}}^2}}} \right]\!\!\cos\!\left( \frac{{m{\text{π}} {{\textit{z}}_{{\rm{Dj}}}}}}{{{h_{\rm{D}}}}}\right)\!\!\cos\! \left(\frac{{m{\text{π}} {{\textit{z}}_{{\rm{Di}}}}(\xi )}}{{{h_{\rm{D}}}}}\right)} \right\}\\ &{\rm{d}}\xi {\rm{d}}\tau \end{split} $$ (24)

    考虑参数$ {G}_{i,j}^{n,k} $由线源沿微元斜边积分所得, 通过调整ξ方向上微元数量使单个微元内该方向上线源函数变化较小, 引入积分中值定理, 以微元中心点$ ({x}_{\text{DI}},{y}_{\text{DI}},{{\textit{z}}}_{\text{DI}}) $处线源与Δξ之积近似代替源函数积分实现式(24)的进一步简化, 获得参数$ {G}_{i,j}^{n,k} $的一重积分形式

    $$\begin{split} & {{G}}_{i,j}^{n,k} = \frac{{2{\text{π}} }}{{{x_{{\rm{eD}}}}{y_{{\rm{eD}}}}}}\int_{t_{\rm{D}}^{{{k - 1}}}}^{t_{\rm{D}}^{{k}}} {\left\{ {1 + \frac{{4{x_{{\rm{eD}}}}}}{{{\text{π}} \Delta {x_{\rm{D}}}}}} \right.} \cdot \sum\limits_{m = 1}^\infty {\frac{1}{m} \cdot } \\ & \exp\!\! \left[ { \!\!- \frac{{{m^2}{{\text{π}} ^2}(t_{\rm{D}}^n \!-\! {\tau _{\rm{D}}})}}{{x_{{\rm{eD}}}^2}}} \right]\left. {\!\!\!\sin\! \left( \frac{{m{\text{π}} \Delta {x_{\rm{D}}}}}{{2{x_{{\rm{eD}}}}}}\!\!\right)\!\cos\! \left( \frac{{m{\text{π}} {x_{Di}}}}{{{x_{{\rm{eD}}}}}}\right)\!\cos\! \left( \frac{{m{\text{π}} {x_{{\rm{DI}}}}}}{{{x_e}}}\right)} \right\} \cdot \\ & \left\{ {1 \!\!+\!\! 2\sum\limits_{m = 1}^\infty {\exp } \left[ { - \frac{{{m^2}{{\text{π}} ^2}({t_{\rm{D}}} - \tau )}}{{y_{{\rm{eD}}}^2}}} \right]\!\cos\! \left( \frac{{m{\text{π}} {y_{{\rm{Dj}}}}}}{{{y_{{\rm{eD}}}}}}\right)\!\cos\!\left( \frac{{m{\text{π}} {y_{{\rm{Di}}}}}}{{{y_{{\rm{eD}}}}}}\right)} \right\} \cdot \\ & \left\{ {1 \!\!+\!\! 2\sum\limits_{m = 1}^\infty {\exp } \left[ { - \frac{{{m^2}{{\text{π}} ^2}({t_{\rm{D}}} - \tau )}}{{h_{\rm{D}}^2}}} \right]\!\cos\!\left( \frac{{m{\text{π}} {{\textit{z}}_{{\rm{Dj}}}}}}{{{h_{\rm{D}}}}}\right)\!\cos\!\left( \frac{{m{\text{π}} {{\textit{z}}_{{\rm{Di}}}}}}{{{h_{\rm{D}}}}}\right)} \right\}{\rm{d}}\tau \\ \end{split} $$ (25)

    (2) 点源代替面源

    参照方案一的思路, 直接对三重积分形式下参数$ {G}_{i,j}^{n,k} $使用积分中值定理, 以微元中心点$ ({x}_{\text{DI}},{y}_{\text{DI}}, $$ {{\textit{z}}}_{\text{DI}}) $处点源与微元面积Δε × Δξ之积近似代替源函数积分实现式(18)的简化, 该条件下需保证裂缝微元的尺寸在一定范围内, 后面将对其精度进行讨论

    $$\begin{split} & G_{i,j}^{n,k} = \frac{{2{\text{π}} }}{{{x_{{\rm{eD}}}}{y_{{\rm{eD}}}}}}\int_{t_{\rm{D}}^{{{k - 1}}}}^{t_{\rm{D}}^{{k}}} {} \\ & \left\{ {1 + 2\sum\limits_{m = 1}^\infty {\exp } \left[ { - \frac{{{m^2}{{\text{π}} ^2}({t_{\rm{D}}} - \tau )}}{{x_{{\rm{eD}}}^2}}} \right]\cos \left(\frac{{m{\text{π}} {x_{{\rm{Dj}}}}}}{{{x_{{\rm{eD}}}}}}\right)\cos \left( \frac{{m{\text{π}} {x_{{\rm{Di}}}}}}{{{x_{{\rm{eD}}}}}}\right)} \right\} \cdot \\ & \left\{ {1 + 2\sum\limits_{m = 1}^\infty {\exp } \left[ { - \frac{{{m^2}{{\text{π}} ^2}({t_{\rm{D}}} - \tau )}}{{y_{{\rm{eD}}}^2}}} \right]\cos \left(\frac{{m{\text{π}} {y_{{\rm{Dj}}}}}}{{{y_{{\rm{eD}}}}}}\right)\cos \left(\frac{{m{\text{π}} {y_{{\rm{Di}}}}}}{{{y_{{\rm{eD}}}}}}\right)} \right\} \cdot \\ & \left\{ {1 \!+\! 2\sum\limits_{m = 1}^\infty {\exp } \left[ { - \frac{{{m^2}{{\text{π}} ^2}({t_{\rm{D}}} - \tau )}}{{h_{\rm{D}}^2}}} \right]\cos\left( \frac{{m{\text{π}} {{\textit{z}}_{{\rm{Dj}}}}}}{{{h_{\rm{D}}}}}\right)\cos\left( \frac{{m{\text{π}} {{\textit{z}}_{{\rm{Di}}}}}}{{{h_{\rm{D}}}}}\right)} \right\}{\rm{d}}\tau \\ \end{split} $$ (26)

    将上述方案中的$ {G}_{i,j}^{n,k} $代入式(21)可实现简化求解条件下油藏流动模型求解矩阵的构建.

    由Peaceman公式[37]可知井筒所在微元内边界条件可由稳态径向流公式处理

    $$p_{{\rm{fwD}}}^n - p_{{\rm{fD}}}^n = \frac{{{w_{\rm{D}}}q_{{\rm{fwD}}}^n\ln ({r_{{\rm{eqD}}}}{\rm{/}}{r_{{\rm{wD}}}})}}{{{C_{{\rm{fD}}}}}}$$ (27)

    式中${r}_{\text{eqD}}=0.14\sqrt{{\left(\Delta {l}_{{\rm{D}}}\right)}^{2}+{\left(\Delta {m}_{{\rm{D}}}\right)}^{2}}$, 表征该微元的等效半径.

    基于Hurst等人[38]的研究, 在模型中考虑井储系数CwD的影响.

    $$1 - \frac{{\gamma {C_{{\rm{wD}}}}}}{{\Delta t_{\rm{D}}^n}}\left( {p_{{\rm{wD}}}^n - p_{{\rm{wD}}}^{n - 1}} \right) = \sum\limits_{i = 1}^M {q_{{{{\rm{fwD}}i}}}^n} $$ (28)

    式中$\gamma = \dfrac{{0.043\;4Br_w^2}}{{l_{{\rm{f/2}}}^2}}$, 为无量纲系数.

    同时, 井筒无限导流条件下, 各裂缝处井筒压力相等, 即

    $$p_{{{{\rm{fwD}}i}}}^n = p_{{{{\rm{fwD}}(i + 1)}}}^n,\;\;i \in (2,3, \cdots, M)$$ (29)

    结合式(27)−(29)建立与井筒内流动相关渗流方程组

    $${\boldsymbol{C}}\left[ {\begin{array}{*{20}{c}} {{{\boldsymbol{p}}_{{\rm{fD}}}}} \\ {{{\boldsymbol{q}}_{{\rm{fD}}}}} \\ {{{\boldsymbol{p}}_{{\rm{fwD}}}}} \\ {{{\boldsymbol{q}}_{{\rm{fwD}}}}} \end{array}} \right] = {{\boldsymbol{R}}_3}$$ (30)

    式中, 矩阵C大小为2M × 2(Nf + M), 表征了式(27)−(29)中第n时间步压力及产量的系数项, 矩阵R3大小为Nf × 1, 由n−1时间步相关参数构成.

    $ {\boldsymbol{A}}=\left[T,-{I}_{{\mathrm{N}}_{\mathrm{f}}},o,a\right] $, 联立方程组(8), (20), (30)即可获得多段压裂大斜度井耦合求解模型

    $$\left[ {\begin{array}{*{20}{c}} {\boldsymbol{A}} \\ {\boldsymbol{B}} \\ {\boldsymbol{C}} \end{array}} \right] \cdot \left[ {\begin{array}{*{20}{c}} {{{\boldsymbol{p}}_{{\rm{fD}}}}} \\ {{{\boldsymbol{q}}_{{\rm{fD}}}}} \\ {{{\boldsymbol{p}}_{{\rm{fwD}}}}} \\ {{{\boldsymbol{q}}_{{\rm{fwD}}}}} \end{array}} \right] = \left[ {\begin{array}{*{20}{c}} {{{\boldsymbol{R}}_1}} \\ {{{\boldsymbol{R}}_2}} \\ {{{\boldsymbol{R}}_3}} \end{array}} \right]$$ (31)

    该方程组中存在2(Nf + M)个方程, 其中未知量包含Nf个裂缝微元的压力$ {p}_{\text{fD}} $, Nf个裂缝微元的产量$ {q}_{\text{fD}} $, M条裂缝处的井底压力$ {p}_{\text{fwD}} $以及M条裂缝向井筒的供给$ {q}_{\text{fwD}} $共计2(Nf + M)个, 未知量数与方程数相等, 方程组可封闭求解.

    本模型中, 考虑多段压裂大斜度井压裂缝为有限导流, 且裂缝面沿任意方向展布, 现有文献相关研究较少. 为验证模型准确性, 分别验证有限导流单一倾斜缝与无限导流多段压裂缝求解的准确性. 针对图4所示部分贯穿倾斜缝(partially penetrating inclined fracture, PPIF), Teng和Li[28]通过建立坐标轴1 (Coordinate 1)实现模型点源函数求解简化. 为验证本文模型三重积分求解的准确性, 将坐标轴1沿z轴旋转45°建立新坐标轴2, 增加模型复杂度.

    图  4  坐标轴转换示意图
    Figure  4.  Schematic diagram of coordinate axis transformation

    模型数据与Teng图4数据相同, 具体无因次参数如下: xeD = yeD = 20, hD = 0.5, θ = 45°, hfD = $ 0.25\sqrt{2} $, CfD = 5, Cs = 10, CwD = 0, wD = 1.0 × 10−5, γ = 1.3 × 10−8.

    坐标轴2条件下微元内局部坐标系如下

    $$\left. {\begin{array}{*{20}{c}} {{{\boldsymbol{e}}_{{\rm{\varepsilon I}}}} = \left( {\dfrac{{\sqrt 2 }}{2}, - \dfrac{{\sqrt 2 }}{2},0} \right)} \\ {{{\boldsymbol{e}}_{\xi {\rm{I}}}} = \left( {\dfrac{{\sqrt 6 }}{6},\dfrac{{\sqrt 6 }}{6}, - \dfrac{{\sqrt 6 }}{3}} \right)} \end{array}} \right\}$$ (32)

    图5所示, 本模型所得不同时间下无因次压力及压力导数曲线与Teng模型结果基本吻合.

    图  5  本文模型与Teng模型[28]结果对比
    Figure  5.  Comparison of the result of our model with Teng model[28]

    为研究多缝条件下模型求解准确性, 选取Al Rbeawi和Tiab[22]2013年提出的无限大地层多段压裂倾斜缝图版(A-1)进行对比验证, 选取裂缝倾角为15°, hD = 1, 无因次段间距DD = 1. 如图6所示.

    图  6  本文模型与Al Rbeawi模型结果对比
    Figure  6.  Comparison of the result of our model with Al Rbeawi model

    可以看出, 不同时间下两模型计算结果基本吻合, 由于本模型采用封闭边界, 当tD > 10压力响应到达封闭边界后的压力及压力导数变化与Al Rbeawi模型不同由于该模型采用三重积分进行源函数计算, 求解时间较长(数小时), 适用性较差, 考虑对点源函数求解进行进一步调整.

    线源、点源简化模型求解主要误差来自于积分中值的选取, 当裂缝微元沿积分方向的长度足够小时, 点源函数沿积分方向上单调变化, 此时选取积分中点作为积分中值可取得较好的积分效果. 图7为不同网格划分下线源求解所得曲线变化对比.可以看出, 当网格划分较为粗糙时, 线源法所得初期压力及压力导数曲线误差较大. 特别的压力导数曲线存在较为明显的回流特征, 考虑是由于中心线所在线源与积分中值相差较大, 各线源之间存在流向各微元中心线源的径向流动所致(图8).

    图  7  网格划分对线源求解效果影响
    Figure  7.  Effect of meshing on line source solution
    图  8  回流段流动示意图
    Figure  8.  Diagram of flow in flowback period

    随着划分微元数目增多, 简化模型所得压力及压力导数曲线与面源积分所得压力及压力导数曲线趋于一致. 大量模拟表明当所划分微元无因次边长在0.15左右时, 采用点源或线源代替面源的方法可获得较高的精度. 选取hfD = $ 0.25\sqrt{2} $的单缝, 划分微元数目15 × 3, 模拟获得相应典型试井曲线如图9所示. 可以看出在tD大于10−7时间内采用点源或线源函数代替面源积分可取得较高计算精度, 相同模型下, 面源积分计算时长为435 443 s, 而线源积分计算时长仅为17 s, 点源积分计算时长仅为22 s, 极大地提升了模型计算效率.

    图  9  简化模型求解效果示意图
    Figure  9.  Diagram of simplified model

    通过分析典型试井曲线中的压力及压力导数曲线特征, 可实现不同流动阶段识别. 选取无因次参数xeD = yeD = 20, hD = 0.5, θ = 45°, hfD = $ 0.25\sqrt{2} $, CfD = 50, Cs = 10, CwD = 5, wD = 1.0 × 10−5, γ = 1.3 × 10−8. 绘制考虑井储的压力及压力导数曲线. 从图10可以看出, 单一倾斜缝典型试井曲线可划分为8个阶段.

    图  10  单一倾斜缝流动阶段划分
    Figure  10.  Identification of the flow regimes of incline fracture

    (1)井筒储集及过渡阶段: 井储阶段压力及压力导数曲线重合, 均为斜率等于1的直线, 该阶段主要流体供给来自井筒;

    (2)裂缝线性流阶段: 压力导数曲线斜率为1/2, 该阶段基质暂未动用, 与井相连的裂缝内部流体近线性流入井筒;

    (3)双线性流阶段: 压力导数曲线斜率为1/4, 该阶段除裂缝内部流动外, 近缝基质向裂缝的线性流动开始出现;

    (4)基质线性流阶段: 压力导数曲线斜率为1/2, 该阶段裂缝内压力趋于稳定, 近缝基质向裂缝的流动仍以线性流为主;

    (5)早期径向流阶段: 压力导数曲线斜率为0, 此时波及范围扩大, 基质向裂缝的供给由线性流转变为径向流;

    (6)椭圆流阶段: 压力导数曲线斜率近似为0.36, 通常认为该阶段为缝周基质径向流向裂缝系统径向流过渡的阶段;

    (7)晚期径向流阶段: 压力导数为斜率等于0的直线, 流体围绕裂缝整体成径向流;

    (8)边界控制流阶段: 流体流动到封闭边界, 压力及压力导数响应主要受边界影响, 压力及压力导数均为斜率等于1的直线.

    为研究裂缝条数对流动阶段影响, 取缝间距DD = 1, 对比裂缝条数分别为1, 2, 3条件下压力及压力导数曲线(图11)可以看出, 随着裂缝条数增加, 主要流动阶段仍为8个, 但与单条缝相比存在一定的差别. 随着裂缝段数增多, 裂缝及近缝区域流动影响增强, 井储阶段提前进入过渡阶段, 相应的双线性流阶段出现时间略有提前; 由于缝间干扰的存在, 裂缝条数越多, 早期径向流阶段越不太明显, 过渡阶段斜率更大, 晚期径向流及边界控制流阶段, 流体相对于裂缝系统进行流动, 不同裂缝条数下压力及压力导数曲线基本重合.

    图  11  多段压裂倾斜缝流动阶段
    Figure  11.  Flow regimes of multi-stage incline fractures

    本文所建模型在处理裂缝过程中主要进行了以下几点考虑: 多段压裂缝、裂缝有限导流、裂缝以任意形态展布以及裂缝向任意方向倾斜, 因此针对上述几点考虑分别进行裂缝导流能力、裂缝倾角、井筒倾角、裂缝纵向高度以及裂缝段间距等参数敏感性分析.

    (1)裂缝导流能力敏感性分析

    图12给出了无因次裂缝导流系数分别为0.5, 5, 50及500时压力及压力导数曲线响应, 可以看出随着裂缝导流能力的升高, 井储阶段及过渡阶段更早结束, 早期及晚期径向流阶段时间变短, 特别是早期径向流阶段基本消失.

    图  12  裂缝导流能力敏感性分析
    Figure  12.  Sensitivity analysis of well storage coefficient

    (2)裂缝倾角敏感性分析

    考虑裂缝纵向高度保持一致, 对比裂缝面倾角分别为15, 30, 45, 60以及75时压力及压力导数曲线变化, 由图13可以看出, 裂缝面倾角主要影响双线性流阶段至椭圆流阶段的压力导数曲线, 当裂缝倾角小于45°时, 倾角影响较小, 随着倾角进一步增大, 双线性流作用时间增长, 早期径向流更为明显, 相应的椭圆流阶段倾角增大.

    图  13  裂缝倾角敏感性分析
    Figure  13.  Sensitivity analysis of well storage coefficient

    (3)裂缝高度敏感性分析

    取裂缝倾角为45°, 考虑裂缝无因次纵向高度在0.2–1.0范围内变化. 由图14可以看出, 当裂缝高度较小时, 裂缝内线性流即双线性流阶段较不明显, 且基质线性流后期存在一定的回流现象, 随着裂缝高度的升高, 随着裂缝高度的增大, 各阶段流态逐渐出现.

    图  14  裂缝高度敏感性分析
    Figure  14.  Sensitivity analysis of fracture height

    (4)裂缝段间距敏感性分析

    取裂缝倾角为15°, 考虑无因次裂缝段间距在1–4范围内变化, 由图15可以看出, 随着裂缝段间距增大, 早期径向流至晚期径向流阶段的特征更为不明显, 表明随着裂缝段间距的增大, 缝间干扰越晚出现且系统径向流出现时间越晚.

    图  15  段间距敏感性分析
    Figure  15.  Sensitivity analysis of fracturing interval

    (1)对于任意方向展布裂缝可通过裂缝离散表征, 耦合裂缝内数值解与基质内面源函数解析解的方法进行求解, 但存在计算效率低的问题. 借助点源、特殊线源替代面源的方法可大大提升计算效率, 且在裂缝微元划分较为精细(微元无因次边长小于0.15)时精度较高.

    (2)倾斜缝典型试井曲线可划分为井筒储集及过渡流、裂缝线性流、双线性流、基质线性流、早期径向流、椭圆流、晚期径向流及边界控制流八个流动阶段. 随着裂缝段数增多, 缝间影响增强, 井储、过渡阶段及径向流阶段特征减弱.

    (3)裂缝倾角及裂缝段间距主要影响裂缝线性流到早期径向流阶段压力及压力导数曲线, 特别的当裂缝高度较小时压力导数存在回流现象.

  • 图  1   多段压裂井三维缝物理模型

    Figure  1.   Physical model of multi-fractured well with 3D fractures

    图  2   椭圆缝离散表征

    Figure  2.   Discrete characterization of elliptical fracture

    图  3   裂缝微元局部坐标系示意图

    Figure  3.   Schematic diagram of local coordinate system of fracture elements

    图  4   坐标轴转换示意图

    Figure  4.   Schematic diagram of coordinate axis transformation

    图  5   本文模型与Teng模型[28]结果对比

    Figure  5.   Comparison of the result of our model with Teng model[28]

    图  6   本文模型与Al Rbeawi模型结果对比

    Figure  6.   Comparison of the result of our model with Al Rbeawi model

    图  7   网格划分对线源求解效果影响

    Figure  7.   Effect of meshing on line source solution

    图  8   回流段流动示意图

    Figure  8.   Diagram of flow in flowback period

    图  9   简化模型求解效果示意图

    Figure  9.   Diagram of simplified model

    图  10   单一倾斜缝流动阶段划分

    Figure  10.   Identification of the flow regimes of incline fracture

    图  11   多段压裂倾斜缝流动阶段

    Figure  11.   Flow regimes of multi-stage incline fractures

    图  12   裂缝导流能力敏感性分析

    Figure  12.   Sensitivity analysis of well storage coefficient

    图  13   裂缝倾角敏感性分析

    Figure  13.   Sensitivity analysis of well storage coefficient

    图  14   裂缝高度敏感性分析

    Figure  14.   Sensitivity analysis of fracture height

    图  15   段间距敏感性分析

    Figure  15.   Sensitivity analysis of fracturing interval

  • [1] 张理, 李志川. 潮流能开发现状、发展趋势及面临的力学问题. 力学学报, 2016, 48(5): 1019-1032 (Zhang Li, Li Zhichuan. Development status, trend and the problems of mechanics of tidal current energy. Chinese Journal of Theoretical and Applied Mechanics, 2016, 48(5): 1019-1032 (in Chinese) doi: 10.6052/0459-1879-15-299
    [2]

    Conti JJ, Holtberg PD, Beamon JA, et al. Annual energy outlook 2020. US Energy Information Administration, 2020, 2

    [3] 康玉柱. 中国非常规油气勘探重大进展和资源潜力. 石油科技论坛, 2018, 37(4): 1-7 (Kang Yuzhu. Significant exploration progress and resource potential of unconventional oil and gas in China. Petroleum Science and Technology Forum, 2018, 37(4): 1-7 (in Chinese) doi: 10.3969/j.issn.1002-302x.2018.04.001
    [4] 孙龙德, 邹才能, 贾爱林等. 中国致密油气发展特征与方向. 石油勘探与开发, 2019, 46(6): 1015-1026 (Sun Longde, Zou Caineng, Jia Ailin, et al. Development characteristics and orientation of tight oil and gas in China. Petroleum Exploration and Development, 2019, 46(6): 1015-1026 (in Chinese)
    [5] 李国欣, 朱如凯. 中国石油非常规油气发展现状、挑战与关注问题. 中国石油勘探, 2020, 25(2): 1-13 (Li Guoxin, Zhu Rukai. Progress, challenges and key issues of unconventional oil and gas development of CNPC. China Petroleum Exploration, 2020, 25(2): 1-13 (in Chinese) doi: 10.3969/j.issn.1672-7703.2020.02.001
    [6] 赵文智, 胡素云, 侯连华等. 中国陆相页岩油类型、资源潜力及与致密油的边界. 石油勘探与开发, 2020, 47(1): 1-10 (Zhao Wenzhi, Hu Suyun, Hou Lianhua, et al. Types and resource potential of continental shale oil in China and its boundary with tight oil. Petroleum Exploration and Development, 2020, 47(1): 1-10 (in Chinese) doi: 10.1016/S1876-3804(20)60001-5
    [7]

    Feng QH, Xu SQ, Xing XD, et al. Advances and challenges in shale oil development: A critical review. Advances in Geo-Energy Research, 2020, 4(4): 406-418 doi: 10.46690/ager.2020.04.06

    [8] 贾承造, 邹才能, 李建忠等. 中国致密油评价标准、主要类型、基本特征及资源前景. 石油学报, 2012, 33(3): 343-350 (Jia Chengzao, Zou Caineng, Li Jianzhong, et al. Assessment criteria, main types, basic features and resource prospects of the tight oil in China. Acta Petrolei Sinica, 2012, 33(3): 343-350 (in Chinese) doi: 10.7623/syxb201203001
    [9] 邹才能, 丁云宏, 卢拥军等. “人工油气藏”理论、技术及实践. 石油勘探与开发, 2017, 44(1): 144-154 (Zou Caineng, Ding Yunhong, Lu Yongjun, et al. Concept, technology and practice of “man-made reservoirs” development. Petroleum Exploration and Development, 2017, 44(1): 144-154 (in Chinese)
    [10] 黄荣樽. 水力压裂裂缝的起裂和扩展. 石油勘探与开发, 1981(5): 62-74 (Huang Rongzun. fracture and propagation of hydraulic fracturing fractures. Petroleum Exploration and Development, 1981(5): 62-74 (in Chinese)
    [11] 邓燕. 大位移井水力压裂裂缝起裂机理研究及应用[硕士论文]. 成都: 西南石油大学, 2003

    (Deng Yan. Research and application of fracture initiation mechanism of hydraulic fracturing in extended reach wells. [Master's Thesis]. Cheng Du: Southwest Petroleum University, 2003 (in Chinese))

    [12]

    Zeng FH, Cheng XZ, Guo JC, et al. Investigation of the initiation pressure and fracture geometry of fractured deviated wells. Journal of Petroleum Science and Engineering, 2018, 165: 412-427 doi: 10.1016/j.petrol.2018.02.029

    [13] 刘京. 基于扩展有限元的页岩水平井压裂裂缝扩展规律研究[硕士论文]. 西安: 西安石油大学, 2019

    (Liu Jing. Fracture Propagation Law of Shale Horizontal Wells Based on Expansion Finite Element Method [Master’s Thesis]. Xi an: Xi’an Shiyou university, 2019 (in Chinese))

    [14]

    Liu ZY, Jin Y, Chen M, et al. Analysis of non-planar multi-fracture propagation from layered-formation inclined-well hydraulic fracturing. Rock Mechanics and Rock Engineering, 2016, 49(5): 1747-1758 doi: 10.1007/s00603-015-0872-1

    [15] 唐慧莹, 张烈辉, 张东旭. 页岩气藏水平井分段压裂缝间应力干扰全三维模拟研究//2018IFEDC油气田勘探与开发国际会议. 2018

    (Tang Huiying, Zhang Liehui, Zhang Dongxu et al. The investigation of stress interference among hydraulic fractures of multi-stage horizontal wells in shale gas reservoirs with fully 3D simulation//2018 the IFEDC Committee. 2018 (in Chinese))

    [16] 张广明, 刘勇, 刘建东等. 页岩储层体积压裂的地应力变化研究. 力学学报, 2015, 47(6): 965-972 (Zhang Guangming, Liu Yong, Liu Jiandong, et al. Research on the geostress change of shale reservoir volume fracturing. Chinese Journal of Theoretical and Applied Mechanics, 2015, 47(6): 965-972 (in Chinese) doi: 10.6052/0459-1879-15-274
    [17]

    Jia P, Cheng LS, Huang SJ, et al. A comprehensive model combining Laplace-transform finite-difference and boundary-element method for the flow behavior of a two-zone system with discrete fracture network. Journal of Hydrology, 2017, 551: 453-469 doi: 10.1016/j.jhydrol.2017.06.022

    [18] 蔡建超, 郭士礼, 游利军等. 裂缝-孔隙型双重介质油藏渗吸机理的分形分析. 物理学报, 2013, 62(1): 228-232 (Cai Jianchao, Guo Shili, You Lijun, et al. Fractal analysis of spontaneous imbibition mechanism in fractured-porous dual media reservoir. Acta Physica Sinica, 2013, 62(1): 228-232 (in Chinese)
    [19]

    Meng Q, Cai J. Recent advances in spontaneous imbibition with different boundary conditions. Capillarity, 2018, 1(3): 19-26 doi: 10.26804/capi.2018.03.01

    [20] 樊冬艳, 姚军, 王子胜等. 基于不同倾角的压裂水平井试井解释. 水动力学研究与进展A辑, 2009, 24(6): 705-712 (Fan Dongyan, Yao Jun, Wang Zisheng, et al. Well testing on fractured horizontal well with different dip angles. Chinese Journal of Hydrodynamics, 2009, 24(6): 705-712 (in Chinese)
    [21]

    AL Rbeawi SJH, Djebbar T. Transient pressure analysis of a horizontal well with multiple inclined hydraulic fractures using type-curve matching//SPE International Symposium and Exhibition on Formation Damage Control. Lafayette, Louisiana, USA: Society of Petroleum Engineers, 2012

    [22]

    AL Rbeawi SJH, Tiab D. Partially Penetrating Hydraulic Fractures: Pressure Responses and Flow Dynamics//All Days. Oklahoma City, Oklahoma, USA: SPE, 2013: SPE-164500-MS

    [23] 贾品, 程林松, 黄世军等. 有限导流压裂定向井不稳定压力分析模型. 石油学报, 2015, 36(4): 496-503 (Jia Pin, Cheng Linsong, Huang shijun, et al. Transient pressure analysis model for finite conductivity fracture directional well. Acta Petrolei Sinica, 2015, 36(4): 496-503 (in Chinese) doi: 10.7623/syxb201504011
    [24] 贾品, 程林松, 黄世军等. 有限导流压裂定向井耦合流动模型. 计算物理, 2014, 31(5): 559-566 (Jia Pin, Cheng Linsong, Huang shijun, et al. A coupling flow model of finite-conductivity fractured directional well. Chinese Journal of Computational Physics, 2014, 31(5): 559-566 (in Chinese) doi: 10.3969/j.issn.1001-246X.2014.05.007
    [25]

    Jia P, Cheng LS, Huang SJ, et al. Transient behavior of complex fracture network. Journal of Petroleum Science and Engineering, 2015, 132: 1-17 doi: 10.1016/j.petrol.2015.04.041

    [26]

    Teng BL, Li HZ. A semi-analytical model for characterizing the pressure transient behavior of finite-conductivity horizontal fractures. Transport in Porous Media, 2018, 123(2): 367-402 doi: 10.1007/s11242-018-1047-9

    [27] 万义钊, 刘曰武, 吴能友等. 基于离散裂缝的多段压裂水平井数值试井模型及应用. 力学学报, 2018, 50(1): 147-156 (Wan Yizhao, Liu Yuewu, Wu Nengyou, et al. A numerical well test model for multi-fractured horizontal wells based on discrete-fracture model and its application. Chinese Journal of Theoretical and Applied Mechanics, 2018, 50(1): 147-156 (in Chinese) doi: 10.6052/0459-1879-17-319
    [28]

    Teng BL, Li HZ. Pressure-transient behavior of partially penetrating inclined fractures with a finite conductivity. SPE Journal, 2019, 24(02): 811-833 doi: 10.2118/194189-PA

    [29]

    Teng BL, Li HZ. A semianalytical model for evaluating the performance of a refractured vertical well with an orthogonal refracture. SPE Journal, 2019, 24(2): 891-911 doi: 10.2118/194216-PA

    [30] 王苏冉. 裂缝性弱挥发碳酸盐岩油藏渗流理论研究及应用[博士论文]. 北京: 中国石油大学(北京). 2020

    (Wang Suran. Study on seepage flow theory and application of fractured weakly volatile carbonate oil reservoirs. [PhD Thesis]. Beijing: China University of Petroleum (Beijing). 2020 (in Chinese))

    [31] 朱光普, 姚军, 樊冬艳等. 页岩气藏压裂水平井试井分析. 力学学报, 2015, 47(6): 945-954 (Zhu Guangpu, Yao Jun, Fan Dongyan, et al. Pressure transient analysis of fractured horizontal well in shale gas reservoir. Chinese Journal of Theoretical and Applied Mechanics, 2015, 47(6): 945-954 (in Chinese) doi: 10.6052/0459-1879-15-229
    [32]

    Wei C, Cheng SQ, Tu K, et al. A hybrid analytic solution for a well with a finite-conductivity vertical fracture. Journal of Petroleum Science and Engineering, 2020, 188: 106900 doi: 10.1016/j.petrol.2019.106900

    [33]

    Wu MT, Luo WJ, Wang XD. A new algorithm for computation of horizontal-well pressure in Laplace domain. Advances in Geo-energy Research, 2018, 2(4): 393-403 doi: 10.26804/ager.2018.04.04

    [34]

    Xu GH, Yin HJ, Yuan HF, et al. Decline curve analysis for multiple-fractured horizontal wells in tight oil reservoirs. Advances in Geo-Energy Research, 2020, 4(3): 296-304 doi: 10.46690/ager.2020.03.07

    [35] 赵欢, 李玮, 强小军等. 致密砂岩储层多裂缝扩展形态及影响因素. 东北石油大学学报, 2020, 44(5): 76-81+122+9 (Zhao Huan, Li Wei, Qiang Xiaojun, et al. Fracture morphology and the influence factors of multi-fracture propagation in tight sand reservoir. Journal of Northeast Petroleum University, 2020, 44(5): 76-81+122+9 (in Chinese)
    [36] 程林松. 高等渗流力学. 北京: 石油工业出版社, 2011: 1-386

    (Cheng Linsong. Advanced Mechanics of Fluids in Porous Media. Beijing: Petroleum Industry Press, 2011: 1-386 (in Chinese))

    [37]

    Peaceman DW. Interpretation of well-block pressures in numerical reservoir simulation with nonsquare grid blocks and anisotropic permeability. Society of Petroleum Engineers Journal, 1983, 23(3): 531-543 doi: 10.2118/10528-PA

    [38]

    Van Everdingen AF, Hurst W. The application of the laplace transformation to flow problems in reservoirs. Journal of Petroleum Technology, 1949, 1(12): 305-324 doi: 10.2118/949305-G

图(15)
计量
  • 文章访问数:  1349
  • HTML全文浏览量:  453
  • PDF下载量:  86
  • 被引次数: 0
出版历程
  • 收稿日期:  2021-04-27
  • 录用日期:  2021-06-24
  • 网络出版日期:  2021-06-25
  • 刊出日期:  2021-08-17

目录

/

返回文章
返回