# Dimensionless Model and Numerical Solution for 1D Horizontal Spontaneous Imbibition

## ZHAO Guozhong

ACTA Scientiarum Naturalium Universitatis Pekinensis - - Contents -

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 自发渗吸模型介绍

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)得到

 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)的假设条件出发, 并通过求解积分方程才能得到。

2 渗吸相关油藏参数的归整

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

S   1 (   1 ); 复合幂函数型( Ⅲ型)[12]: (S m 1  1)1m 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)

 ( Sn )   ln Sn +1。然后, 用端点微移法找到符合式(16)的  (Sn)。为

 ( 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)1m ; n n对于Ⅳ型, 直接设 ( S n)  ln S 。n对于指定的小量  0, 做以下定义可符合式(20)的要求:

1  ( S n)   ((1  ) S ) ( Sn [0, 1] )。 (22)  ( ) n

  ( S ) , n n max   ( Sn )

 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))的解就是确定的。数值求解该模型的好处是, 一旦得到对应的数值解, 其他油藏参数的变化对渗吸特征的影响可以通过对无量纲时间与真实时间之间的转换直接得到, 无须再一一求解原问题。

4.1 渗吸速率的计算结果

[4]早期特征意义更大, 所以受到较多关注。Ferno等指出渗吸实验得到的速率结果明显低于理论值, 因此本文得到的数值解可能更可信。

4.2 毛管回压的存在对同逆并发渗吸采出程度及其构成的影响

4.3 渗吸参数对渗吸特征量的影响

fr Tr50 =16.604。另外, 当 >4时, 水相前缘到达右端之后, 才能采出50%的可采储量。依据计算结果, 可得出毛管压力曲线幂律指数 和油相相对渗透率曲线幂律指数no 的变化对同逆并发渗吸特征量的影响为:  和 no 的增大也不利于渗吸采油的认识。

5 主要认识和结论