ACTA Scientiarum Naturalium Universitatis Pekinensis
Dimensionless Model and Numerical Solution for 1D Horizontal Spontaneous Imbibition
ZHAO Guozhong
Exploration & Development Institute of Daqing Oilfield Ltd, Daqing 163712; E-mail: zhaoguozh@petrochina.com.cn
Abstract The reservoir parameters affecting imbibition performance are normalized. The dimensionless model for 1D horizontal spontaneous imbibition in consideration of capillary back pressure is derived. In the model, the normalized water saturation and dimensionless water pressure are chosen as the primary unknown variables. The capillary pressure function only appears in one PDE in the model which is convenient to be solved numerically. The numerical solution can be obtained by the fully implicit and conservative difference scheme with the saturation function value between two adjacent points to be set on the value at the arithmetic average saturation for the two points. The solution of this model with relative permeability and capillary pressure curves matching power law is discussed. The assumption used by prior researchers to obtain the analytical solution for 1D incompressible spontaneous imbibition problem is pointed out to be only partially applicable for purely countercurrent case. More numerical results for co-current imbibition cases are also obtained and discussed, and some preliminary cognition about the variation of imbibition performance along with reservoir parameters is concluded. Key words spontaneous imbibition; dimensionless model; numerical solution; capillary back pressure
1 自发渗吸模型介绍
渗吸是裂缝型油藏水驱开采的关键机理之一。近年来, 石油行业积极尝试用人工压裂造缝方法开采特低、超特低渗透油藏。此过程中, 外来液对油层的自发渗吸作用对开采效果有重要的影响。因
此, 自发渗吸问题逐渐得到较多实验[1–4]、理论[5–7]和计算[8–10]研究的关注。渗吸的机理是, 两相流体通过多孔介质流动时, 相间压力差(即毛管压力)及其随相饱和度的变化对流体流动起控制作用。自发渗吸分为逆向渗吸和同向渗吸两种方式, 这两种自发渗吸方式的产生条件、采油能力和时效性是研究
渗吸采油模式、预测产能动态、优化开发设计的基础。有关两相流体通过多孔介质流动的一般数学模型及其数值求解方法是可以考虑毛管压力作用的,比如ECLIPSE。但是, 为了研究纯自发渗吸作用对开采特征的影响, 还需要对一般模型进行简化和适当处理(如考虑特殊的边界条件、彻底忽略压缩性等), 才能准确地反映自发渗吸的物理过程, 揭示渗吸特征的影响因素。考虑到均匀多孔介质内不可压、不混溶油水两相流体的一维水平流动, 则达西速度、连续性方程及约束条件为
Kk Pw
w rw , (1) w x Kk P ro o , (2) o o x
Sw w
0, (3) t x So o 0, (4) t x P P P (S ), (5) o w c w Sw So 1, (6)
,, K, 其中 和 分别表示介质的绝对渗透率和孔隙度; vk , P和 S 分别表示相达西速度、相相对r渗透率、相黏度、相压力和相饱和度; 下角标 w 和o 分别对应水相和油相; P 为毛管压力。c将式(1)和(2)分别代入式(3)和(4), 再与式(5)和(6)联立, 消去变量 Pw 和 So 后, 在自发渗吸的初始边界条件下, 便构成以P ( x , t) 和 S ( x , t) 为未知变o w量的偏微分方程定解问题。该问题是非线性的, 一[8]般得不到解析解。Pooladi-darvish等 曾对自发渗吸问题进行数值求解。在特殊的边界条件假设下,通过叠代求解一个关于饱和度函数的积分方程便能够得到解析解。为此, 设总流速为 , 由t w o式(1), (2)和(5)得到
将式(7)代入式(3), 并考虑不可压条件, 可以得到
S v d f S S w t w w D w 0, (9) t d Sw x x x 其中D定义为
K dp D o f kro c 。 (10) w dsw式(9)为一维不可压渗吸问题的饱和度方程。 0 t时为纯逆向渗吸, 0时为同向、逆向渗吸并t存。纯同向渗吸只能在特殊的实验条件下才有可能发生, 在油藏条件下出现的可能性很小。Mcwhorter等[5]基于对某些渗吸实验结果中产油率数据的观察, 从以下假设条件出发:
( x , t ) F ( x , t ) (0, t ) F ( x , t ) At 1 2 , (11) w w其中 F为未知待定的关于位置和时间的函数, 通过叠代求解一个关于F ( S ) 的积分方程, 就能够得到
w纯逆向或纯同向渗吸的解析解。但是, 对于纯同向渗吸的情况, 只适合[0, ) 的问题域, 且积分方程
[6]的适定性不能保证。Chen 等 针对纯逆向渗吸情况的研究得到相同的结果。纯逆向渗吸的解在文献[8]中也得到验证, 表明此解的可靠性较高, 也说明假设条件(式(11))对纯逆向渗吸有一定的合理性。
[7] Nooruddin 等 给出同向、逆向渗吸并发情况的解析解, 但仍然要从式(11)的假设条件出发, 并通过求解积分方程才能得到。
以上工作存在两方面的局限性: 1) 式(11)的假设条件是观察有限个纯逆向渗吸实验结果得到的近似计算公式, 是否对所有纯逆向渗吸都适用? 对逆向、同向渗吸并发的情况是否适用? 目前没有可靠的论证; 2) 入口端的油相压力边界条件等价于毛管压力为零, 这在物理学上不具一般性, 无法解释自发渗吸实验[4]中发现的毛管回压现象, 所以基于该假设条件的解析解可能存在一定的偏差。因此, 数值求解方法应是研究自发渗吸问题继续采用的手段之一。
对于较成熟的油藏模拟求解技术, 数值求解一维自发渗吸问题的方法从表面上看十分简单, 但目前油藏模拟器(尤其是商业化的油藏模拟器)大多不能模拟不可压问题, 且只能处理第二类边界条件(块中心网格系统), 那么是否可以利用已有油藏模拟器来研究自发渗吸问题？Schembre等[9]和li等[10]
在有关渗吸问题的研究中使用了商业油藏模拟器ECLIPSE, 他们都在问题域的入口边界网格外附加了一个虚拟网格。这样做是利用已有模拟器数值求解自发渗吸问题的必由之路, 可以得到该问题大致的解。但是, 在这样直接模拟计算中, 边界条件不能严格满足, 且岩石和流体的压缩性必然对收敛性及得到的解产生影响。在使用商业油藏模拟器的过程中, 压缩系数不可无限小, 一般要大于10–15 kpa–1才能正常运行, 更不允许出现零值。但是, 文献[9– 10]均未给出压缩系数的取值, 其影响的大小是未知的。所以, 将商业油藏模拟器得到的自发渗吸模型的解作为标准是不合适的。本文从与渗吸相关油藏参数的归整及模型的无因次化入手, 给出考虑毛管回压的一维自发渗吸问题无因次模型及其数值求解方案, 通过对数值解的讨论, 指出前人在求解析解时所用假设条件的局限性, 并分析渗吸参数对开采特征的影响。
2 渗吸相关油藏参数的归整
从式(9)可以看出, 问题的非线性是由两个关于含水饱和度的函数 f 和D 引起的。其中, f 由两
w w相的相对渗透率曲线和流度比确定; D 除与 f 具
w有同样的影响因素外, 还取决于毛管压力曲线。这些因素确定后, 在一定的初始边界条件下, 解是确定的。无论是用解析方法先求式(9)的解, 还是用数值方法直接求原方程组的解, 对描述这些因素的油藏参数进行综合归整, 对模型的无因次化及求解都十分有帮助。
为此, 设两相的相对渗透率曲线以及油相对水相的流度比分别为
k k k w( S n) ( S [0, 1], kw [0, 1]) , (12) rw rwm n k k k o( S n) ( S [0, 1], ko [0, 1]) , (13) ro rom n krom o
, (14) krwm w
其中, krwm 和 krom 分别为水相、油相最大相对渗透率, k 和 ko 分别为水相、油相归一化相对渗透率, w Sn 为归一化水相饱和度: Sw Smn Sn , (15) S Smn mx Smx 和 Smn分别为水相最大、最小饱和度。 相间毛管压力差及其随饱和度的变化对流体的流动方向和流速的影响是产生渗吸现象的关键作用机理。吸入毛管压力曲线对渗吸问题的求解至关重要, 它与岩石孔隙分布、两相流体的性质以及两相流体岩石相互作用性质密切相关。前人通过归纳实验资料, 得到一些关于毛管压力与湿相归一化饱和度的经验公式, 其中经常使用的有简单幂函数型(Ⅰ型)[7]: (1 Sn ) ( 0 ); 简单幂函数型(Ⅱ型)[11]:
S 1 ( 1 ); 复合幂函数型( Ⅲ型)[12]: (S m 1 1)1m n n (0 m 1 ); 对数型( Ⅳ型)[8]: ln Sn ; 等等。除Ⅱ型外, 其他公式在归一化饱和度为1处都为0。然而,广泛运用的压汞法测得的退汞毛管压力曲线显示,吸入曲线末点毛管压力值一般都大于零。自发渗吸实验也发现, 渗吸入口端存在毛管回压现象[2 4]。此外, Ⅱ, Ⅲ和Ⅳ型在归一化饱和度Sn = 0处是无界的。由此可知, 直接使用这些公式都有局限性, 也不便进行参数化分析。为了解决端点无界和毛管回压描述问题, 本文从几个可调的特征量出发, 通过端点微移法, 给出可以近似地反映上述经验关系的统一表达方式, 以便分析自发渗吸问题解的特征与吸入毛管压力曲线特征的关系, 从而认识自发渗吸规律, 指导渗吸实验的数值模拟拟合及特征参数的反求。假设吸入毛管压力曲线的起始压力为P , 末
cm点压力为 P (0, P ), 二者之比 1。将各型的
ce cm关系归纳为如下的一般形式: P ( S ) P ( S )( S [0, 1], [1, ]) , (16)
cm c n n n
其中 ( Sn )来自上述任意一种类型的表达式。例如,对于Ⅰ型, 取下式可直接满足要求:
( S n) ( 1)(1 S n) 1; (17)
对于Ⅱ型, 可先设 ( S n) Sn 1 ;对于Ⅲ型, 先设 (S ) (S m 1 1)1 m 1; n n对于Ⅳ型, 先设
( Sn ) ln Sn +1。然后, 用端点微移法找到符合式(16)的 (Sn)。为
此, 考虑待定的小量 0, 做以下定义:
( S n) ((1 ) S )( Sn [0, 1] ) , (18) n可以看出, 若 满足以下关系: 1( ), (19)则 ( Sn )可满足式(16)的要求, 于是式(18)定义的函数即为寻找目标。对于Ⅱ, Ⅲ和Ⅳ型经验关系(不限于这3 种类型), 式(19)是可解的。统一归纳为式(16)形式表达式的好处是, 只要给出 P , P 和 cm ce中的任意两个特征量, 吸入毛管压力曲线就能得到描述。这样, 不仅可以解决上述问题, 还能灵活地选择 ( Sn ), 各自的内部特征量(如 , , m等)的变化对毛管压力的影响也可以进行描述。另外, 对于Ⅰ , Ⅲ和Ⅳ型经验关系, 由于原型无毛管回压, 所以可直接构造如下形式的毛管压力曲线: P ( S ) P ( S )( S [0, 1], [0, 1]) 。 (20) c n cm n n对于Ⅰ型, 下式可满足要求: ( S n) (1 Sn ) ; (21)对于Ⅲ型, 直接设
( S ) ( S m 1 1)1m ; n n对于Ⅳ型, 直接设 ( S n) ln S 。n对于指定的小量 0, 做以下定义可符合式(20)的要求:
1 ( S n) ((1 ) S ) ( Sn [0, 1] )。 (22) ( ) n
下文将式(16)和(20)统一归纳为式(16)的形式,式(20)在形式上相当于式(16)中的 不出现, 并将 ( Sn )的值域变为[0, 1]。3 渗吸模型的无因次化将 ( Sn )归一化为 ( Sn )
( S ) , n n max ( Sn )
则由式(16)可统一求得:
将式(27)和(30)联立, 便是以水相无量纲压力和归一化饱和度为未知变量的油水两相不可压一维均质水平渗流问题的无因次模型。该模型只在一个偏微分方程(27)中含有毛管压力项, 便于处理。如图 1 所示, 设岩样初始孔隙压力P 0, qo cou i和qococ分别表示自发逆向和同向渗吸速率, 则自发渗吸无因次化模型的初始条件为
s ( x ,0) 0, p ( x ,0) 。 (31) max 纯逆向自发渗吸(图 1(a))时, 左端与水相接触,右端封闭。考虑左端水相的连续性, 应满足以下边界条件:
p (0, t ) 0, s (0, t ) 1, (32) p s (1, t ) 0, (1, t ) 0。 (33) x x应当指出, 左端油相可以流出, 但考虑毛管回压时, 压力并不连续, 因此选择水相压力作为主未知变量, 左端油相压力边界条件已隐含在水相饱和度边界条件中。同向、逆向自发渗吸并存(图1(b))时, 出口无水相流出, 且出口油相与入口水相水力压头相等,油相压力(即水相压力与毛管压力之和)保持初始值不变, 所以除满足式(32)外, 还应满足以下边界条件: ( s (1, t )) v (1, t ) 0, p (1, t ) 0 , (34) w max 对式(27)和(30)采用全隐式守恒型差分格式, 两邻点间的饱和度函数取其饱和度算数平均值对应的值(数值实验表明, 按传统模拟器采用流度上游权格式对本文讨论的一维不可压无量纲渗吸问题不收敛), 可构造自发渗吸模型的差分方程组(见附录),然后用牛顿型非线性代数方程组数值解法叠代求解, 在每个牛顿叠代步中还要求解一个块三对角线性系统, 块的阶数为2。 4 数值算例及结果分析可以看出, 给定 , , kw(s), ko(s)和 (s)后, 无 因次模型(式(27)和(30))的解就是确定的。数值求解该模型的好处是, 一旦得到对应的数值解, 其他油藏参数的变化对渗吸特征的影响可以通过对无量纲时间与真实时间之间的转换直接得到, 无须再一一求解原问题。
本文以Intel FORTRAN和PETSC (http://www. mcs.anl.gov/petsc, 其KSP构件可用于求解块三对角线性系统)为开发工具, 实现对无因次模型的数值求解。相对渗透率曲线采用的幂律型为
其中 nw 和 no 分别为水、油相归一化相对渗透率幂指数。所用毛管压力曲线是式(17)和(21), 分别对应考虑毛管回压和不考虑毛管回压。
本文算例取网格数 N 101 , 无量纲时间步长 t 9 105 , 计算过程中压力方程和饱和度方程的最大误差均小于9 107 , 终止条件为水饱和度前缘到达右端且可采储量采出程度高于50%。将 nw 2.5, no=2.0, =3.0, =6 和=0.4 的算例作为基本算例, 分别对纯逆向、同逆并发两种渗吸方式以及是否考虑毛管回压两种情况进行 4次数值求解。图2 给出同逆并发渗吸、考虑毛管回压时的水相 归一化饱和度和无量纲压力分布。图中横坐标 xd 为到入口端的无量纲距离, 纵坐标 tdsq 为归一化无量纲时间的平方根。
为了分析模型参数对渗吸特征的影响, 针对同逆并发且带毛管回压的情况, 我们变更 , , n o, nw , 等参数, 进行一系列算例的数值求解。
4.1 渗吸速率的计算结果
图 3给出计算得到的无量纲渗吸速率与无量纲
时间的关系。在纯逆向渗吸情况下, 考虑与不考虑毛管回压的无量纲渗吸速率结果相同。这是因为两种情况的差分方程完全相同, 初始水相压力只相差一个因子 /( 1) , 不影响饱和度的计算结果。但是, 由于两种情况的无量纲时间也相差这一个因子,所以换算到原模型后, 结果会有一定的差异。从图3 可以看出, 同逆并发渗吸的早期(在无量纲时间0.1 之前), 渗吸速率与纯逆向渗吸基本上相同, 说明主要是逆向渗吸在起作用; 此后二者明显偏离,说明同向渗吸开始起主导作用, 使得同逆并发渗吸的水相饱和度前缘推进更快。为了分析渗吸速率与 时间的关系, 在图 3 中画出一条斜率为–0.5 的参考虚直线。
可以看到, 同逆并发情况下不存在斜率为–0.5的虚直线与任意渗吸速率线匹配。对于纯逆向渗吸, 在无量纲时间 0.1 之后至水相饱和度前缘到达右端之前, 存在一条斜率为–0.5 的虚直线与渗吸速率线近似较好, 但在渗吸的早期, 尤其在无量纲时间 0.001 之前, 也不存在这样的近似直线。这表明,式(11)的假设条件仅仅对于纯逆向渗吸的无限作用流动期的后段有较高的可靠度。对于同逆并发渗吸, 该假设的偏差太大, 由此得到的解析解可信度有限。另外, 在渗吸的早期, 基于该假设条件得到的渗吸速率一定偏高。由于在采油实践中, 渗吸的
[4]早期特征意义更大, 所以受到较多关注。Ferno等指出渗吸实验得到的速率结果明显低于理论值, 因此本文得到的数值解可能更可信。
4.2 毛管回压的存在对同逆并发渗吸采出程度及其构成的影响
图 4给出考虑和不考虑毛管回压两种情况下同逆并发渗吸可采储量采出程度及其构成的对比。考虑毛管回压时, 同逆并发渗吸的总采出程度及其同向采出部分均高于不考虑毛管回压的对应值, 总采出程度中的逆向采出部分则相反。这是考虑毛管回压得到的新认识。另外, 同逆并发渗吸的总采出程度中的同向采出量占主导地位, 逆向采出量在大部分时段几乎不增长, 但考虑毛管回压后逆向采出的油能超过可采储量的10%。这一认识与以往的研究也略有差别。
4.3 渗吸参数对渗吸特征量的影响
水相饱和度前缘到达右端的无量纲时间、该时间可采储量采出程度和采出50%可采储量所用无量纲时间是自发渗吸的重要特征量, 分别用Tfr , fr 和Tr50 表示。这3个特征量对渗吸采油的开发决策和方案设计至关重要。由无量纲时间的定义(式(29))可知, 油藏绝对渗透率、孔隙度、水相最大相对渗透率、水相黏度、水相饱和度范围、最大毛管压力、岩石沿渗吸方向的长度等油藏参数并不直接出现在无因次模型中, 它们对自发渗吸特征的影响程度可以在得到数值解后通过简易的变换得到, 这里不做讨论。
计算结果显示, 在基本算例条件下, 纯逆向渗吸 Tfr 为 16.575, fr 为 0.340, Tr50 为 37.000, 在同逆并发并且考虑毛管回压的情况下, 这 3 个值分别为8.5127, 0.456 和 9.6415, 表明同逆并发渗吸的采油速度远远大于纯逆向渗吸。在渗吸采油实践中, 应尽量将渗吸方式从纯逆向转化为同逆并发, 达到在有限的时间内获得更高采收率的目的。
图5显示最大毛管压力与回压之比 的变化对同逆并发渗吸特征量的影响。这里已将无量纲时间换算为与同逆并发渗吸、无毛管回压相同的尺度。可以看出, 随 的增加, Tfr 和 Tr50 都增长, 而fr 下降。因此, 的增大不利于渗吸采油。但是, 当 >10 后, 影响程度明显变弱。事实上, 当时, 同逆并发渗吸特征量应趋于不考虑毛管回压情 况下的对应值, 分别为T =11.347, fr = 0.401 和
fr Tr50 =16.604。另外, 当 >4时, 水相前缘到达右端之后, 才能采出50%的可采储量。依据计算结果, 可得出毛管压力曲线幂律指数 和油相相对渗透率曲线幂律指数no 的变化对同逆并发渗吸特征量的影响为: 和 no 的增大也不利于渗吸采油的认识。
图6显示油相与水相流度之比 的变化对同逆并发渗吸特征量的影响。可以看出, 随 增加, Tfr和 Tr50 都减小, 而fr 增加。因此, 的增大有利于
渗吸采油。但是, 当>0.5 时, 影响程度明显变弱。另外, 当 >0.7 时, 水相前缘到达右端之前, 就可采出 50%的可采储量。
图 7显示水相相对渗透率曲线幂律指数nw 的变化对同逆并发渗吸特征量的影响。随nw 增加, 3个特征量都增加。如果更关注采油时间, nw的增大不利于渗吸采油; 如果更关注采出程度, nw 的增大有利于渗吸采油。值得重视的是, 当 nw >3时, 水相前缘到达右端之前就可采出50%的可采储量, 且随nw 增大, Tr50 和 Tfr 的差距变大, 所以 nw的增大有利于渗吸采油。另外, 有限时间内采出更多的油, 说明水相饱和度前缘更陡, 其分布可能接近水油段塞流。这一点值得通过调整水相相对渗透率曲线拟合自发渗析实验结果时参考。
5 主要认识和结论
本文考虑毛细管回压, 建立了基于 4 种类型毛管压力与湿相归一化饱和度经验关系的, 两相不可压、不混溶流体通过一维均匀多孔介质的自发渗吸流动统一无因次模型。该模型选择水相归一化饱和度和无量纲压力作为主未知变量, 采用全隐式守恒型差分格式, 对两邻点间相对渗透率和毛管压力取平均饱和度下的值可以得到该问题的数值解。依据幂律型相对渗透率和毛管压力曲线情形下本文模型的数值解, 得出前人在求取不考虑毛管回
压的一维不可压自发渗吸问题解析解过程中, 所用的关键性假设条件仅对纯逆向渗吸的无限作用流动期的后段有一定的可靠度。对于同逆并发渗吸, 该假设的偏差太大, 从而得到的解析解不可信。另外,在渗吸的早期, 前人所用假设条件会导致渗吸速率偏高。毛管回压的存在对同逆并发渗吸采出程度及其构成有明显影响。考虑毛管回压时, 同逆并发渗吸的总采出程度及其同向采出部分均高于不考虑毛管回压的对应值; 总采出程度中的逆向采出部分则相反。考虑毛管回压后, 同逆并发渗吸的逆向采出量可能超过可采储量的10%。最大毛管压力与毛管回压之比、毛管压力曲线幂律指数和油相相对渗透率曲线幂律指数这 3 个油藏参数的增大都会减缓水相推进速度, 且使采出程度下降, 不利于渗吸采油, 本文定义的油相与水相流度之比的增大有利于渗吸采油。当水相相对渗透率曲线幂律指数大于 3 时, 水相前缘到达右端之前,就可采出 50%的可采储量, 有利于渗吸采油。