NUMERICAL SIMULATION FOR COUPLING FLOW AND GEOMECHANICS IN EMBEDDED DISCRETE FRACTURE MODEL BASED ON XFEM-MBEM
-
摘要: 离散缝网的表征与模拟是目前国内外研究的热点. 在非常规油气开发过程中, 由于地应力场的存在会对裂缝的流动属性产生显著影响, 若将裂缝视为静态对象, 与矿场数据会出现极大偏差, 因此要基于动态裂缝做更深入的研究. 本文针对致密油藏应力场−渗流场耦合力学问题, 提出了一种高效的混合数值离散化方法, 其中采用扩展有限元法 (XFEM) 求解岩石的弹性形变, 采用了混合边界元法 (MBEM) 精确计算基岩与裂缝间的非稳态窜流, 这两种数值格式是完全耦合的, 并对整体计算格式的时间项进行了全隐式求解, 可准确表征致密油藏开采过程中的裂缝变形及流体流动机理. 此外, 本文采用了嵌入式离散裂缝前处理算法显式表征大尺度水力压裂缝, 并考虑了支撑剂的作用; 采用了双孔有效应力原理和双重介质隐式裂缝表征方法, 可捕捉基质与小尺度天然裂缝的动态信息; 由此, 本文所提出的混合模型综合表征了基质−天然裂缝−水力压裂缝共同组成的致密油藏复杂渗流环境, 并通过几个实例论证了模型的准确性, 研究表明: 对致密油藏压裂水平井进行产能评价时, 应力场所引起渗流参数的改变及裂缝开度降低的影响不可忽略. 本文研究可为非常规油气资源的开发提供理论指导.Abstract: The characterization and simulation of discrete fracture networks is a hot topic at home and abroad. In the development process of unconventional oil or gas reservoir, the in-situ stress field will have a significant impact on the flow properties of fractures. If fractures are regarded as static objects, there will be a great deviation from the field data. Therefore, more in-depth research should be done based on dynamic fractures. In this paper, an efficient hybrid numerical discretization method is proposed to solve the coupled mechanical problems of coupling geomechanics and fluids flow in tight oil reservoirs. The extended finite element method (XFEM) is used to solve the elastic deformation of rock, and the mixed boundary element method (MBEM) is adopted to accurately calculate the unsteady flux between matrix and fracture. The two numerical schemes are fully-coupled and the time-terms in overall calculation scheme is solved by the fully-implicit method, which can accurately and efficiently simulate the mechanism of fracture deformation and fluids flow in the development of tight oil reservoirs. In addition, the embedded pre-treatment is used to characterize the large-scale hydraulic fracture, and the effect of proppant is considered. The dynamic information of matrix and small-scale natural fracture can be captured by using the double-porosity effective stress principle and the characterization method of implicit fracture in dual-media. Therefore, the hybrid model proposed in this paper comprehensively characterizes the complex system composed of matrix, natural fractures and hydraulic fractures. The accuracy of proposed model is demonstrated by several examples in this paper. The study shows that the influence of the flow parameters change and the fracture aperture reduction caused by stress-field can not be ignored when evaluating the productivity of fractured horizontal well in tight oil reservoirs. This work can provide theoretical guidance for the development of unconventional oil and gas resources.
-
引 言
非常规裂缝性储层中, 岩石基质提供油气储存空间, 人工裂缝、诱导裂缝和天然裂缝提供其流动通道[1]. 目前油田实际历史数据已表明, 地质力学特性对裂缝性油藏的产量起重要作用, 会导致一些关键渗流参数发生改变, 如基岩的孔隙度和渗透率; 此外, 致密油藏以“缝控储量”为主[2], 岩石形变会引起裂缝开度变化, 这对于压裂水平井产能评价的影响极大. 传统的应力敏感模型无法揭示深层次的科学道理[3-4], 为了准确模拟此类裂缝性多孔介质中的裂缝动态行为和流体流动机理, 需要综合考虑流体和岩石两方面的力学性质.
岩石力学中有效应力的概念最早由Terzaghi[5]和Biot[6]提出, 并由Biot[7]发展为较完善的三维固结理论, 这是研究应力场与渗流场耦合的早期基础. 为了将上述的理论基础应用于传统的渗流模型, Barenblatt和Zheltov[8]首次提出了双重孔隙度/双重渗透率的概念, 基岩与裂缝为两种重叠的、不同的连续介质, 其中基岩完全被裂缝所围绕, 用于表征裂缝性多孔介质中的流体和岩石形变. 后续, Bai[9]提出了关于双重介质更为严格的物理解释即双孔有效应力理论, 用于分析不同尺度的裂缝和基质的形变. 针对多孔介质流固耦合数学模型的建立, Wilson和Aifantis[10]、Beskos和Aifantis[11]和Khaled等[12]建立了基于有限元法的双重介质流固耦合模型, 该模型的主要问题在于其忽视了基质与裂缝之间压差导致的耦合流动作用; Lewis和Schrefler[13]对于油藏中的流固耦合问题进行了研究, 分别考虑了单相流、多相流、线弹性变形以及弹塑性变形, 并给出了相应的有限元计算格式. 在上述早期的研究中, 采用连续介质力学理论表征裂缝性油藏, 将不连续的裂缝系统粗化为连续介质, 能在一定程度上刻画基岩和裂缝的流体流动及力学特征, 但对于裂缝的几何形状及传导率具有较强的约束, 更适用于描述小尺度裂缝, 对于大尺度裂缝处理精度较低, 无法满足以体积压裂为开发手段的非常规油气藏的矿场实际需求.
近年来, 基于离散裂缝的流固耦合理论受到了广泛关注, 此类模型将裂缝离散化显式表征, 能准确刻画大尺度裂缝对渗流场和应力场的影响, 根据前处理算法的不同, 可划分为离散裂缝模型[14] (DFM) 与嵌入式离散裂缝模型[15-16] (EDFM), 其中嵌入式裂缝模型最早由Li和Lee[17]开发, 旨在提高天然裂缝建模的精确程度. 与使用复杂非结构网格在空间上匹配裂缝几何结构的DFM方法相比, EDFM只需使用一组固定的结构化矩形网格, 裂缝在前处理中与网格分开并显式精确刻画, 这是EDFM相对于DFM的关键优势[18]. 流固耦合研究方面, 由于裂缝开度与油藏尺度存在极大差异, 需对裂缝进行降维处理避免裂缝周围局部加密以减小网格量, 而降维后裂缝两侧的位移场存在间断性[19]. 由于EDFM采用非匹配性网格, 使得其在表征裂缝力学间断效应方面极具挑战性. 目前主要处理方法有: 有限元法、位移不连续法和扩展有限元法[20] (XFEM). 其中, XFEM与EDFM具有相似的优点, 两者均基于结构化网格, 分别相较于传统的基于有限元法的应力场、渗流场计算, 避免了网格重构和非结构化网格剖分带来的困难. 基于此, Ren等[21-22]、Yan等[23]、Zeng等[24]先后实现了将EDFM与用于应力场计算的扩展有限元法的耦合研究; Ren等[25]进一步开发了pEDFM-XFEM流固耦合数值模拟器, 能更精确地模拟多相流情况下的岩石和流体的力学性质.
在上述前人的研究中, 仍亟待解决的主要问题有两点: 致密基质的渗透率极低, 大规模体积压裂之后储层中的大尺度水力压裂缝与大量被激活的、小尺度的天然裂缝共存, 对此类“人工裂缝性油藏”, 单一数学模型无法进行准确表征, 需要更一步发展综合数学模型; 目前现有数值模拟方法在计算包含裂缝单元的基质网格内的压力分布时, 采用了线性分布假设, 这导致了早期及非稳态流动计算精度的不足. 本文以上述突出的科学问题为导向, 建立了致密油藏基质-天然裂缝-水力压裂缝综合数学模型, 基于单孔模型表征基质渗流参数随着压力、应变场的改变, 基于双重介质模型表征天然裂缝渗流参数的变化, 基于嵌入式离散裂缝模型表征水力压裂缝的动态特征; 并提出了XFEM-MBEM的混合数值离散化方法, 用于计算基岩和裂缝间的非稳态窜流, 提高了模拟的早期精度, 可准确表征致密油藏开采过程中的裂缝变形及流体流动机理, 此研究为致密油气开发领域提供了理论基础.
1. 流固耦合数学模型
1.1 物理模型及基本假设
本文采用了黑油模型假设, 考虑基质和流体的压缩性, 油水两相流动均符合达西定律, 等温渗流, 岩石变形考虑为微小线弹性变形, 力学方程中采用常规弹性力学符号约定, 即拉伸为正, 压缩为负.
1.2 准静态力学平衡方程
图1为弹性静力平衡状态下的多孔介质, 中间存在一条水力压裂缝,
$ {p}_{\mathrm{F}} $ 为作用裂缝内表面的缝内孔隙压力,$ {p}_{\mathrm{p}} $ 为支撑剂颗粒的支撑力,$ {\boldsymbol{n}}_{\mathrm{F}}^{+} $ 和$ {\boldsymbol{n}}_{\mathrm{F}}^{-} $ 表示裂缝表面两侧$ {\varGamma }_{\mathrm{F}} $ 上的两个单位法向量,$ {\boldsymbol{n}}_{\mathrm{t}} $ 表示基岩应力外边界$ {{\varGamma }}_{\mathrm{t}} $ 上的单位外法向量,$ \boldsymbol{t} $ 为应力外边界上的牵引力,$ {\bar{\boldsymbol{u}}} $ 为位移外边界$ {\varGamma }_{\mathrm{u}} $ 上的平均固体位移. 由此可建立准静态应力平衡方程式(1) ~ 式(4)$$ \nabla \cdot {\boldsymbol{\sigma }}\left( {{\boldsymbol{u}},{p_{\rm{m}}}} \right) + {\rho _{\rm{b}}}{\boldsymbol{g}} = 0,\;{\text{ on }}{\varOmega} $$ (1) $$ {\boldsymbol{\sigma }} \cdot {{\boldsymbol{n}}_{\rm{t}}} = {\boldsymbol{t}},\;{\text{ on }}{{\varGamma} _{\rm{t}}} $$ (2) $$ {\boldsymbol{u}} = \overline {\boldsymbol{u}},\; {\text{ on }}{{\varGamma} _{\rm{u}}} $$ (3) $$ {\boldsymbol{\sigma }} \cdot {\boldsymbol{n}}_{\rm{F}}^ + = - {\boldsymbol{\sigma }} \cdot {\boldsymbol{n}}_{\rm{F}}^ - {\text{ = }}\left( {{p_{\rm{F}}} + {p_{\rm{p}}}} \right){\boldsymbol{n}}_{\rm{F}}^ - ,\;{\text{ on }}{{\varGamma} _{\rm{F}}} $$ (4) 式中,
$ \boldsymbol{\sigma } $ 为总应力张量, Pa;$ \boldsymbol{u} $ 为岩石固体位移, m;$ {P}_{\mathrm{m}} $ 为基岩孔隙压力, Pa;$ {\mathit{\rho }}_{\mathrm{b}} $ 为多孔介质基岩密度, kg/m3;$ \boldsymbol{g} $ 为重力加速度, N/kg.1.3 油水渗流控制方程
基质体积单元中的油水两相质量平衡方程
$$ {\left[ {\frac{\partial }{{\partial t}}\frac{{({\rho _{\rm{o}}}\phi _{\rm{m}}^ * {s_{\rm{o}}})}}{{{B_{\rm{o}}}}} - \nabla \cdot \left( {{\rho _{\rm{o}}}\frac{{k{k_{\rm{ro}}}}}{{{\mu _{\rm{o}}}{B_{\rm{o}}}}}\nabla {p_{\rm{o}}}} \right)} \right]^m} = {\left[ {{\rho _{\rm{o}}}{q_{\rm{o}}}} \right]^{mF}} + {\left[ {{\rho _{\rm{o}}}{q_{\rm{o}}}} \right]^{mw}} $$ (5) $$ {\left[ {\frac{\partial }{{\partial t}}\frac{{({\rho _{\rm{w}}}\phi _{\rm{m}}^ * {s_{\rm{w}}})}}{{{B_{\rm{w}}}}} - \nabla \cdot \left( {{\rho _{\rm{w}}}\frac{{k{k_{\rm{rw}}}}}{{{\mu _{\rm{w}}}{B_{\rm{w}}}}}\nabla {p_{\rm{w}}}} \right)} \right]^m} = {\left[ {{\rho _{\rm{w}}}{q_{\rm{w}}}} \right]^{mF}} + {\left[ {{\rho _{\rm{w}}}{q_{\rm{w}}}} \right]^{mw}} $$ (6) 离散裂缝单元中的油水两相质量平衡方程
$$ {\left[ {\frac{\partial }{{\partial t}}\frac{{({\rho _{\rm{o}}}\phi _{\rm{F}}^ * {s_{\rm{o}}})}}{{{B_{\rm{o}}}}} - \nabla \cdot \left( {{\rho _{\rm{o}}}\frac{{k{k_{\rm{ro}}}}}{{{\mu _{\rm{o}}}{B_{\rm{o}}}}}\nabla {p_{\rm{o}}}} \right)} \right]^F} = {\left[ {{\rho _{\rm{o}}}{q_{\rm{o}}}} \right]^{mF}} + {\left[ {{\rho _{\rm{o}}}{q_{\rm{o}}}} \right]^{Fw}} $$ (7) $$ {\left[ {\frac{\partial }{{\partial t}}\frac{{({\rho _{\rm{w}}}\phi _{\rm{F}}^ * {s_{\rm{w}}})}}{{{B_{\rm{w}}}}} - \nabla \cdot \left( {{\rho _{\rm{w}}}\frac{{k{k_{\rm{rw}}}}}{{{\mu _{\rm{w}}}{B_{\rm{w}}}}}\nabla {p_{\rm{w}}}} \right)} \right]^F} = {\left[ {{\rho _{\rm{w}}}{q_{\rm{w}}}} \right]^{mF}} + {\left[ {{\rho _{\rm{w}}}{q_{\rm{w}}}} \right]^{Fw}} $$ (8) 附加方程
$$ \qquad\qquad\qquad {s_{\rm{o}}} + {s_{\rm{w}}} = 1 $$ (9) $$ \qquad\qquad\qquad {p_{\rm{w}}} = {p_{\rm{o}}} - {p_{\rm{c}}} $$ (10) 式中,
$ k $ 为渗透率, m2;$ {k}_{\mathrm{r}\mathrm{o}} $ 为油相相对渗透率, 小数;$ {k}_{\mathrm{r}\mathrm{w}} $ 为水相相对渗透率, 小数;$ \mu $ 为流体黏度, mPa·s;$ s $ 为流体饱和度, 小数;$ B $ 为流体体积系数, 小数;$ p $ 为流体压力, Pa; 要补充的是, 上述式中下标o和w分别表示油相和水相, 上标m和F分别表示基岩和水力压裂缝, 上标w表示井筒.$ {\left[\rho q\right]}^{mF} $ 表示基岩与裂缝之间的窜流量;$ {\left[\rho q\right]}^{mw} $ 表示基岩与井筒之间的流量;$ {\phi }^{*} $ 表示为Lagrangen孔隙度, 后续进行详细阐述.1.4 基质-裂缝传质方程
EDFM中存在三类非相邻连接关系 (NNRs): 位于不同基质网格中的相邻裂缝单元之间的连接; 位于同一基质网格中的相邻裂缝单元之间的连接; 基质网格与所包含裂缝单元之间的连接. 其中, 传导率可表示为流度与其几何因子的乘积, 以油相为例, 流度和几何因子可表示为
$$ {\lambda _{ij}} = \dfrac{{{k_{{\rm{ro}}({\rm{upstream}})}}}}{{\left( {\dfrac{{{\mu _{{\rm{o}},i}} + {\mu _{{\rm{o}},j}}}}{2}} \right)\left( {\dfrac{{{B_{{\rm{o}},i}} + {B_{{\rm{o}},j}}}}{2}} \right)}} $$ (11) $$ {G_{ij}} = \frac{{{G_i}{G_j}}}{{{G_i} + {G_j}}} $$ (12) $$ T{I_{ij}} = {\lambda _{ij}} \cdot {G_{ij}} $$ (13) 式中,
$ {\lambda }_{ij} $ 表示单元$ i $ 和$ j $ 之间的油相流度;$ {k}_{\mathrm{r}\mathrm{o}\left(\mathrm{u}\mathrm{p}\mathrm{s}\mathrm{t}\mathrm{r}\mathrm{e}\mathrm{a}\mathrm{m}\right)} $ 表示利用上游权格式取值的油相相对渗透率;$ {G}_{ij} $ 为单元$ i $ 和$ j $ 之间的几何因子;$ {TI}_{ij} $ 表示单元$ i $ 和$ j $ 之间的传导系数.基质与裂缝之间传质项的处理是基于EDFM中的第三类连接关系, 原始EDFM方法[26]在计算包含裂缝单元的基质网格内的压力分布时, 采用了Lee等[27]及Praditia等[28]提出的压力与到裂缝面的垂向距离成正比的线性分布假设如下
$$ \qquad\qquad\qquad {G_i} = \frac{{{k_i}A}}{{\left\langle d \right\rangle }} \qquad$$ (14) $$ \qquad\qquad\qquad \left\langle d \right\rangle = \dfrac{{\displaystyle\int {\left| {{\boldsymbol{n}} \cdot x} \right|{\rm{d}}S} }}{S} $$ (15) 式中
$ {k}_{i} $ 表示单元$ i $ 的渗透率, m2;$ A $ 代表接触面积, m2;$ \left\langle{d}\right\rangle $ 代表基质网格中心点与裂缝面之间的距离, m;$ S $ 为所在基质网格的面积, m2.但在实际物理情景中由于裂缝与基质的渗透率级别相差极大, 这种求取“平均法向距离”的近似方法在处理非稳态传质时会导致一定的误差. 因此, 本文基于Fang等[29]与Cao等[30]所提出的混合边界元法, 利用稳态渗流控制方程的边界积分方程推导得到双孔基质网格向裂缝网格传质量的新近似格式, 以油相为例见式(16), EDFM中MBEM流量处理的连接关系示例见图2所示.
![]()
图 2 双孔系统EDFM中MBEM流量项处理的连接关系示例Figure 2. Example of connection list for MBEM handling flux in EDFM with dual-porosity system$$ \begin{split} & {\left[ {{q_{{\rm{omF}}}}} \right]_{({n_{\rm{F}}} \times 1)}} = \\ &\quad \frac{{k{k_{{\rm{ro}}}}h}}{{{\mu _{\rm{o}}}{B_{\rm{o}}}}}{\left( {{{{{\boldsymbol{G}}_{\rm{F}}}}_{({n_{\rm{F}}} \times 4)}}{\boldsymbol{G}}_{(4 \times 4)}^{ - 1}{{ {{{\boldsymbol{G}}_{{\rm{mF}}}}} }_{\left. {(4 \times {n_{\rm{F}}}} \right)}} - {{ {{{\boldsymbol{G}}_{\rm{FF}}}} }_{({n_{\rm{F}}} \times {n_{\rm{F}}})}}} \right)^{ - 1}} \cdot \\ &\quad \Bigg\{ {\left[ {{{\left( {\frac{{\partial {{\boldsymbol{G}}_{\rm{F}}}}}{{\partial {\boldsymbol{n}}}}} \right)}_{({n_{\rm{F}}} \times 4)}} - {{{{\boldsymbol{G}}_{\rm{F}}}}_{({n_{\rm{F}}} \times 4)}}{\boldsymbol{G}}_{(4 \times 4)}^{ - 1}{{\left( {\frac{{\partial {\boldsymbol{G}}}}{{\partial {\boldsymbol{n}}}}} \right)}_{(4 \times 4)}}} \right]} \cdot \Bigg.\\ &\quad \Bigg. { {{{\boldsymbol{G}}_{{\rm{conv}}}}} _{(4 \times 5)}}{ {{{\boldsymbol{p}}_{\rm{o}}}} _{(5 \times 1)}} - { {{{\boldsymbol{p}}_{\rm{F}}}} _{({n_{\rm{F}}} \times 1)}}\Bigg\} \end{split} $$ (16) 式中,
$ {q}_{\mathrm{o}\mathrm{m}\mathrm{F}} $ 为从裂缝单元流入基质网格的油相流量, m3; 下标$ {n}_{\mathrm{F}} $ 代表裂缝单元个数;$ h $ 为基质网格厚度, m;$ \boldsymbol{G} $ ,$ {\boldsymbol{G}}_{\mathrm{m}\mathrm{F}} $ ,$ {\boldsymbol{G}}_{\mathrm{F}} $ ,$ {\boldsymbol{G}}_{\mathrm{F}\mathrm{F}} $ ,$ \dfrac{\partial \boldsymbol{G}}{\partial {\boldsymbol{n}}} $ 等矩阵形式见附录.1.5 井方程
经典Peaceman公式[31]可用于计算井筒与基质或裂缝之间的流量. 以油相为例, 对于压裂水平井井筒与裂缝单元相连接
$$ \left[ {{\rho _{\rm{o}}}q_{\rm{o}}} \right]^{Fw} = {\rho _{\rm{o}}}\dfrac{{{{\left( {k{k_{{\rm{ro}}}}} \right)}^F}}}{{{\mu _{\rm{o}}}{B_{\rm{o}}}}} \cdot \dfrac{{ {\dfrac{{2{\text{π}} {w_{\rm{F}}}}}{{\cos \theta }}} }}{{\ln \left( {\dfrac{{{r_{\rm{e}}}}}{{{r_{\rm{w}}}}}} \right)}}{\left( {{p^F} - {p^w}} \right)_{\rm{o}}} $$ (17) $$ {r_{\rm{e}}} = 0.14{\left( {{l_1}^2 + {l_2}^2} \right)^{1/2}} $$ (18) 式中,
$ {w}_{\mathrm{F}} $ 为水力压裂缝的开度, m;$ {r}_{\mathrm{e}} $ 为等效补给井径, m;$ {r}_{\mathrm{w}} $ 为油井井径, m;$ {l}_{1} $ 和$ {l}_{2} $ 分别为裂缝矩形单元的两个边长, m.2. 多孔介质本构关系
2.1 修正的Biot方程
为准确描述非常规储层中基质与小尺度裂缝间的流动和岩石形变, 可采用由Bai[32]提出的双孔有效应力原理, 对基质和裂缝分别建立具有严格物理意义的有效应力关系
$$ {{\boldsymbol{\sigma }}_{\rm{m}}} = {{\boldsymbol{\sigma '}}_{\rm{m}}} - {\alpha _{\rm{m}}}{p_{\rm{m}}}{\boldsymbol{I}} $$ (19) $$ {{\boldsymbol{\sigma }}_{\rm{f}}} = {{\boldsymbol{\sigma '}}_{\rm{f}}} - {p_{\rm{f}}}{\boldsymbol{I}} $$ (20) $$ {\boldsymbol{\sigma }} = {{\boldsymbol{\sigma }}_{\rm{m}}} = {{\boldsymbol{\sigma }}_{\rm{f}}} $$ (21) $$ {\boldsymbol{\varepsilon }} = {{\boldsymbol{\varepsilon }}_{\rm{m}}} + {{\boldsymbol{\varepsilon }}_{\rm{f}}} $$ (22) 式中,
$ {\boldsymbol{\sigma }}_{\mathrm{m}} $ 表示基质总应力张量, Pa;$ {\boldsymbol{\sigma }}_{\mathrm{f}} $ 为裂缝总应力张量, Pa;$ {\boldsymbol{\sigma }}_{\mathrm{m}}' $ 为基质有效应力, Pa;$ {\boldsymbol{\sigma }}_{\mathrm{f}}' $ 表示裂缝有效应力, Pa;$ {\alpha }_{\mathrm{m}} $ 表示毕渥系数, 小数;$ \boldsymbol{\varepsilon } $ 为总应变张量, m;$ {\boldsymbol{\varepsilon }}_{\mathrm{m}} $ 为基质应变张量, m;$ {\boldsymbol{\varepsilon }}_{\mathrm{f}} $ 为裂缝应变张量, m; 上述式中下标f表示小尺度裂缝.Ren根据式(19) ~ 式(22)及胡克定律, 指出双重介质有效应力的表达式为[22]
$$ \begin{split} & {\boldsymbol{\sigma '}} = {{\boldsymbol{D}}_{{\rm{mf}}}}:{\boldsymbol{\varepsilon }} - {{\boldsymbol{D}}_{{\rm{mf}}}}:{{\boldsymbol{C}}_{\rm{m}}}:{\alpha _{\rm{m}}}{p_{\rm{m}}}{\boldsymbol{I}} \\ &\qquad - {{\boldsymbol{D}}_{{\rm{mf}}}}:{{\boldsymbol{C}}_{\rm{f}}}:{p_{\rm{f}}}{\boldsymbol{I}} \end{split} $$ (23) $$ {\boldsymbol{\varepsilon}} = \frac{1}{2}\left( {{\nabla ^{\text{T}}}{\boldsymbol{u}} + \nabla {\boldsymbol{u}}} \right) $$ (24) 式中,
$ {{\boldsymbol{D}}}_{\mathrm{m}\mathrm{f}} $ 为弹性刚度张量, Pa;$ {{\boldsymbol{C}}}_{\mathrm{m}} $ 为基质柔度张量, Pa;$ {{\boldsymbol{C}}}_{\mathrm{f}} $ 为裂缝柔度张量, Pa.2.2 孔渗参数表征方法
孔隙度是关于体积应变和流体压力的函数[25]. 考虑线弹性微小形变下的基岩Lagrangen孔隙度和大尺度裂缝Lagrangen孔隙度可分别由下式表征
$$ {\phi }_{{\rm{m}}}^{*}={\phi }_{{\rm{m}}}^{0}+{\alpha }_{{\rm{m}}}{\epsilon}_{V}+\frac{1}{M}\left({p}_{{\rm{m}}}-{p}_{{\rm{m}}}^{0}\right) $$ (25) $$ {\phi }_{{\rm{F}}}^{*}={\phi }_{{\rm{F}}}^{0}\left(1+{\epsilon}_{{\rm{V}}}\right) $$ (26) 式中
$ {\phi }_{\mathrm{m}}^{0} $ 为基岩原始孔隙度, 小数;$ {\phi }_{\mathrm{F}}^{0} $ 为大尺度裂缝原始孔隙度, 小数;$ {\phi }_{\mathrm{m}}^{0} $ 为基岩原始孔隙度, 小数;$ {p}_{\mathrm{m}} $ 为基岩流体压力, Pa;$ {p}_{\mathrm{m}}^{0} $ 为基岩原始流体压力, Pa;$ 1/M $ 物理意义同$ {\alpha }_{\mathrm{m}};{\epsilon}_{\mathrm{V}} $ 为体积应变, 即总应变张量的迹.本文中将小尺度天然裂缝和基岩视为双重连续介质, 因此小尺度裂缝表征方法与基质相同. 对于大尺度水力压裂缝, 生产过程中由于缝内流体压力损失会导致裂缝面的闭合, 然而当存在支撑剂时, 支撑剂发生形变会对裂缝面施加支撑力阻止其闭合. 根据胡克定律, 由支撑剂引起的支撑力可表示为
$$ {p_{\rm{p}}} = \left\{ {\begin{array}{*{20}{c}} {{E_{\rm{p}}}\dfrac{{ - {\boldsymbol{w}} \cdot {{\boldsymbol{n}}_{\rm{F}}}}}{{w_{\rm{F}}^0}},}&{{\boldsymbol{w}} \cdot {{\boldsymbol{n}}_{\rm{F}}} < 0} \\ {0,}&{{\text{else }}} \end{array}} \right. $$ (27) 式中,
$ {E}_{\mathrm{p}} $ 为支撑剂弹性模量, Pa;${ \boldsymbol{w} }$ 为裂缝开度的位移, m;$ {w}_{\mathrm{F}}^{0} $ 为水力压裂缝的原始开度, m.此过程中, 水力压裂缝的实际导流能力依赖于裂缝渗透率和开度的更新迭代[33]
$$ {k_{\rm{F}}} = k_{\rm{F}}^0{\left( {\frac{{{w_{\rm{F}}}}}{{w_{\rm{F}}^0}}} \right)^2} $$ (28) 式中,
$ {k}_{\mathrm{F}}^{0} $ 为水力压裂缝的原始渗透率, m2.3. 耦合模型的数值离散
3.1 有限体积法离散流动方程
XFEM与EDFM具有相似的优点, 两者均基于结构化网格, 网格划分容易且物理意义清晰. 因此本文采用正交矩形网格对基质进行几何离散, 针对EDFM前二类非相邻连接关系, 采用满足局部物质守恒且简单易行的块中心有限体积方法来获取渗流控制方程的离散格式, 仍以油相为例, 此时为公式简洁, 方程中去除可能存在的点源/汇项, 对式(5)的两侧对时间和控制体积进行积分
$$ \int_t^{t + \Delta t} {\int\limits_{\Delta V} {\nabla \cdot \left( {{\rho _{\rm{o}}}\frac{{k{k_{{\rm{ro}}}}}}{{{\mu _{\rm{o}}}{B_{\rm{o}}}}}\nabla {p_{\rm{o}}}} \right)} } {\rm{d}}V{\rm{d}}t = \int\limits_{\Delta V} {\int_t^{t + \Delta t} { {\frac{\partial }{{\partial t}}\left( {\frac{{\phi _{\rm{m}}^ * {s_{\rm{o}}}}}{{{B_{\rm{o}}}}}} \right)} } } {\rm{d}}t{\rm{d}}V $$ (29) 式中左侧利用散度定理将体积分转换为面积分, 并用符合物理意义的相邻网格之间的法向流量之和近似; 右侧对时间的积分精确求解, 对体积分采用矩形法估计, 由此式(29)可以写为
$$ \int_t^{t + \Delta t} {\sum\limits_{j = 1}^n {\left[ {{\lambda _{ij}}{G_{ij}}\left( {{p_{{\rm{o}}i}} - {p_{{\rm{o}}j}}} \right)} \right]} {\rm{d}}t} = \Delta {V_i}\left[ {{{\left( {\frac{{\phi _{\rm{m}}^ * {s_{\rm{o}}}}}{{{B_{\rm{o}}}}}} \right)}^{t + \Delta t}} - {{\left( {\frac{{\phi _{\rm{m}}^ * {s_{\rm{o}}}}}{{{B_{\rm{o}}}}}} \right)}^t}} \right] $$ (30) 对于相邻基质网格如图3, 式(30)取隐式格式可改写为
$$ \sum\limits_{j = 1}^n {\left[ {{\lambda _{ij}}{G_{ij}}\left( {p_{{\rm{o}}i}^{t + \Delta t} - p_{{\rm{o}}j}^{t + \Delta t}} \right)} \right]} = \frac{{\Delta {V_i}}}{{\Delta t}}\left[ {{{\left( {\frac{{\phi _{\rm{m}}^ * {s_{\rm{o}}}}}{{{B_{\rm{o}}}}}} \right)}^{t + \Delta t}} - {{\left( {\frac{{\phi _{\rm{m}}^ * {s_{\rm{o}}}}}{{{B_{\rm{o}}}}}} \right)}^t}} \right] $$ (31) 式中,
$ n $ 表示与该基质网格相邻的网格数量,$ {\Delta V}_{i} $ 为该基质网格的体积,$ \Delta t $ 为相邻时间步.根据式(29) ~ 式(31)的推导过程可以获得油水相的全隐式离散方程, 但对于裂缝与基质连接的第三类非相邻连接关系即含传质项, 需要进行单独考虑. 这里将已推导的多相流情况下EDFM第三类连接中基质网格与裂缝网格之间传质量的近似格式(16), 耦合到基质与裂缝之间多相流方程的有限体积离散格式中, 这也是混合边界元EDFM的核心思想. 仍以油相为例, 对于基质系统
$$ \begin{split} & \sum\limits_{j = 1}^n {\left[ {{\lambda _{o,ij}}{T_{ij}}\left( {p_{oi}^{t + \Delta t} - p_{oj}^{t + \Delta t}} \right)} \right]} - \sum\limits_{k = 1}^{{n_{\rm{F}}}} {q_{{\rm{omF}},k}^{t + \Delta t}} =\\ &\qquad \frac{{\Delta {V_i}}}{{\Delta t}}\left[ {{{\left( {\frac{{\phi _{\rm{m}}^ * {s_{\rm{o}}}}}{{{B_{\rm{o}}}}}} \right)}^{t + \Delta t}} - {{\left( {\frac{{\phi _{\rm{m}}^ * {s_{\rm{o}}}}}{{{B_{\rm{o}}}}}} \right)}^t}} \right] \end{split} $$ (32) 同理, 对于裂缝系统
$$ \begin{split} & \sum\limits_{j = 1}^n {\left[ {{\lambda _{o,ij}}{T_{ij}}\left( {p_{oi}^{t + \Delta t} - p_{oj}^{t + \Delta t}} \right)} \right]} + q_{{\rm{omF}},i}^{t + \Delta t} =\\ &\qquad \frac{{\Delta {V_i}}}{{\Delta t}}\left[ {{{\left( {\frac{{\phi _{\rm{m}}^ * {s_{\rm{o}}}}}{{{B_{\rm{o}}}}}} \right)}^{t + \Delta t}} - {{\left( {\frac{{\phi _{\rm{m}}^ * {s_{\rm{o}}}}}{{{B_{\rm{o}}}}}} \right)}^t}} \right] \end{split} $$ (33) 式中,
$ {q}_{\mathrm{o}\mathrm{m}\mathrm{F},k}^{t+\Delta t} $ 表示第$ i $ 个基质网格流入第$ k $ 个裂缝单元的油相流量, m3; 同理,$ {q}_{\mathrm{o}\mathrm{m}\mathrm{F},i}^{t+\Delta t} $ 表示基质网格向第$ i $ 个裂缝单元的油相流量, m3.3.2 扩展有限元离散力学方程
XFEM是在常规的有限元位移模式中根据单位分解的思想加引进一个跳跃函数和裂尖渐进位移场以反映位移不连续性的数值方法. XFEM中的不同节点和单元的类型见图4, 单元位移的近似解可以写成
$$ \begin{split} & {\boldsymbol{u}} = \sum\limits_{i \in I} {{N_i}} {{\boldsymbol{u}}_i} + \sum\limits_{i \in L} {\sum\limits_{j \in {\mathcal{F}_L}(i)} {{N_i}} } \left( {{H_j} - H_j^i} \right){{\boldsymbol{a}}_{i,j}} + \\ &\qquad \sum\limits_{i \in M} {\sum\limits_{j \in {\mathcal{F}_M}(i)} {{N_i}} } \left( {{J_j} - J_j^i} \right){{\boldsymbol{c}}_{i,j}} + \\ &\qquad \sum\limits_{i \in K} {\sum\limits_{j \in {\mathcal{F}_K}(i)} {\sum\limits_{l = 1}^4 {{N_i}} } } \left( {{F_{jl}} - F_{jl}^i} \right){\boldsymbol{b}}_{i,j}^l \end{split} $$ (34) 式中,
$ I $ 表示网格中所有节点的集合,$ L $ 表示包含所有裂缝段所在单元的节点集合,$ M $ 为裂缝与裂缝相交的节点集合,$ K $ 表示裂缝尖端所在单元的节点集合;$ {N}_{i} $ 为节点位移形函数,$ {\mathit{u}}_{i} $ 为节点位移矢量,$ {\boldsymbol{a}}_{i,j} $ ,$ {\boldsymbol{c}}_{i,j} $ ,$ {\boldsymbol{b}}_{i,j}^{l} $ 为单元增强节点的额外自由度,$ H $ 表示Heaviside函数,$ J $ 表示裂间连接函数,$ F $ 表示裂尖渐进函数, 三者的表达式如下所示$$ H\left( {{\xi _j}({\boldsymbol{x}})} \right) = \left\{ {\begin{array}{*{20}{l}} {1,}&{\forall {\xi _j} > 0} \\ { - 1,}&{\forall {\xi _j} < 0} \end{array}} \right. \qquad\qquad\qquad\qquad$$ (35) $$ J\left( {{\xi _{j1}}({\boldsymbol{x}}),{\xi _{j2}}({\boldsymbol{x}})} \right) = \left\{ {\begin{array}{*{20}{l}} {1,}&{\forall {\xi _{{j_1}}} \cdot {\xi _{{j_2}}} > 0} \\ {0,}&{\forall {\xi _{{j_1}}} \cdot {\xi _{{j_2}}} = 0} \\ { - 1,}&{\forall {\xi _{{j_1}}} \cdot {\xi _{{j_2}}} < 0} \end{array}} \right.\qquad\qquad $$ (36) $$ {F_{jl}}(r,\theta ) = \left\{ {\sqrt r \sin \frac{\theta }{2},\sqrt r \cos \frac{\theta }{2},\sqrt r \sin \frac{\theta }{2}\sin \theta ,\sqrt r \cos \frac{\theta }{2}\sin \theta } \right\} $$ (37) 将式(19)和式(20)代入式(1), 考虑边界条件式(2) ~ 式(4), 得到方程(1)的等效积分弱形式见式(38), 弱形式方程包含不同尺度裂缝以及基岩流体压力的相互作用
$$ \begin{split} & \int_\varOmega \delta {\boldsymbol{\varepsilon }}:{D_{{\rm{mf}}}}:{\boldsymbol{\varepsilon }}{\rm{d}}\varOmega = \int_\varOmega \delta {\boldsymbol{\varepsilon}} :\left( {\alpha _{\rm{m}}^ * {p_{\rm{m}}}{\boldsymbol{I}}} \right){\rm{d}}\varOmega + \\ & \int_\varOmega \delta {\boldsymbol{\varepsilon}} :\left( {{p_F}{\boldsymbol{I}}} \right){\rm{d}}\varOmega + \int_{{\varGamma _{\rm{t}}}} \delta {\boldsymbol{u}} \cdot {\boldsymbol{t}}{\rm{d}}{\varGamma _{\rm{t}}} + \\ & \int_{{\varGamma _{\rm{F}}}} {\left[\kern-0.15em\left[ {\delta {\boldsymbol{u}}} \right]\kern-0.15em\right] \cdot \left( {\bar p + {p_{\rm{p}}}} \right) \cdot {{\boldsymbol{n}}_{\rm{F}}}} {\rm{d}}{\varGamma _{\rm{F}}} \end{split} $$ (38) 式中,
$ \delta {\boldsymbol{\varepsilon }} $ 表示试应变,$ \delta \boldsymbol{u} $ 表示试位移,$ ⟦\delta {\boldsymbol{u}}⟧={\delta {\boldsymbol{u}}}^{+}-\delta {\boldsymbol{u}}^{-} $ .裂缝开度位移可以表示为
$$ {\boldsymbol{w}} = {{\boldsymbol{n}}_{\rm{F}}} \cdot \left[\kern-0.15em\left[ {\boldsymbol{u}} \right]\kern-0.15em\right]{{\boldsymbol{n}}_{\rm{F}}} $$ (39) 采用标准伽辽金有限元法对系统整体进行离散, 可获得一个线性平衡方程组如下
$$ {\boldsymbol{Ku}} = {\boldsymbol{f}} $$ (40) 式中, K表示整体刚度矩阵, f表示整体力向量.
3.3 整体计算格式及求解策略
本文分别采用有限体积法和扩展有限元法求解渗流控制方程及力学平衡方程, 主要变量包括基岩及裂缝中的流体压力、含水饱和度以及分别设置在网格中心和角点的节点位移, 对时间项的离散采用一阶向后欧拉全隐式差分格式, 对完全耦合的全局方程组采用Newton-Raphson迭代求解, 并利用自动微分算法计算迭代过程中的雅可比矩阵
$$ \left[ {\begin{array}{*{20}{l}} {{{\boldsymbol{A}}_{{\rm{pp}}}}}&{{{\boldsymbol{C}}_{{\rm{ps}}}}}&{{{\boldsymbol{C}}_{{\rm{pu}}}}} \\ {{{\boldsymbol{C}}_{{\rm{sp}}}}}&{{{\boldsymbol{A}}_{{\rm{ss}}}}}&{{{\boldsymbol{C}}_{{\rm{su}}}}} \\ {{{\boldsymbol{C}}_{{\rm{up}}}}}&{{{\boldsymbol{C}}_{{\rm{us}}}}}&{{{\boldsymbol{A}}_{{\rm{uu}}}}} \end{array}} \right]\left( {\begin{array}{*{20}{c}} {\delta {\boldsymbol{p}}} \\ {\delta {\boldsymbol{s}}} \\ {\delta {\boldsymbol{u}}} \end{array}} \right) = - \left( {\begin{array}{*{20}{c}} {{{\boldsymbol{R}}_{\rm{p}}}} \\ {{{\boldsymbol{R}}_{\rm{s}}}} \\ {{{\boldsymbol{R}}_{\rm{u}}}} \end{array}} \right) $$ (41) 式中系数矩阵是非对称的, 其中反对角项是主变量之间的耦合项.
4. 模型验证分析
4.1 连续性介质流固耦合算例
此算例选取了经典的一维弹性土固结的孔弹性问题, 将本文模型与Terzaghi[5]所提出的解析解进行对比. 假设等温多孔介质由单相流体和固体组成, 表现为线性孔弹性特征, 模型的宽度为5 m, 深度为10 m, 顶部为排水边界, 其余边界均不排水, 底部为零位移边界, 左右边界位移在水平方向上约束, 于顶部施加20 MPa的恒定载荷. 岩土的孔隙度为0.4, 渗透率为100 mD, 弹性模量为20 GPa, 泊松比为0.2, 孔隙比为0.4, 一维模型数值解与解析解对比结果如图5所示, 从图中可看出不同时间下本文提出的流固耦合数值解与Terzaghi解析解的结果较为吻合, 误差在合理范围之内, 能证明该模型的准确性.
![]()
图 5 一维模型数值解与解析解的对比Figure 5. Comparison between numerical solution and analytical solution of one dimensional model4.2 离散裂缝单相流算例
此算例是通过多物理场耦合模拟的商业模拟器COMSOL的标准有限元 (SFEM) 验证本文建立的流动模块. 物理模型为矩形油藏, 油藏中心有一条水力压裂缝, 设置裂缝半长为基础长度将整个区域无因次化, 裂缝无因次长度为4, 油藏在X方向的长度为10且Y方向长度为10, 见图6. 基于拉氏空间下无因次渗流方程, 将本文修正后的EDFM和有限元法与Blasingame精确解[34]进行对比, 其中裂缝无因次导流能力为200, 有限元网格裂缝采用局部加密如图7所示, 裂缝加密的网格个数为20. 本文模型采用了非匹配性的矩形网格和混合边界元法, 不对裂缝进行局部加密, 分析离散裂缝压力动态曲线中不同数值算法的适应性.
![]()
图 8 不同数值解与解析解的对比Figure 8. Comparison of different numerical solutions and analytical solution模拟结果如图8所示, 模拟无因次时间从10−9 ~ 10, 覆盖了有限导流能力的几个阶段, 可以看出修正后的EDFM方法与精确解无论从早期还是晚期都较为吻合, 而有限元模型在早期与解析解的差距很大, 中后期曲线重合, 这在图9的误差分析中也可以看出. 原因在于有限元法处理网格累积项时, 采用线性插值积分, 在一个时间步内是稳态过程, 而本文模型采用混合边界元法处理时采用了非稳态渗流的基本解, 能精确描述初期的流动特征, 根据对比结果可以说明本文模型的流动模块具有较高的早期精度.
4.3 裂缝性介质应力强度因子算例
此算例是通过经典的Tada解析解[35]验证本文建立的力学模块计算裂缝尖端应力强度因子 (SIFs) 的准确性. 考虑一个平面应变模型 (10 m×10 m), 上边界为施加均匀的拉应力, 下边界为固定位移的滚轮边界, 模型内部中心存在一条倾斜裂缝, 裂缝半长l为1 m, 如图10所示. 解析解中应力强度因子是外部应力和裂缝几何尺寸和倾角的函数, 因此为方便对比, 定义无因次应力强度因子:
$ {K}_{\mathrm{I}}^{*}={K}_{\mathrm{I}}∕\left(\sigma \sqrt{{\text{π}} l}\right) $ 及$ {K}_{\mathrm{I}\mathrm{I}}^{*}={K}_{\mathrm{I}\mathrm{I}}∕\left(\sigma \sqrt{{\text{π}} l}\right) $ , 在裂缝尺寸保持不变的条件下, 计算不同裂缝倾角无因次应力强度因子数值解与解析解的结果见图11, 对比显示本文模型的数值解与经典解析解的结果几乎保持一致, 具有较高的可信度.![]()
图 11 不同裂缝倾角无因次应力强度因子数值解与解析解的对比Figure 11. Comparison of numerical and analytical solutions for dimensionless SIFs with different fracture dip angles5. 典型实例分析
如图12所示的致密油藏, 考虑分为压裂改造区(SRV) 与非压裂改造区 (USRV), SRV内部存在一口水平井和7条水力压裂缝, 水力压裂缝中考虑支撑剂的作用. 采用平面应变模型, 油藏边界条件如图12所示. SRV内考虑微裂缝的影响, 采用双孔模型; USRV不考虑微裂缝的影响, 采用单孔模型. 致密油藏的基本参数见表1, 设置该水平井为定压生产1000 d, 井底流压为10 MPa, 并将模拟的结果与商业模拟器COMSOL的标准有限元 (SFEM) 作对比, 图13为两种数值模拟器的网格剖分, 其中COMSOL-SFEM中对SRV内部及裂缝尖端进行了局部加密处理, 本文模型XFEM-MBEM中采用非匹配性的矩形网格. 图14为两种数值模拟器在开发1000 d时刻的压力场对比, 图15和图16分别为X-位移场与Y-位移场的对比, 从上述场图可看出两种数值模拟器的结果基本一致. 以图12中SRV内部最中间的那条水力压裂缝为研究目标, 分析其不同时刻的开度分布见图17, 图中能看出XFEM-MBEM与COMSOL-SFEM不同时刻裂缝开度的计算结果较为吻合, 可以说明XFEM-MBEM流固耦合模拟的准确性. 此外, 从图中可以随着油藏的不断开发, 裂缝呈不断闭合的趋势, 定压生产下靠近井筒附近的裂缝段闭合更加剧烈. 图18为该模型在表1参数条件下考虑流固耦合作用与不考虑流固耦合作用的生产动态曲线对比, 可以看出应力场所引起的渗流参数改变及裂缝开度降低对产能的影响很大, 因此对致密油藏压裂水平井进行产能评价时, 地质力学效应的影响不可忽略.
表 1 致密油藏基本参数Table 1. Basic parameters of tight oil reservoirProperties Value Properties Value reservoir scale/m 750×340 SRV scale/m 590×220 initial reservoir pressure/MPa 25 reservoir thickness/m 10 matrix porosity 0.12 matrix permeability/mD 0.1 matrix Young’s modulus/GPa 50 matrix Poisson’s ratio 0.3 matrix Biot’s coefficient 0.8 initial water saturation 0.3 fracture porosity 0.4 fracture permeability/mD 5000 initial fracture aperture/m 0.001 fracture proppant Young's modulus/MPa 20 ![]()
图 13 两种数值模拟器的网格剖分对比Figure 13. Comparison of mesh generation between two numerical simulators![]()
图 14 两种数值模拟器开发1000 d的压力场对比(单位: MPa)Figure 14. Comparison of pressure distribution maps over 1000 days of two simulators (unit: MPa)![]()
图 15 两种数值模拟器开发1000 d的X-位移场对比(单位: mm)Figure 15. Comparison of X-displacement distribution maps over 1000 days of two simulators (unit: mm)![]()
图 16 两种数值模拟器开发1000 d的Y-位移场对比(单位: mm)Figure 16. Comparison of Y-displacement distribution maps over 1000 days of two simulators (unit: mm)![]()
图 18 考虑与不考虑流固耦合作用的生产动态曲线对比Figure 18. Comparison of production dynamic curveswith and without flow-geomechanics coupled effect6. 结 论
(1)本文提出了基于XFEM-MBEM的嵌入式离散裂缝模型流固耦合数值模拟方法, 采用正交的非匹配型网格进行几何离散, 其优势在于裂缝独立于计算网格, 与常规离散裂缝模型基于裂缝形状的匹配型网格相比, 网格划分难度大大降低, 能避免当裂缝间距或夹角非常小时, 网格划分质量不佳等问题. 本文通过几个实例分别验证了其力学模块、流动模块、耦合模块的准确性.
(2)传统EDFM方法对于第三类非相邻连接关系在计算包含裂缝单元的基质网格内的压力分布时采用了线性分布假设, 导致了开发早期及非稳态流动计算精度的不足. 局部网格加密技术可以减小包含裂缝网格的基质网格尺寸, 从而有效解决这个问题, 但会对前处理算法带来巨大的复杂性. 本文采用MBEM方法获取了双孔基质网格向裂缝网格传质量的新近似格式, 能在不加密情况下提高早期计算精度.
(3) SRV改造区内的微裂缝对提高储层流体的动用程度有重要影响. 然而, 用EDFM和XFEM进行显式表征是较为困难的, 且会增加计算耗时. 本文采用双孔有效应力原理耦合双重介质模型隐式表征SRV区内部的天然裂缝, 能有效模拟具有复合分区特征的此类裂缝性油藏问题.
(4)致密油藏开发过程中随着缝内流体被采出, 裂缝会出现不同程度的闭合, 从而导致了整体导流能力降低, 油藏采收率将降低. 考虑到非常规储层原油的产量对裂缝的依赖性如此之高. 因此, 考虑开发过程中流固耦合作用引起的参数变化对产量的评价具有重要意义.
附录
$$ {\boldsymbol{G}} = \left[ {\begin{array}{*{20}{l}} {{c_{{\rm{M}}1,1}}}&{{c_{{\rm{M}}1,2}}}&{{c_{{\rm{M}}1,3}}}&{{c_{{\rm{M}}1,4}}} \\ {{c_{{\rm{M}}2,1}}}&{{c_{{\rm{M}}2,2}}}&{{c_{{\rm{M}}2,3}}}&{{c_{{\rm{M}}2,4}}} \\ {{c_{{\rm{M}}3,1}}}&{{c_{{\rm{M}}3,2}}}&{{c_{{\rm{M}}3,3}}}&{{c_{{\rm{M}}3,4}}} \\ {{c_{{\rm{M}}4,1}}}&{{c_{{\rm{M}}4,2}}}&{{c_{{\rm{M}}4,3}}}&{{c_{{\rm{M}}4,4}}} \end{array}} \right] \tag{A1} $$ $$ {c_{{\rm{M}}i,j}} = \int_{{S_{{\rm{b}}j}}} {{G_{\rm{H}}}} \left( {{x_{{\rm{m}}i}},{y_{{\rm{m}}i}};x,y} \right){\rm{d}}s(x,y) \tag{A2} $$ $$ \partial {\boldsymbol{G}}/\partial{\boldsymbol{ n }}= \left[ {\begin{array}{*{20}{l}} {{\rm{d}}{c_{{\rm{M}}1,1}}}&{{\rm{d}}{c_{{\rm{M}}1,2}}}&{{\rm{d}}{c_{{\rm{M}}1,3}}}&{{\rm{d}}{c_{{\rm{M}}1,4}}} \\ {{\rm{d}}{c_{{\rm{M}}2,1}}}&{{\rm{d}}{c_{{\rm{M}}2,2}}}&{{\rm{d}}{c_{{\rm{M}}2,3}}}&{{\rm{d}}{c_{{\rm{M}}2,4}}} \\ {{\rm{d}}{c_{{\rm{M}}3,1}}}&{{\rm{d}}{c_{{\rm{M}}3,2}}}&{{\rm{d}}{c_{{\rm{M}}3,3}}}&{{\rm{d}}{c_{{\rm{M}}3,4}}} \\ {{\rm{d}}{c_{{\rm{M}}4,1}}}&{{\rm{d}}{c_{{\rm{M}}4,2}}}&{{\rm{d}}{c_{{\rm{M}}4,3}}}&{{\rm{d}}{c_{{\rm{M}}4,4}}} \end{array}} \right] \tag{A3} $$ $$ {{{\rm{d}}}}{c_{{\rm{M}}i,j}} = \left\{ {\begin{array}{*{20}{l}} {\displaystyle\int_{{S_{{\rm{b}}j}}} {\frac{{\partial {G_{\rm{H}}}\left( {{x_{{\rm{m}}i}},{y_{{\rm{m}}i}};x,y} \right)}}{{\partial {\boldsymbol{n}}}}} \cdot {\rm{d}}s(x,y) - \lambda \left( {{x_{{\rm{m}}i}},{y_{{\rm{m}}i}}} \right),{\text{ }}i = j} \\ {\displaystyle\int_{{S_{{\rm{b}}j}}} {\frac{{\partial {G_{\rm{H}}}\left( {{x_{{\rm{m}}i}},{y_{{\rm{m}}i}};x,y} \right)}}{{\partial {\boldsymbol{n}}}}} \cdot {\rm{d}}s(x,y),{\text{ }}i \ne j} \end{array}} \right. \tag{A4} $$ $$ {{\boldsymbol{G}}_{{\rm{mF}}}} = \left[ {\begin{array}{*{20}{c}} {{c_{{\rm{FM}}1,1}}}&{{c_{{\rm{FM}}1,2}}}& \cdots &{{c_{{\rm{FM}}1,{n_{\rm{F}}}}}} \\ {{c_{{\rm{FM}}2,1}}}&{{c_{{\rm{FM}}2,2}}}& \cdots &{{c_{{\rm{FM}}2,{n_{\rm{F}}}}}} \\ \vdots & \vdots & \ddots & \vdots \\ {{c_{{\rm{FM}}4,1}}}&{{c_{{\rm{FM}}4,2}}}& \cdots &{{c_{{\rm{FM}}4,{n_{\rm{F}}}}}} \end{array}} \right] \tag{A5} $$ $$ {c_{{\rm{FM}}i,j}} = \int_{{S_{{\rm{F}}j}}} {{G_{\rm{H}}}} \left( {{x_{{\rm{m}}i}},{y_{{\rm{m}}i}};x,y} \right){\rm{d}}s(x,y) \tag{A6} $$ $$ {{\boldsymbol{G}}_{\rm{F}}} = \left[ {\begin{array}{*{20}{c}} {{c_{{\rm{MF}}1,1}}}&{{c_{{\rm{MF}}1,2}}}& \ldots &{{c_{{\rm{MF}}1,4}}} \\ {{c_{{\rm{MF}}2,1}}}&{{c_{{\rm{MF}}2,2}}}& \ldots &{{c_{{\rm{MF}}2,4}}} \\ \vdots & \vdots & \ddots & \vdots \\ {{c_{{\rm{MF}}{n_F},1}}}&{{c_{{\rm{MF}}{n_F},2}}}& \ldots &{{c_{{\rm{MF}}{n_F},4}}} \end{array}} \right] \tag{A7} $$ $$ {c_{{\rm{MF}}i,j}} = \int_{{S_{{\rm{b}}j}}} {{G_{\rm{H}}}} \left( {{x_{{\rm{F}}i}},{y_{{\rm{F}}i}};x,y} \right){\rm{d}}s(x,y) \tag{A8}$$ $$ \partial {{\boldsymbol{G}}_{\rm{F}}}/\partial {\boldsymbol{n}} = \left[ {\begin{array}{*{20}{c}} {{\rm{d}}{c_{{\rm{MF}}1,1}}}&{{\rm{d}}{c_{{\rm{MF}}1,2}}}& \cdots &{{\rm{d}}{c_{{\rm{MF}}1,4}}} \\ {{\rm{d}}{c_{{\rm{MF}}2,1}}}&{{\rm{d}}{c_{{\rm{MF}}2,2}}}& \cdots &{{\rm{d}}{c_{{\rm{MF}}2,4}}} \\ \vdots & \vdots & \ddots & \vdots \\ {{\rm{d}}{c_{{\rm{MF}}{n_{\rm{F}}},1}}}&{{\rm{d}}{c_{{\rm{MF}}{n_{\rm{F}}},2}}}& \cdots &{{\rm{d}}{c_{{\rm{MF}}{n_{\rm{F}}},4}}} \end{array}} \right] \tag{A9} $$ $$ {{{\rm{d}}}}{c_{{\rm{MF}}i,j}} = \int_{{S_{{\rm{b}}j}}} {\frac{{\partial {G_{\rm{H}}}\left( {{x_{{\rm{F}}i}},{y_{{\rm{F}}i}};x,y} \right)}}{{\partial {\boldsymbol{n}}}}} \cdot {\rm{d}}s(x,y) \tag{A10} $$ $$ {{\boldsymbol{G}}_{{\rm{FF}}}} = \left[ {\begin{array}{*{20}{c}} {{c_{{\rm{FF}}1,1}}}&{{c_{{\rm{FF}}1,2}}}& \ldots &{{c_{{\rm{FF}}1,{n_F}}}} \\ {{c_{{\rm{FF}}2,1}}}&{{c_{{\rm{FF}}2,2}}}& \ldots &{{c_{{\rm{FF}}2,{n_F}}}} \\ \vdots & \vdots & \ddots & \vdots \\ {{c_{{\rm{FF}}{n_F},1}}}&{{c_{{\rm{FF}}{n_F},2}}}& \cdots &{{c_{{\rm{FF}}{n_F},{n_F}}}} \end{array}} \right] \tag{A11} $$ $$ {c_{{\rm{FF}}i,j}} = \int_{{S_{{\rm{F}}j}}} {{G_{\rm{H}}}} \left( {{x_{{\rm{F}}i}},{y_{{\rm{F}}i}};x,y} \right){\rm{d}}s(x,y) \tag{A12} $$ -
![]()
图 2 双孔系统EDFM中MBEM流量项处理的连接关系示例
Figure 2. Example of connection list for MBEM handling flux in EDFM with dual-porosity system
![]()
图 5 一维模型数值解与解析解的对比
Figure 5. Comparison between numerical solution and analytical solution of one dimensional model
![]()
图 8 不同数值解与解析解的对比
Figure 8. Comparison of different numerical solutions and analytical solution
![]()
图 11 不同裂缝倾角无因次应力强度因子数值解与解析解的对比
Figure 11. Comparison of numerical and analytical solutions for dimensionless SIFs with different fracture dip angles
![]()
图 13 两种数值模拟器的网格剖分对比
Figure 13. Comparison of mesh generation between two numerical simulators
![]()
图 14 两种数值模拟器开发1000 d的压力场对比(单位: MPa)
Figure 14. Comparison of pressure distribution maps over 1000 days of two simulators (unit: MPa)
![]()
图 15 两种数值模拟器开发1000 d的X-位移场对比(单位: mm)
Figure 15. Comparison of X-displacement distribution maps over 1000 days of two simulators (unit: mm)
![]()
图 16 两种数值模拟器开发1000 d的Y-位移场对比(单位: mm)
Figure 16. Comparison of Y-displacement distribution maps over 1000 days of two simulators (unit: mm)
![]()
图 18 考虑与不考虑流固耦合作用的生产动态曲线对比
Figure 18. Comparison of production dynamic curveswith and without flow-geomechanics coupled effect
表 1 致密油藏基本参数
Table 1 Basic parameters of tight oil reservoir
Properties Value Properties Value reservoir scale/m 750×340 SRV scale/m 590×220 initial reservoir pressure/MPa 25 reservoir thickness/m 10 matrix porosity 0.12 matrix permeability/mD 0.1 matrix Young’s modulus/GPa 50 matrix Poisson’s ratio 0.3 matrix Biot’s coefficient 0.8 initial water saturation 0.3 fracture porosity 0.4 fracture permeability/mD 5000 initial fracture aperture/m 0.001 fracture proppant Young's modulus/MPa 20 
下载: 导出CSV
-
[1] 孙龙德, 邹才能, 贾爱林等. 中国致密油气发展特征与方向. 石油勘探与开发, 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) [2] 姚军, 孙致学, 张凯等. 非常规油气藏开采中的工程科学问题及其发展趋势. 石油科学通报, 2016, 01: 128-142 (Yao Jun, Sun Zhixue, Zhang Kai, et al. Scientific engineering problems and development trends in unconventional oil and gas reservoirs. Petroleum Science Bulletin, 2016, 01: 128-142 (in Chinese) [3] 王友净, 宋新民, 田昌炳等. 动态裂缝是特低渗透油藏注水开发中出现的新的开发地质属性. 石油勘探与开发, 2015, 42(2): 222-228 (Wang Youjing, Song Xinmin, Tian Changbing, et al. Dynamic fractures are an emerging new development geological attribute in water-flooding development of ultra-low permeability reservoirs. Petroleum Exploration and Development, 2015, 42(2): 222-228 (in Chinese) doi: 10.11698/PED.2015.02.12 [4] 张景, 王颖, 范希彬等. 致密砂砾岩储层应力敏感性评价研究. 中国海上油气, 2020, 32(3): 105-110 (Zhang Jing, Wang Ying, Fan Xibin, et al. Evaluation study on the stress sensitivity of tight glutenite reservoir. China Offshore Oil and Gas, 2020, 32(3): 105-110 (in Chinese) [5] Terzaghi K. Theoretical Soil Mechanics. New York: Wiley online Science, 1943: 9-49
[6] Biot MA. General theory of three-dimensional consolidation. Journal of Applied Physics, 1941, 12(2): 155-164 doi: 10.1063/1.1712886
[7] Biot MA. Theory of elasticity and consolidation for a porous anisotropic solid. Journal of Applied Physics, 1955, 26(2): 182-185 doi: 10.1063/1.1721956
[8] Barenblatt G, Zheltov I. Basic flow equations for homogeneous fluids in naturally fractured rocks. Dokl Akad Nauk, 1960, 13: 545-548
[9] Bai M, Abousleiman Y, Cui L, et al. Dual-porosity poroelastic modeling of generalized plane strain. International Journal of Rock Mechanics and Mining Sciences, 1999, 36(8): 1087-1094
[10] Wilson RK, Aifantis EC. On the theory of consolidation with double porosity. International Journal of Engineering Science, 1982, 20(9): 1009-1035 doi: 10.1016/0020-7225(82)90036-2
[11] Beskos DE, Aifantis EC. On the theory of consolidation with double porosity—II. International Journal of Engineering Science, 1986, 24(11): 1697-1716 doi: 10.1016/0020-7225(86)90076-5
[12] Khaled MY, Beskos DE, Aifantis EC. On the theory of consolidation with double porosity—III A finite element formulation. International Journal for Numerical and Analytical Methods in Geomechanics, 1984, 8(2): 101-123 doi: 10.1002/nag.1610080202
[13] Lewis RW, Schrefler BA. The Finite Element Method in the Deformation and Consolidation of Porous Media. New York: John Wiley & Sons, 1987: 1-337
[14] Jiang JM, Younis RM. Numerical study of complex fracture geometries for unconventional gas reservoirs using a discrete fracture-matrix model. Journal of Natural Gas Science & Engineering, 2015, 26: 1174-1186
[15] Jiang JM, Yang J. Coupled fluid flow and geomechanics modeling of stress-sensitive production behavior in fractured shale gas reservoirs. International Journal of Rock Mechanics & Mining Sciences, 2018, 101: 1-16
[16] Rao X, Cheng LS, Cao RY, et al. An efficient three-dimensional embedded discrete fracture model for production simulation of multi-stage fractured horizontal well. Engineering Analysis with Boundary Elements, 2019, 106: 473-492
[17] Li L, Lee SH. Efficient field-scale simulation of black oil in a naturally fractured reservoir through discrete fracture networks and homogenized media. SPE Res. Eval. & Eng., 2008, 11(4): 750-758
[18] Wang B, Fidelibus C. An open-source code for fluid flow simulations in unconventional fractured reservoirs. Geosciences (Switzerland)
, 2021, 11(2): 106 [19] Melenk JM, Babuka I. The partition of unity finite element method: Basic theory and applications. Computer Methods in Applied Mechanics and Engineering, 1996, 139(1-4): 289-314 doi: 10.1016/S0045-7825(96)01087-0
[20] Moës N, Dolbow J, Belytschko T. A finite element method for crack growth without remeshing. International Journal for Numerical Methods in Engineering, 1999, 46: 131-150 doi: 10.1002/(SICI)1097-0207(19990910)46:1<131::AID-NME726>3.0.CO;2-J
[21] Ren GT, Jiang JM, Younis RM. A fully coupled XFEM-EDFM model for multiphase flow and geomechanics in fractured tight gas reservoirs. Procedia Computer Science, 2016, 80: 1404-1415 doi: 10.1016/j.procs.2016.05.449
[22] Ren GT, Jiang JM Younis RM. Fully-coupled XFEM-EDFM hybrid model for geomechanics and flow in fractured reservoirs//SPE Reservoir Simulation Conference, Montgomery, Texas, USA, 2017. Paper No. SPE-182726-MS
[23] Yan X, Huang ZQ, Yao J, et al. An efficient hydro-mechanical model for coupled multi-porosity and discrete fracture porous media. Computational Mechanics, 2018, 62: 943-962 doi: 10.1007/s00466-018-1541-5
[24] Zeng QD, Yao J, Shao JF. Study of hydraulic fracturing in an anisotropic poroelastic medium via a hybrid EDFM-XFEM approach. Computers and Geotechnics, 2019, 105: 51-68
[25] Ren G, Jiang J, Younis RM. A model for coupled geomechanics and multiphase flow in fractured porous media using embedded meshes. Advances in Water Resources, 2018, 122: 113-130
[26] Jansen G, Benoît V, Miller SA. THERMAID—a matlab package for thermo-hydraulic modeling and fracture stability analysis in fractured reservoirs, 2018, arXiv:1806.10942
[27] Lee SH, Lough MF, Jensen CL. Hierarchical modeling of flow in naturally fractured formations with multiple length scales. Water Resour. Res., 2001, 37(3): 443-455 doi: 10.1029/2000WR900340
[28] Praditia T, Helmig R, Hajibeygi H. Multiscale formulation for coupled flow-heat equations arising from single-phase flow in fractured geothermal reservoirs. Computational Geosciences, 2018, 22(5): 1305-1322 doi: 10.1007/s10596-018-9754-4
[29] Fang SD, Cheng LS, Ayala LF. A coupled boundary element and finite element method for the analysis of flow through fractured porous media. Journal of Petroleum Science and Engineering, 2017, 152: 375-390 doi: 10.1016/j.petrol.2017.02.020
[30] Cao RY, Fang SD, Jia P, et al. An efficient embedded discrete-fracture model for 2D anisotropic reservoir simulation. Journal of Petroleum Science and Engineering, 2018, 174: 115-130
[31] Peaceman DW. Interpretation of well-block pressures in numerical reservoir simulation. Society of Petroleum Engineers Journal, 1978, 18(3): 183-194 doi: 10.2118/6893-PA
[32] Bai M. On equivalence of dual-porosity poroelastic parameters. Journal of Geophysical Research, 1999, 104(B5): 10461-10466 doi: 10.1029/1999JB900072
[33] Moinfar A, Sepehrnoori K, Johns RT. et al. Coupled geomechanics and flow simulation for an embedded discrete fracture model//SPE Reservoir Simulation Symposium, the Woodlands, Texas, 2013
[34] Blasingame TA, Poe BD. Semianalytic solutions for a well with a single finite-conductivity vertical fracture. Paper SPE, 1993. Paper No. SPE-163666-MS
[35] Tada H, Paris PC, Irwin GR. The Stress Analysis of Cracks Handbook. New York: ASME Press, 2000: 10-100
计量
- 文章访问数: 1364
- HTML全文浏览量: 506
- PDF下载量: 132
