Characterization of Key Tight Oil Parameters and Mass Transfer of Counter-Current Imbibition in Fractured Tight Oil Reservoirs
-
摘要: 在致密油藏水平井体积压裂开采过程中，压裂液通过缝网与基质接触并发生逆向渗吸作用，由于接触面积很大渗吸作用不可忽视；但目前关于表征致密储层的渗吸作用，从而研究渗吸对水平井体积压裂生产过程影响的研究尚未深入.为了解决以上问题，首先利用毛管束模型，通过考虑致密储层中边界层的特征，建立了解析的渗流参数计算表达式，用以计算致密储层的渗透率、毛管力、相渗曲线这3个关键渗流参数；同时，基于以上关键渗流参数和渗吸控制方程建立了适用于致密储层的渗吸速度计算模型；然后，将渗吸项作为源汇项加入到考虑缝网的双孔单渗模型中.最后，在真实水平井体积压裂开采过程中，耦合渗吸作用.研究表明，相比于不考虑边界层特征的致密油藏，边界层的存在将大幅度减弱储层的渗吸能力，同时也说明了在致密储层中，边界层的存在是不可忽视的，如果在渗吸计算中忽视致密储层的边界层特征会严重高估渗吸对致密储层产能的影响.Abstract: During the process of tight oil exploration, counter current imbibition effect is significantly different due to the presence of complex fracture network and flow characteristics in tight oil reservoirs. But at present, there is no proper model to simulate counter-current imbibition in fractured reservoir considering the characteristics of tight oil formation. In order to solve this problem, PEBI grids are used to match the complex fracture network, natural fractures and matrix are idealized as dual-porosity medium, rate of mass transfer of imbibition between matrix and fractures is treated as source or sink term in dual porosity model. A new semi analytical model for the calculation of mass transfer function of counter-current imbibition in the presence of complex fracture network is established by using radial integration boundary element method (RIBEM). In addition, to reflect the flow characteristics of tight oil, relative permeability and capillary pressure curve with the effect of boundary layer considered, and mixed wettability have also used in the mass transfer model. Besides, we show the capacity and practical application of the model with a field example from tight oil reservoir. From simulated results, it is concluded that counter-current imbibition plays an important role in improving the oil recovery and the existence of boundary layer reduces the contribution of imbibition to oil production dramatically.
-
0. 引言
根据邹才能等(2012)的研究成果：我国致密储层由两个关键参数定义，孔隙度小于10%，孔喉直径小于1 μm或是渗透率小于1 mD的储层可称之为致密储层.在致密储层的微纳米尺度孔喉中，存在着一种不容忽视的现象——边界层现象(李洋等，2011)，边界层是在多孔介质内壁上黏附的一层不可动润湿相流体，因为边界层在致密储层的微纳米喉道中所占的比例不能忽略，导致流体在致密储层中渗流时的真实流动半径小于储层原始的喉道半径(田虓丰等，2016).
与北美致密储层相比，我国鄂尔多斯、长庆等油田的致密储层，渗透率更低，含油饱和度更低(吴松涛等，2015；任大忠等，2016)，从孔喉分布概率分布曲线上看，微纳米孔喉所占的比例更高，意味着我国致密储层受边界层的影响更加严重.
渗吸是在毛管力作用下，岩石吸入润湿相流体排出非湿相流体的过程，对于致密储层水平井周围的SRV(stimulated reservoir volume)区域，无论是衰竭开采还是水吞吐或注水开采方式，都会有大量的润湿相流体通过缝网进入地层与基质大面积接触，并且在毛管力的作用下发生大面积的渗吸；同时，基于毛管力的表征公式：
$$ {p_{\rm{c}}} = \frac{{2\delta \cos \theta }}{r}, $$ (1) 其中，pc为毛管压力，10-1 MPa；δ为表面张力，N/m；θ为接触角，度；r为毛细管半径，μm.
由于致密储层的微小喉道半径，导致致密储层渗吸过程毛管力较普通油藏大，具有更强的渗吸驱动力，综合以上两点，在致密储层开发过程中，渗吸作用不容忽视.
目前关于致密储层渗吸速度的研究，从实验方面，当前学者主要使用更为先进的手段，如运用核磁共振技术(韦青等，2016)来监测致密储层的渗吸过程，得到了致密岩心的渗吸流动下限、喉道半径以及渗吸过程的主流喉道半径，但考虑致密储层流动特征的渗吸速度表达式尚未得出；而从理论方面，当前渗吸速度的计算主要来源于3类方法，第1类是以Aronofsky et al.(1958)、Ma et al.(1997)为代表的标度模型方法，第2类是以Kashchiev and Firoozabadi(2003)为代表的解析方法，第3类是以Hassan(2004)、Bagherinezhad and Pishvaie(2014)为代表的数值模拟方法；第3类数值方法，克服了第1类方法理论依据薄弱以及第2类方法假设条件过于苛刻不适用于实际渗吸模拟的缺点，另外为了表征裂缝性油藏的渗吸作用，Xu et al.(2017)通过引入径向积边界元方法，提出了简单高效且适用于双重介质的渗吸速度计算模型，同时分形模型在计算渗吸速度方面也取得了较好的效果(蔡建超和郁伯铭，2012；蔡建超等，2013)，然而当前所有的渗吸速度计算模型都没有考虑致密储层的流动特征，因此，总的来看当前缺乏考虑致密储层流动特征的渗吸速度计算模型.
为了解决以上关于表征致密储层流动特征的问题，本文首先基于微纳米喉道中的边界层理论，表征了致密储层中具有有效流动能力的喉道半径(有效喉道半径)特征，在此基础上利用毛管束模型，提出了反映致密油藏流体真实渗流能力的有效渗透率模型；其次，利用统计得到的致密储层毛管压力模型，通过考虑渗吸过程饱和度的变化顺序，提出了致密储层的毛管压力计算模型；再次，利用已建立的有效渗透率模型，提出了致密储层的相对渗透率计算方法；同时，将以上得到的关键渗吸参数，带入到已有的渗吸速度计算模型中，得到了关于致密储层渗吸速度的表征模型，并利用此渗吸速度表征模型，通过与实验对比的方法验证了本文提出的关键参数计算方法的正确性；最后，通过实例计算阐明了本模型在宏观数值模拟过程中的具体用法，同时通过模拟，研究了渗吸对致密储层产能以及油水分布的影响.
1. 致密储层渗吸机理及模型
1.1 致密储层中的渗吸机理
在渗吸过程中，若水的吸入方向与油的排出方向相同，则为顺向渗吸，否则为逆向渗吸，两种渗吸机理依赖于边界条件和重力与毛管力的比值Bond数的倒数，即NB-1(Schechter et al., 1994)：
$$ N_{\rm{B}}^{ - 1} = C\frac{\delta }{{\Delta \rho gH}}\sqrt {\frac{\varphi }{K}} , $$ (2) 其中C为结构常数，无因次；δ为表面张力，N/m；Δρ为油水密度差，kg/m3；g为重力加速度，取9.8 m/s2；H为储层厚度，m；φ为孔隙度，无因次；K为渗透率，mD，NB-1反映基质渗吸过程中毛细管力和重力贡献的相对大小及渗吸机理.当NB-1较大时，毛管力起支配作用.当NB-1趋于零时，重力起支配作用；在致密储层中，因重力作用引起的垂向上的流动可以忽略；同时，由于致密储层中不可忽视边界层的存在，由公式(1) 可以知道，致密储层渗吸过程毛管力比普通中低渗油藏大，所以在致密储层中毛管力起绝对支配作用；另外，如图 1所示，在致密储层体积压裂的开发方式下，裂缝极为发育，极大地增加了基质块与裂缝中润湿相流体的接触面积；综合以上因素，致密储层开发过程中发生着大规模的自发逆向渗吸.
图 1 存在复杂缝网的逆向渗吸Fig. 1. Counter-current imbibition in the presence of fracture network1.2 致密储层特征及渗吸流动特征
通过对我国鄂尔多斯盆地致密岩心的统计，鄂尔多斯致密岩心的喉道分布范围在20 nm~1.2 μm之间，绝大部分的喉道半径分布都在20.0 nm~0.7 μm之间，渗透率数量级大小在10-1 mD.在这种微纳米尺度的喉道中，边界层的存在不容忽视，边界层对致密储层中渗吸过程的影响如图 2所示，相比于普通中低渗油藏喉道中的渗吸过程(图 2a)，致密油藏喉道中因存在边界层导致了有效流动喉道半径reff减小(图 2b)，在极端情况下，当边界层厚度h等于原始喉道半径rori时，边界层的存在甚至会导致喉道内的所有流体都变成不可动流体.为了定量表征边界层对流动的影响，根据Meng et al.(2017)以及田虓丰(2015)的研究，可以定量表征边界层厚度，其表达式为：
图 2 边界层对致密储层流动的影响a.存在边界层时流体在喉道中的流动；b.不存在边界层时流体在喉道中的流动Fig. 2. The influence of boundary layer on the flow of tight throat$$ \begin{array}{l} h = \\ \left\{ \begin{array}{l} r0.25763{{\rm{e}}^{ - 0.261r}}{\left( {\nabla p} \right)^{ - 0.419}}\mu ,\nabla p < 1{\rm{MPa/m;}}\\ r0.25763{{\rm{e}}^{ - 0.261r}}\mu ,\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\nabla p > 1{\rm{MPa/m;}} \end{array} \right. \end{array} $$ (3) 其中，r为孔隙的喉道半径，μm；h为孔隙中的边界层厚度，μm.
需要注意的是公式(3) 是通过实验拟合得到的，具有一定的适用范围，应用时喉道壁面需完全润湿，压力梯度不得小于0.02 MPa，同时喉道半径应大于0.2 μm，管内的流体需是牛顿流体.
1.3 考虑缝网的渗吸速度计算模型
在微观上致密储层受边界层的影响，而从宏观上看，如图 1所示，在致密储层的水平井分段多簇体积压裂开发方式下，近水平井地带会形成复杂的缝网改造区，在表征致密储层中的流体流动以及渗吸过程时都不可避免地需要考虑缝网的存在.Jia et al.(2015)、方思冬等(2017)在致密储层缝网的相关研究中认为，改造区内的缝网系统以导流能力较强的人工压裂缝、导流能力较弱的次生裂缝及天然裂缝纵横交织形成，在数值模拟处理过程中，水力压裂缝因其导流能力较强，需对其显式描述；导流能力较弱的次生裂缝和自然裂缝因其数目众多，准确地描述这些为数众多的次生裂缝以及天然裂缝既不现实同时也降低了计算的效率，因此通常使用双孔单渗介质近似描述次生裂缝、天然缝以及致密储层的基质.
在双孔单渗模型的基础上，基于Xu et al.(2017)的研究，考虑逆向渗吸作用时，整个裂缝-基质系统的流动可由式(4)~(6) 描述：
$$ 裂缝介质中油相:\nabla \cdot \left( {K\nabla {p_{\rm{o}}}} \right) + {q_{\rm{w}}} + {q_{{\rm{imo}}}} = {\mu _{\rm{o}}}{C_{\rm{t}}}\varphi \frac{{\partial {p_{\rm{o}}}}}{{\partial t}}, $$ (4) $$ 裂缝介质中水相:\nabla \cdot \left( {K\nabla {p_{\rm{w}}}} \right) + {q_{\rm{w}}} + {q_{{\rm{imw}}}} = {\mu _{\rm{w}}}{C_{\rm{t}}}\varphi \frac{{\partial {p_{\rm{w}}}}}{{\partial t}}, $$ (5) $$ 基质块中渗吸:\frac{\partial }{{\partial {x_i}}}\left( {D\left( {{S_{\rm{w}}}} \right)\frac{{\partial {S_{\rm{w}}}}}{{\partial {x_i}}}} \right) = \varphi \frac{{\partial {S_{\rm{w}}}}}{{{\partial _t}}}, $$ (6) 其中，K为天然裂缝的渗透率，mD；po，pw为油、水相的压力，10-1 MPa；qo，qw为源汇相，m3/d；qimo，qimw为渗吸速度项，m3/d；μo，μw为油相和水相的粘度，mPa·s；φ为孔隙度，无因次；Ct为综合压缩系数，1/MPa；t为时间，d；D(Sw)为毛管力扩散系数，m2/s.
笔者将自发逆向渗吸引起的基质与裂缝间的流体交换作为源汇项qimo和qimw分别引入到裂缝系统中的油相控制方程(式(4) 和水相控制方程(式(5))中.从式(4) 和(5) 可以看出，只要准确地描述了渗吸项qimo和qimw，便可以搞清渗吸在体积压裂开发过程中的影响.
自发逆向渗吸的控制方程如式(6) 所示，渗吸速度由式(6) 中的毛管力扩散系数D(Sw)决定，D(Sw)具体形式如式(7) 所示：
$$ D\left( {{S_{\rm{w}}}} \right) = \varphi \left( {\frac{1}{{{\lambda _{\rm{o}}}}} + \frac{1}{{{\lambda _{\rm{w}}}}}} \right){\left( {\frac{{\partial {{p'}_{\rm{c}}}}}{{\partial {S_{\rm{w}}}}}} \right)^{ - 1}}, $$ (7) $$ {\lambda _{\rm{o}}} = \frac{{k{k_{{\rm{ro}}}}}}{{{\mu _{\rm{o}}}}},{\lambda _{\rm{w}}} = \frac{{k{k_{{\rm{rw}}}}}}{{{\mu _{\rm{w}}}}}, $$ (8) 其中，λo、λw分别为油相与水相的流度，mD/mpa·s；kro、krw分别为油相与水相的相对渗透率，小数.
从毛管力扩散系数的表达式(7) 可以看出，基质块中油水的流度以及毛管力曲线共同决定着毛管力扩散系数，而流度项又由渗透率、相对渗透率共同决定.因此为了计算致密储层的渗吸特征，只要得到致密油储层中考虑边界层的绝对渗透率、相对渗透率以及毛管力曲线，便可以得到相应致密油储层的毛管力扩散系数；而渗吸速度的计算则可以采用前文提到的Xu et al.(2017)的径向积边界元方法(radial integration boundary element method，简称RIBEM)模型来计算，现在问题的关键便是合理计算致密储层的绝对渗透率、相对渗透率以及毛管力曲线.
2. 致密储层关键渗流参数的计算
本文基于毛管束模型，建立了致密储层关键渗流参数的解析表达式，并为致密储层渗吸模型的提出提供基础.所用的毛管束模型遵循的基本假设如下：(1) 毛管束模型由半径从rmin到rmax连续分布的毛细管组成，毛细管长度等于毛管束模型的长度；(2) 流体微可压缩；(3) 流体粘度不变；(4) 单根喉道中的流动符合泊肃叶定律；(5) 毛管束模型的流动符合达西公式；(6) 每根毛细管中只有一种流体；(7) 流动是等温油水两相渗流.
2.1 致密储层渗透率的计算
如图 2b所示，单根喉道内, 考虑边界层以后，喉道内有效喉道半径变小，考虑边界层的毛细管层流流量为：
$$ {q_i} = \frac{{{\rm{\pi }}r_i^4{{\left( {1 - \frac{{{h_i}}}{{{r_i}}}} \right)}^4}\nabla p}}{{8\mu }}, $$ (9) 其中，qi为第i根喉道的流量，10-9 cm3/s；ri为第i根喉道半径，μm；∇p为压力梯度，MPa/m；μ为流体粘度，mPa·s.
将毛管束模型中所有的圆管层流流量进行累加，得到毛管束模型的总流量为：
$$ Q = \sum\limits_{{r_{\rm{c}}}}^{{r_{\max }}} {\frac{{{n_i}{\rm{\pi }}r_i^4{{\left( {1 - \frac{{{h_i}}}{{{r_i}}}} \right)}^4}\nabla p}}{{8\mu }}} , $$ (10) 其中，Q为毛管束流量，10-9 cm3/s；rmax为最大喉道半径，μm；rc为临界喉道半径，μm.由于压汞实验终止时的压力不一定能够满足最小喉道的压力，因此当rc＜rmin时，rc=rmin；ni为半径为ri的喉道数量，无因次.
因为毛管束模型的基本假设是毛细管是直管，因此毛管束上任意截面的面孔率与孔隙度相等，半径为ri的喉道所占体积比例与面积比例相等.将式(10) 与达西定律(式(11))联立求解得到有效渗透率Keff：
$$ Q = \frac{{KA\nabla p}}{\mu }, $$ (11) $$ {K_{{\rm{eff}}}} = \sum\limits_{{r_{\rm{c}}}}^{{r_{\max }}} {\frac{{\tau \varphi {S_i}r_i^2{{\left( {1 - \frac{{{h_i}}}{{{r_i}}}} \right)}^2}}}{8}} , $$ (12) 式(12) 中的Si和φ分别由式(13) 和式(14) 表示：
$$ {S_i} = \frac{{{n_i}{\rm{\pi }}r_i^2}}{{\sum\limits_{{r_{{\rm{min}}}}}^{{r_{\max }}} {{n_i}{\rm{\pi }}r_i^2} }}, $$ (13) $$ \varphi = \frac{{\sum\limits_{{r_{\rm{c}}}}^{{r_{\max }}} {{n_i}{\rm{\pi }}r_i^2} }}{A}, $$ (14) 其中，Keff为有效渗透率，D；A为毛管束截面积，cm2；Si为半径为ri的喉道占岩心孔隙体积的比例，无因次.
使用毛管力与半径的关系式(1)，将式(12) 变为式(15)：
$$ {K_{{\rm{eff}}}} = \frac{{\tau \varphi {{\left( {{\delta _{{\rm{Hg}}}}\cos {\theta _{{\rm{Hg}}}}} \right)}^2}}}{2}\sum\limits_{{r_{\rm{e}}}}^{{r_{\max }}} {\frac{{{S_i}}}{{p_{{\rm{cHg}}i}^2}}} , $$ (15) 其中：
$$ {p_{{\rm{cHg}}i}} = \frac{{2{\delta _{{\rm{Hg}}}}\cos {\theta _{{\rm{Hg}}}}}}{{{r_i} - {h_i}}}, $$ (16) 其中，pc为毛管压力，MPa；δ为表面张力，N/m；θ为接触角，度；r为毛细管半径，μm；δHg为汞气表面张力，N/m；θHg为汞气接触角，度；pcHgi为半径为ri的毛管的进汞毛管压力，MPa.
对于任意毛细管，式(15) 中δcosθ/pcHg相同，只考虑累加项即可，与流体无关；因此毛管力一直采用汞气毛管力，将进汞过程分为M个阶段，在第nj个阶段，压力由pcHg1上升到pcHg2，进汞饱和度由SHg1增加到SHg2，进汞过程由第n根毛管到第m根毛管，根据式(15)，在nj阶段内，渗透率可以表示为：
$$ {K_{{\rm{eff}}}} = \frac{{\tau \varphi {{\left( {{\delta _{{\rm{Hg}}}}\cos {\theta _{{\rm{Hg}}}}} \right)}^2}}}{2}\sum\limits_{i = n}^m {\frac{{{S_i}}}{{p_{{\rm{cHg}}i}^2}}} , $$ (17) 其中，n为第nj个阶段进汞开始时的喉道编号，m为第nj个阶段进汞结束时的喉道编号.
对式(17) 中的累加部分进行变换得：
$$ \begin{array}{l} \sum\limits_{i = n}^m {\frac{{{S_i}}}{{p_{{\rm{cHg}}i}^2}}} = \frac{{{S_{{\rm{Hg}}n}} - {S_{{\rm{Hg}}n - 1}}}}{{p_{{\rm{cHg}}n}^2}} + \frac{{{S_{{\rm{Hg}}n + 1}} - {S_{{\rm{Hg}}n}}}}{{p_{{\rm{cHg}}n + 1}^2}} + \cdots + \\ \frac{{{S_{{\rm{Hg}}m - 1}} - {S_{{\rm{Hg}}m - 2}}}}{{p_{{\rm{cHg}}m - 1}^2}} + \frac{{{S_{{\rm{Hg}}m}} - {S_{{\rm{Hg}}m - 1}}}}{{p_{{\rm{cHg}}m}^2}} = \\ \frac{{\Delta {S_{{\rm{Hg}}n}}}}{{p_{{\rm{cHg}}i}^2}} + \frac{{\Delta {S_{{\rm{Hg}}n + 1}}}}{{p_{{\rm{cHg}}n + 1}^2}} + \cdots + \frac{{\Delta {S_{{\rm{Hg}}m - 1}}}}{{p_{{\rm{cHg}}m - 1}^2}} + \frac{{\Delta {S_{{\rm{Hg}}m}}}}{{p_{{\rm{cHg}}m}^2}}, \end{array} $$ (18) 其中，SHgi为压力pcHgi时的汞饱和度，小数；ΔSHgi为压力由pcHgi-1变为pcHgi时的进汞饱和度变化，小数.
由于进汞过程中压力在逐渐增加，因此pcHgn＜pcHgm，则式(18) 有如下关系：
$$ \frac{{\Delta {S_{{\rm{Hg}}nj}}}}{{p_{{\rm{cHg}}m}^2}} \le \sum\limits_{i = n}^m {\frac{{{S_i}}}{{p_{{\rm{cHg}}m}^2}}} \le \frac{{\Delta {S_{{\rm{Hg}}nj}}}}{{p_{{\rm{cHg}}n}^2}}, $$ (19) 将式(18) 中累积项简化为：
$$ \sum\limits_{i = n}^m {\frac{{{S_i}}}{{p_{{\rm{cHg}}i}^2}}} \approx \frac{{\Delta {S_{{\rm{Hg}}nj}}}}{{p_{{\rm{cHgave}}}^2}}, $$ (20) 式(20) 及(19) 中，ΔSHgnj为第nj阶段的汞饱和度变化，ΔSHgnj=ΔSHgn+ΔSHgn+1+…+ΔSHgm-1+ΔSHgm.另外式(20) 中：
$$ {p_{{\rm{cHgave}}}} = \frac{{{p_{{\rm{cHg}}n}} + {p_{{\rm{cHg}}m}}}}{2}, $$ (21) 其中，pcHgave为第nj阶段的平均进汞压力，MPa.将式(20) 代入式(15) 得：
$$ {K_{{\rm{eff}}}} \approx \frac{{\tau \varphi {{\left( {{\delta _{{\rm{Hg}}}}\cos {\theta _{{\rm{Hg}}}}} \right)}^2}}}{2}\sum\limits_{nj = 1}^{nj = M} {\frac{{\Delta {S_{{\rm{Hg}}nj}}}}{{p_{{\rm{cHgave}}}^2}}} , $$ (22) 其中，M为整个进汞过程细分的阶段数.
当进汞过程被细分成无穷多阶段时，即M→∞时，ΔSHgnj→0，pcHgave→pcHgn=pcHgm，此时的渗透率不再是近似渗透率(式(17))，而是真实渗透率：
$$ K = \frac{{\tau \varphi {{\left( {{\sigma _{{\rm{Hg}}}}\cos {\theta _{{\rm{Hg}}}}} \right)}^2}}}{2}\int_0^{{S_{{\rm{im}}}}} {\frac{1}{{p_{{\rm{cHg}}}^2}}{\rm{d}}{S_{{\rm{Hg}}}}} , $$ (23) 其中，Sim为不可动流体饱和度，小数.
不可动流体由两部分构成，一是边界层流体，二是不连通孔喉中的流体.毛管束模型中边界层流体饱和度为：
$$ {S_{\rm{b}}} = \sum\limits_{{r_{\min }}}^{{r_{\max }}} {\left[ {{{\left( {\frac{{{h_i}}}{{{r_i}}}} \right)}^2} \cdot {\alpha _i}} \right]} . $$ (24) 压汞曲线中不能进汞的部分是由于毛管力的作用导致汞无法进入微细孔喉，因此在油水渗流的过程这部分孔喉中的流体也无法参与流动，定义这部分流体为不连通孔喉中的流体，则其饱和度为：
$$ {S_{\rm{s}}} = 1 - {S_{{\rm{Hgmax}}}}, $$ (25) 其中，SHgmax为最大进汞饱和度，小数.
将式(24) 和式(25) 相加得到不可动流体饱和度为，由此可以计算岩心的绝对渗透率：
$$ {S_{{\rm{im}}}} = {S_{\rm{s}}} + {S_{\rm{b}}}. $$ (26) 2.2 致密储层毛管力曲线的计算
田虓丰(2015)统计发现，致密储层压汞曲线可以划分为两种模式，模式一与模式二进汞压力与进汞饱和度分别满足关系式(27) 和(28)：
$$ 模式一:\;\;\;\;\;\;\;\;\;{p_{{\rm{cHg}}}} = {p_{{\rm{c1}}}}{{\rm{e}}^{{p_{{\rm{c2}}}}/\left( {1 - {S_{{\rm{Hg}}}}} \right)}}, $$ (27) $$ 模式二:\;\;\;\;\;\;\;\;\;{p_{{\rm{cHg}}}} = {p_{{\rm{c1}}}}{{\rm{e}}^{{p_{{\rm{c2}}}}{S_{{\rm{Hg}}}}}}, $$ (28) 其中，pc1，pc2为常数.
模式一压汞曲线进汞段更加陡峭，喉道分布的非均质性严重，喉道分布范围更大.模式二压汞曲线进汞段更加平直，喉道分布的非均质性较弱，喉道分布范围更集中.以上计算的是汞的毛管压力和饱和度，再使用式(29)，将汞的毛管压力转化为实际的油水毛管压力：
$$ {p_{{\rm{cw}}}} = \frac{{{\sigma _{\rm{w}}}\cos {\theta _{\rm{w}}}}}{{{\sigma _{{\rm{Hg}}}}\cos {\theta _{{\rm{Hg}}}}}}{p_{{\rm{cHg}}}}. $$ (29) 2.3 相对渗透率曲线的计算
将式(27)、式(28) 代入到式(23) 中，得到考虑两相边界层以及混合润湿条件下的相有效渗透率计算公式：
$$ \begin{array}{l} K = \frac{{\tau \cdot \varphi ' \cdot {{\left( {{\sigma _{{\rm{Hg}}}} \cdot \cos {\theta _{{\rm{Hg}}}}} \right)}^2}}}{2} \cdot \\ \int_{{S_{{\rm{nwi}}}}}^{1 - {S_{{\rm{wi}}}}} {\frac{1}{{\left[ {{{p}_{{\rm{cHg}}}}'/\left( {1 - c} \right)} \right]}}{\rm{d}}{S_{{\rm{nw}}}}} , \end{array} $$ (30) 其中，Snwi为不可动非湿相饱和度，小数；Swi为不可动湿相饱和度，小数.
当非湿相饱和度为Snw时，非湿相和湿相有效渗透率为：
$$ {K_{{\rm{nw}}}} = \frac{{\tau \cdot \varphi ' \cdot {{\left( {{\sigma _{{\rm{Hg}}}}\cos {\theta _{{\rm{Hg}}}}} \right)}^2}}}{2} \cdot \int_{{S_{{\rm{nwi}}}}}^{{S_{{\rm{nw}}}}} {\frac{1}{{p_{{\rm{cHg}}}^{'2}}}{\rm{d}}{S_{{\rm{nw}}}}} , $$ (31) $$ {K_{\rm{w}}} = \frac{{\tau \cdot \varphi ' \cdot {{\left( {{\sigma _{{\rm{Hg}}}}\cos {\theta _{{\rm{Hg}}}}} \right)}^2}}}{2} \cdot \int_{{S_{{\rm{nw}}}}}^{1 - {S_{{\rm{wi}}}}} {\frac{1}{{p_{{\rm{cHg}}}^{'2}}}{\rm{d}}{S_{{\rm{nw}}}}} , $$ (32) 其中，Knw为当非湿相饱和度为Snw，考虑边界层以及混合润湿情况下非湿相有效渗透率，103 mD；Kw为当非湿相饱和度为Snw，考虑边界层以及混合润湿情况下湿相有效渗透率，103 mD.
将式(31) 和式(32) 除以式(33) 得到当非湿相饱和度为Snw时，非湿相和湿相相对渗透率为：
$$ {K_{{\rm{nw}}}} = \frac{{{K_{{\rm{nw}}}}}}{K},K = \frac{{{K_{\rm{w}}}}}{K}. $$ (33) 总之，在已知岩心孔喉分布的前提下，通过以上得到的表达式，可以得到致密储层存在边界层情况下的绝对渗透率、相对渗透率以及毛细管压力计算表达式，将以上参数带入式(7) 中，便得到考虑边界层的自发逆向渗吸模型，使用Xu et al.(2017)的RIBEM便可计算基质块的渗吸速度.
3. 参数计算方法的验证与分析
为了验证第2节所提的流动参数计算方法的正确性，首先选用两端敞开四周封闭的致密岩心进行渗吸实验(图 3)，然后分两种情况进行数值模拟，第一种情况为渗吸时考虑边界层的存在，第二种情况为渗吸时不考虑边界层的存在.两种情况都采用Xu et al.(2017)的RIBEM数值模拟方法，最后通过与实验结果对比验证所提流动参数计算方法的正确性.
图 3 致密岩心OEO渗吸实验示意Fig. 3. Schematic of tight core imbibition experiments with OEO boundary conditions岩心基本参数如表 1所示，流动参数的计算步骤如下，首先利用高压压汞实验获得的喉道分布曲线(图 4a不考虑边界层的毛管压力曲线)建立平行毛管束模型；再利用第2节的理论计算出相应的反应致密储层的流动参数，包括毛管力曲线、相渗曲线、以及相应的渗透率；将上述计算好的渗流参数代入毛管力扩散系数表达式(7) 建立相应的自发逆向渗吸数学模型；最后，利用已有的数值模拟方法针对一定边界条件下的致密岩心进行渗吸计算.
表 1 岩心基本参数Table Supplementary Table Basic parameters of imbibition experiment under different boundary conditions岩心编号 岩心直径(mm) 岩心长度(mm) 孔隙度(%) 束缚水饱和度(%) 模拟油密度(kg/m3) 地层水密度(kg/m3) S102 24.96 52.6 9.3 41.1 807 1 058 考虑边界层的影响, 计算得到的毛管力曲线如图 4所示，可以看出当考虑边界层时，其毛管力曲线在整体上要高于不考虑边界层的致密岩心.进一步观察表明，随着含水饱和度的减小毛管力曲线间的差异越来越大，在含水饱和度较大的时候，毛管力曲线间的差别不大，而在含水饱和度较小的时候，毛管力曲线的差别较大；出现这种现象的原因是高含水饱和度区域对应着较大孔喉半径，而在大孔喉中边界层虽然存在，但所占比例较小，因此较大喉道半径的孔喉分布区域中的毛管力曲线相差较小，而在较小喉道半径的孔喉分布区域中，边界层在喉道中所占比例较大，基于公式(1) 可知，在边界层的微纳米喉道中，毛管力明显增加.
考虑边界层的影响, 计算得到的相渗曲线如图 5所示, 可以看出边界层的存在明显减弱了各相流体的流动能力，出现这种现象的原因是：边界层的存在，缩小了油水的流动通道，更容易使得原本迂回联通的流动通道被边界层截断，形成不连续的流动通道，从而减小油相以及水相的流动能力.另外不考虑边界层时其渗透率为7.56 mD，考虑边界层后其渗透率下降到1.28 mD，说明边界层的存在很大程度上减小了岩心的流动能力.
利用以上参数计算得到的相应毛管力扩散系数如图 6所示，从图 6可以明显看出考虑边界层以后，致密岩心的毛管力扩散系数影响范围变窄，渗吸强度也大幅度减弱，综合来说便是边界层的存在大幅度减弱了致密储层的渗吸能力.
最终从计算结果与实验结果的对比来看(图 7)，如果不考虑边界层的影响，计算得到的采收率及渗吸速度要明显高于考虑边界层的致密储层的采收率以及渗吸速度，并且考虑边界层以后，采收率以及渗吸速度与实验结果对比效果较好，说明了本文计算致密储层参数方法的正确性.
图 7 模型计算结果与实验数据对比Fig. 7. Comparison between the calculated results and experimental data4. 实例计算
为了研究在油藏尺度下，渗吸对致密油藏开发效果的影响，采用鄂尔多斯油田的一口真实的多级压裂水平井，水平井筒周围水力压裂缝的展布如图 8a所示，使用PEBI网格贴合人工水力压裂缝对致密油藏进行网格剖分，对应的网格剖分结果如图 8b所示；将第3节中考虑边界层与未考虑边界层的情况，应用到PEBI网格模型中，通过对比不同时刻油藏的含水饱和度和产量，评价渗吸对致密储层开发效果的影响.
本节模拟一口水平井的衰竭开采，地层在模拟之前已被注入压裂夜，初始含水饱和度如图 9a所示，衰竭前已关井20 d，使得储层中基质块与裂缝系统中的压裂液充分接触，发挥逆向渗吸作用，衰竭开采一共持续50 d，与模拟有关的其他参数如表 2所示，最终模拟的饱和度分布结果如图 9b和9c所示，采油速度及累计产量如图 10所示.
图 9 不同情况下饱和度分布a.初始时刻含水饱和度；b.50 d时不考虑边界层的含水饱和度；c.50 d时考虑边界层的含水饱和度Fig. 9. Saturation distribution under the condition of without considering flow characteristics of tight oil reservoir表 2 油藏基本参数Table Supplementary Table Basic parameters of reservoirs参数名称 数值 初始油藏压力(10-1 MPa) 250 井底压力(10-1 MPa) 5 致密油藏基质渗透率(考虑边界层)(mD) 1.28 致密油藏基质渗透率(不考虑边界层)(mD) 7.56 油藏天然裂缝渗透率(mD) 100 水力压裂缝的渗透率(mD) 2 000 基质的孔隙度 0.093 裂缝介质的空隙度 0.001 水力压裂缝的孔隙度 0.5 水力压裂缝的平均半长(m) 150 
图 10 考虑致密储层流动特征情况下的采油速度及累计采油量对比a.产油速度对比；b.累计采油量对比Fig. 10. Comparison of oil recovery rate and cumulative oil production under tight reservoir flow characteristics通过图 9b与9c的对比可知，模拟到了第50 d时，水平井周围裂缝介质内的含水饱和度明显不同.不考虑致密储层边界层存在的储层，因其渗吸能力较强，在裂缝介质内吸走了更多的水，排出了更多的油，导致水平井周围裂缝介质内的含水饱和度明显低于考虑边界层的致密储层.
另外，从图 10a的产油速度对比图中可以看出，不考虑致密储层边界层特征的油藏其产油速度较高；同时，从图 10b的累积采油量对比图中可以更明显地看出，不考虑边界层的致密储层，其最终的累积采油量比考虑致密储层流动特征储层的最终累积采油量要高出一倍还多，说明如果在渗吸计算中忽视致密储层的边界层特征会严重高估渗吸对致密储层产能的影响.
5. 结论
(1) 本文通过考虑致密储层的边界层特征，建立了关于致密储层的渗透率、相渗曲线流动参数的解析表达式.
(2) 致密储层考虑边界层后，会大幅度减小储层的渗透率，从而减小储层的渗吸能力.
(3) 致密储层考虑边界层后，其毛管压力在原先的平缓段稍大于不考虑边界层的情况，同时边界层的存在使毛管力曲线原先的平缓段提前上翘，末端上翘段的毛管力明显大于未考虑边界层存在的情况.
(4) 致密储层考虑边界层后，缩小了油水原有流动通道的半径，使得原本迂回联通的流动通道更容易被截断形成不连续的流动通道，从而减小油相以及水相的流动能力，在相对渗透率曲线上则表现为油相及水相的相对渗透率因为边界层的存在而明显减小.
(5) 致密储层渗吸的强弱，主要由岩心的绝对渗透率决定；渗吸的作用范围，主要由相渗曲线的油水同流区决定；对于致密储层，虽然边界层的存在增大了毛管力，但同时也减小了储层的渗透率，缩小了有效渗吸的范围；综合来讲，边界层的存在一般会减弱致密储层的渗吸能力，从实例计算结果来看，考虑了边界层的致密储层其最终采收率仅为不考虑边界层致密储层的1/2.
(6) 从油藏尺度来看，渗吸作用对开采井近井地带的含油饱和度具有明显影响，渗吸能力越强近井地带的饱和度变化也就越强.对于致密油藏，因为边界层的存在使得致密储层渗吸能力大幅下降，导致在真实致密储层中因渗吸引起的饱和度分布上的变化是十分有限的；同时，从实例计算的结果来看，考虑边界层后累计采油量约为不考虑边界层情况的3/8，因此如果在渗吸计算中忽视致密储层的边界层特征会严重高估渗吸对致密储层产能的影响.
-
图 1 存在复杂缝网的逆向渗吸
Fig. 1. Counter-current imbibition in the presence of fracture network
图 2 边界层对致密储层流动的影响
a.存在边界层时流体在喉道中的流动；b.不存在边界层时流体在喉道中的流动
Fig. 2. The influence of boundary layer on the flow of tight throat
图 3 致密岩心OEO渗吸实验示意
Fig. 3. Schematic of tight core imbibition experiments with OEO boundary conditions
图 7 模型计算结果与实验数据对比
Fig. 7. Comparison between the calculated results and experimental data
图 9 不同情况下饱和度分布
a.初始时刻含水饱和度；b.50 d时不考虑边界层的含水饱和度；c.50 d时考虑边界层的含水饱和度
Fig. 9. Saturation distribution under the condition of without considering flow characteristics of tight oil reservoir
图 10 考虑致密储层流动特征情况下的采油速度及累计采油量对比
a.产油速度对比；b.累计采油量对比
Fig. 10. Comparison of oil recovery rate and cumulative oil production under tight reservoir flow characteristics
表 1 岩心基本参数
Table 1. Basic parameters of imbibition experiment under different boundary conditions
岩心编号 岩心直径(mm) 岩心长度(mm) 孔隙度(%) 束缚水饱和度(%) 模拟油密度(kg/m3) 地层水密度(kg/m3) S102 24.96 52.6 9.3 41.1 807 1 058 
下载: 导出CSV
表 2 油藏基本参数
Table 2. Basic parameters of reservoirs
参数名称 数值 初始油藏压力(10-1 MPa) 250 井底压力(10-1 MPa) 5 致密油藏基质渗透率(考虑边界层)(mD) 1.28 致密油藏基质渗透率(不考虑边界层)(mD) 7.56 油藏天然裂缝渗透率(mD) 100 水力压裂缝的渗透率(mD) 2 000 基质的孔隙度 0.093 裂缝介质的空隙度 0.001 水力压裂缝的孔隙度 0.5 水力压裂缝的平均半长(m) 150
下载: 导出CSV
-
[1] Aronofsky, J.S., Masse, L., Natanson, S.G., 1958.A Model for the Mechanism of Oil Recovery from the Porous Matrix due to Water Invasion in Fractured Reservoirs.Society of Petroleum Engineers, 213:17-19. http://www.onepetro.org/general/SPE-932-G [2] Bagherinezhad, A., Pishvaie, M.R., 2014.A New Approach to Counter-Current Spontaneous Imbibition Simulation Using Green Element Method.Journal of Petroleum Science and Engineering, 119:163-168.doi: 10.1016/j.petrol.2014.05.004 [3] Cai, J.C., Yu, B.M., 2012.Advances in Studies of Spontaneous Imbibition in Porous Media.Advances in Mechanics, 42(6):735-754 (in Chinese with English abstract). http://en.cnki.com.cn/Article_en/CJFDTOTAL-LXJZ201206006.htm [4] Cai, J.C., Guo, S.L., You, L.J., et al., 2013.Fractal Analysis of Spontaneous Imbibition Mechanism in Fractured-Porous Dual Media Reservoir.Acta Physica Sinica, 62(1):220-224 (in Chinese with English abstract). http://en.cnki.com.cn/Article_en/CJFDTotal-WLXB201301035.htm [5] Fang, S.D., Cheng, L.S., Ayala, L.F., 2017.A Coupled Boundary Element and Finite Element Method for the Analysis of Flow through Fractured Porous Media.Journal of Petroleum Science and Engineering, 152:375-390.doi: 10.1016/j.petrol.2017.02.020 [6] Hassan, S., 2004.Analysis Scaling and Simulation of Counter-Current Imbibition (Dissertation).London Imperial College, London. [7] Jia, P., Cheng, L.S., Huang, S.J., et al., 2015.Transient Behavior of Complex Fracture Networks.Journal of Petroleum Science and Engineering, 132:1-17.doi: 10.1016/j.petrol.2015.04.041 [8] Kashchiev, D., Firoozabadi, A., 2003.Analytical Solutions for 1D Countercurrent Imbibition in Water-Wet Media.SPE Journal, 8(4):401-408.doi: 10.2118/87333-pa [9] Li, Y., Lei, Q., Liu, X.G., 2011.Characteristics of Micro Scale Nonlinear Filtration.Petroleum Exploration & Development, 38(3):336-340 (in Chinese with English abstract). http://en.cnki.com.cn/Article_en/CJFDTOTAL-SKYK201103018.htm [10] Ma, S.X., Morrow, N.R., Zhang, X.Y., 1997.Generalized Scaling of Spontaneous Imbibition Data for Strongly Water-Wet Systems.Journal of Petroleum Science and Engineering, 18(3-4):165-178.doi: 10.1016/s0920-4105(97)00020-x [11] Meng, Q., Liu, H.Q., Wang, J., 2017.A Critical Review on Fundamental Mechanisms of Spontaneous Imbibition and the Impact of Boundary Condition, Fluid Viscosity and Wettability.Advances in Geo-Energy Research, 1(1):1-17. [12] Ren, D.Z., Sun, W., Huang, H., et al., 2016.Formation Mechanism of Chang 6 Tight Sandstone Reservoir in Jiyuan Oilfield, Ordos Basin.Earth Science, 41(10):1735-1744 (in Chinese with English abstract). http://www.en.cnki.com.cn/Article_en/CJFDTotal-DQKX201610009.htm [13] Schechter, D.S., Zhou, D., Orr, F.M., 1994.Low IFT Drainage and Imbibition.Journal of Petroleum Science and Engineering, 11(4):283-300.doi: 10.1016/0920-4105(94)90047-7 [14] Tian, X.F., 2015.Non-Linear Flow Behavior and Application in Ultra-Low Permeability Reservoirs (Dissertation).China University of Petroleum, Beijing (in Chinese with English abstract). [15] Tian, X.F., Cheng, L.S., Cao, R.Y., et al., 2016.Characteristics of Boundary Layer in Micro and Nano Throats of Tight Sandstone Oil Reservoirs.Chinese Journal of Computational Physics, 33(6):717-725 (in Chinese with English abstract). [16] Wei, Q., Li, Z.P., Bai, R.T., et al., 2016.An Experimental Study on the Effect of Microscopic Pore Structure on Spontaneous Imbibition in Tight Sandstones.Petroleum Drilling Techniques, 44(5):109-116 (in Chinese with English abstract). [17] Wu, S.T., Zou, C.N., Zhu, R.K., et al., 2015.Reservoir Quality Characterization of Upper Triassic Chang 7 Shale in Ordos Basin.Earth Science, 39(11):1810-1823 (in Chinese with English abstract). http://www.en.cnki.com.cn/Article_en/CJFDTotal-DQKX201511004.htm [18] Xu, Z.Y., Cheng, L.S., Cao, R.Y., et al., 2017.Simulation of Counter-Current Imbibition in Single Matrix and Field Scale Using Radical Integral Boundary Element Method.Journal of Petroleum Science and Engineering, 156:125-133.doi: 10.1016/j.petrol.2017.05.016 [19] Zou, C.N., Zhu, R.K., Wu, S.T., et al., 2012.Types, Characteristics, Genesis and Prospects of Conventional and Unconventional Hydrocarbon Accumulations: Taking Tight Oil and Tight Gas in China as an Instance.Acta Petrolei Sinica, 33(2):173-187 (in Chinese with English abstract). doi: 10.1038/aps.2011.203 [20] 蔡建超, 郭士礼, 游利军, 胡祥云, 2013.裂缝-孔隙型双重介质油藏渗吸机理的分形分析.物理学报, 62(1): 220-224. http://www.cnki.com.cn/Article/CJFDTOTAL-WLXB201301035.htm [21] 蔡建超, 郁伯铭, 2012.多孔介质自发渗吸研究进展.力学进展, 42(6): 735-754. doi: 10.6052/1000-0992-11-096 [22] 李洋, 雷群, 刘先贵, 2011.微尺度下的非线性渗流特征.石油勘探与开发, 38(3): 336-340. http://www.cnki.com.cn/Article/CJFDTOTAL-SKYK201103018.htm [23] 任大忠, 孙卫, 黄海, 等, 2016.鄂尔多斯盆地姬塬油田长6致密砂岩储层成因机理.地球科学, 41(10): 1735-1744. http://earth-science.net/WebPage/Article.aspx?id=3375 [24] 田虓丰, 2015. 超低渗透油藏非线性渗流特征与应用(博士学位论文). 北京: 中国石油大学. [25] 田虓丰, 程林松, 曹仁义, 等, 2016.致密油藏微纳米喉道中的边界层特征.计算物理, 33(6): 717-725. http://www.cnki.com.cn/Article/CJFDTOTAL-JSWL201606011.htm [26] 韦青, 李治平, 白瑞婷, 等, 2016.微观孔隙结构对致密砂岩渗吸影响的试验研究.石油钻探技术, 5: 109-116. http://www.cnki.com.cn/Article/CJFDTOTAL-SYZT201605024.htm [27] 吴松涛, 邹才能, 朱如凯, 等, 2015.鄂尔多斯盆地上三叠统长7段泥页岩储集性能.地球科学, 39(11): 1810-1823. http://earth-science.net/WebPage/Article.aspx?id=3188 [28] 邹才能, 朱如凯, 吴松涛, 等, 2012.常规与非常规油气聚集类型、特征、机理及展望——以中国致密油和致密气为例.石油学报, 33(2): 173-187. doi: 10.7623/syxb201202001 期刊类型引用(9)
1. 王强，赵金洲，胡永全，赵超能，张祯祥. 岩心尺度静态自发渗吸的数值模拟. 石油学报. 2022(06): 860-870 . 
百度学术2. 刘凯，赵洋，王维波，杨红. 界面张力与润湿性对化学渗吸的影响——以延长油田杏子川低渗油藏为例. 石油地质与工程. 2022(04): 123-126 . 百度学术3. 邸士莹，程时清，白文鹏，魏操，汪洋，秦佳正. 致密油藏动态裂缝扩展机理及应用. 力学学报. 2021(08): 2141-2155 . 百度学术4. 曹仁义，程林松，杜旭林，时俊杰，杨晨旭. 致密油藏渗流规律及数学模型研究进展. 西南石油大学学报(自然科学版). 2021(05): 113-136 . 百度学术5. 刘凯，高振东，王成俊，高怡文，孟选刚. 低渗透油藏渗吸采油技术研究进展. 油田化学. 2021(04): 747-753 . 百度学术6. 徐中一，方思冬. 边界层对致密油藏逆向渗吸的数值模拟. 石油与天然气地质. 2020(03): 638-646 . 百度学术7. 阮迪，屈亚光，沈畅，刘德华. 裂缝性致密砂岩油藏渗吸驱油效果影响因素. 科学技术与工程. 2019(29): 118-123 . 百度学术8. 唐洪明，朱柏宇，王茜，张烈辉，赵峰，王俊杰. 致密砂岩气层水锁机理及控制因素研究. 中国科学:技术科学. 2018(05): 537-547 . 百度学术9. 苏玉亮，韩秀虹，王文东，盛广龙，杨正明，何英，王志远. 致密油体积压裂耦合渗吸产能预测模型. 深圳大学学报(理工版). 2018(04): 345-352 . 百度学术其他类型引用(5)
-
下载:
百度学术
点击查看大图
计量
- 文章访问数: 4186
- HTML全文浏览量: 1850
- PDF下载量: 13
- 被引次数: 14
