稠油开采中非等温压强扩散方程的数学建模
Mathematical Modeling of Non-Isothermal Pressure Diffusion in Heavy Oil Recovery
摘要: 在常见的稠油开采技术建模中，往往忽略压强影响或者统一用等温压强扩散方程替代非等温压强扩散方程进行计算求解。本文针对蒸汽辅助重力泄油(SAGD)、温溶剂萃取采油(warm VAPEX)两种采油技术分别建立考虑温度影响的压强扩散方程和同时考虑温度和浓度影响的压强扩散方程。并以SAGD为例对比等温与非等温压强扩散方程计算得到的压强分布，发现随着温度升高，两者差异逐渐增大，说明在有温度作用的采油技术中，更应该建立非等温压强扩散模型。
Abstract: In common heavy oil recovery technology modeling, the impact of pressure is often neglected, or the isothermal pressure diffusion equation is used as a substitute for the non-isothermal pressure diffusion equation for calculation and solution. In this paper, pressure diffusion equations that account for the temperature effect and pressure diffusion equations that consider both temperature and concentration effects are established for two oil extraction techniques: Steam Assisted Gravity Drainage (SAGD) and warm VAPEX (Warm Solvent Extraction for Oil). Taking SAGD as an example, the pressure distributions obtained from the isothermal and non-isothermal pressure diffusion equations are compared. It is found that as the temperature increases, the differences between the two distributions gradually grow, indicating that in oil extraction techniques affected by temperature, a non-isothermal pressure diffusion model should be established.
文章引用：梁芳, 王乔, 张艳. 稠油开采中非等温压强扩散方程的数学建模[J]. 石油天然气学报, 2024, 46(4): 472-480. https://doi.org/10.12677/jogt.2024.464059

1. 引言

我国稠油资源丰富，但由于稠油粘度高密度大，不易流动，因此在开采中面临很多挑战[1]。关于稠油开采，目前提出的技术有蒸汽辅助重力泄油(SAGD)、温溶剂萃取采油(warm VAPEX)等，任何采油技术在具体油藏中应用之前，都要进行模拟和预测，其中数学建模作为快速有效的预测方法，在油藏工程中应用广泛。

现有文章关于SAGD或warm VAPEX建模时，往往忽略压强影响，但事实上由于压差的存在，会导致对流速度增大，若直接忽略压强影响，计算得到的泄油速率会比实际偏小。此外个别文章在考虑压强影响时，建立的均为恒温下的压强扩散模型[2] [3]，实际上无论是在SAGD还是warm VAPEX中均应建立非恒温压强扩散方程。另外由于warm VAPEX同时存在传热、传质机制，因此在建立压强扩散方程时除了考虑温度还应考虑到溶剂的影响。

本文创新地详细推导了SAGD中非等温压强扩散方程，并且推导了warm VAPEX对应的考虑温度和溶剂浓度变化的压强扩散方程。以SAGD为例，通过在同一条件下分别计算等温扩散模型与本文模型，发现两者的压强分布差异会随着温度升高而增大，说明在有温度作用的采油技术中，有必要建立非等温压强扩散模型。

2. 数学建模

图1为腔室的横截面图。本文忽略平行于井的压强变化，只关注垂直于边界层的压强扩散。图2为腔室和稠油接触的过渡区物理模型图，其中以垂直过渡区边界方向为x轴建立坐标系。

Figure 1. Cross section of SAGD/warm VAPEX

图1. SAGD/warm VAPEX横截面图

Figure 2. Physical model of transition zone

图2. 过渡带物理模型图

2.1. SAGD压强扩散推导

由质量守恒得到连续方程为：

∂ ( ρ o ϕ ) ∂ t + ∂ ( ρ o V ) ∂ x = 0 (1)

其中 ϕ 为孔隙度，无量纲，本文假设孔隙度为常数； ρ o 表示稠油密度， kg / m 3 ； V 表示对流速度，可以根据达西公式计算得到

V = − k k r o μ o ( ∂ p ∂ x ) (2)

其中 k 为绝对渗透率，达西； k r o 为油相相对渗透率，无量纲； p 为压强， Pa ； μ o 为稠油粘度， Pa ⋅ s ，并且粘度为温度函数：

ln ln ( μ o ) = − 3.6261 ln ( T + 273.15 ) + 22.8339 (3)

稠油饱和度和温度之间的关系采用Sharma和Gates的假设，而相对渗透率和沥青饱和度之间的相关性参考Corey的相关性[4]：

S o = S o r + ( S o i − S o r ) ( 1 − T − T r T s − T r ) (4)

k r o = k r o , r w ( S o − S o r 1 − S w c − S o r ) m (5)

其中 T 为温度，℃； S o i 和 S o r 分别是初始和残余油饱和度，无量纲； S w c 是原生水饱和度，无量纲； k r o , r w 是最小含水饱和度下的石油相对渗透率，无量纲； T r 和 T s 分别为油藏温度和注入温度，℃； m 是Corey系数，无量纲。

非等温条件下的状态方程：

ρ o = ρ o i [ 1 + c f ( p − p 0 ) − β ( T − T 0 ) ] (6)

因此：

ρ o ϕ = ϕ ρ o i [ 1 + c f ( p − p 0 ) − β ( T − T 0 ) ] (7)

其中 c f 为稠油压缩系数， Pa − 1 ； β 为稠油的热膨胀系数， ˚ C − 1 ；带入连续方程等号左侧第一项得：

∂ ( ρ o ϕ ) ∂ t = ϕ ρ o i ( ∂ p ∂ t c f − ∂ T ∂ t β ) (8)

同样地，

ρ o V = ρ o i [ 1 + c f ( p − p 0 ) − β ( T − T 0 ) ] ⋅ [ − k k r o μ o ( ∂ p ∂ x ) ] (9)

将上式带入连续方程左侧第二项得到：

∂ ( ρ o V ) ∂ x = ρ o i { ( ∂ p ∂ x c f − ∂ T ∂ x β ) ⋅ [ − k k r o μ o ( ∂ p ∂ x ) ] − [ 1 + c f ( p − p 0 ) − β ( T − T 0 ) ] ⋅ [ ∂ ∂ x ( k k r o μ o ) ( ∂ p ∂ x ) + k k r o μ o ∂ 2 p ∂ x 2 ] } (10)

由于压缩系数和热膨胀系数很小 [ 1 + c f ( p − p 0 ) − β ( T − T 0 ) ] ≈ 1 ，上式简化为：

∂ ( ρ o V ) ∂ x = ρ o i { ( ∂ p ∂ x c f − ∂ T ∂ x β ) ⋅ [ − k k r o μ o ( ∂ p ∂ x ) ] − [ ∂ ∂ x ( k k r o μ o ) ( ∂ p ∂ x ) + k k r o μ o ∂ 2 p ∂ x 2 ] } (11)

压缩系数和热膨胀系数远远小于1 [5]，上式可继续化简为：

∂ ( ρ o V ) ∂ x = − ρ o i [ ∂ ∂ x ( k k r o μ ) ( ∂ p ∂ x ) + k k r o μ ∂ 2 p ∂ x 2 ] (12)

因此，连续方程转化为：

ϕ ρ o i ( ∂ p ∂ t c f − ∂ T ∂ t β ) = ρ o i [ ∂ ∂ x ( k k r o μ o ) ( ∂ p ∂ x ) + k k r o μ o ∂ 2 p ∂ x 2 ] (13)

上式即为非等温条件下的压强扩散方程。

2.2. Warm VAPEX压强扩散推导

由质量守恒得到连续方程为：

∂ ( ρ ϕ ) ∂ t + ∂ ( ρ V ) ∂ x = 0 (14)

其中 ϕ 为孔隙度，无量纲，本文假设孔隙度为常数； ρ 表示稠油和溶剂混合流体的密度， kg / m 3 ； V 表示对流速度，可根据达西公式计算得到：

V = − k k r o μ ( ∂ p ∂ x ) (15)

其中 k 为绝对渗透率，达西； k r o 为油相相对渗透率，无量纲； μ 为稠油和溶剂混合物的粘度， Pa ⋅ s ； p 表示压强， Pa ；石油饱和度和浓度之间的关系采用Sharma和Gates的假设，而相对渗透率和沥青饱和度之间的相关性参考Corey的相关性[6]：

S o = S o r + ( S o i − S o r ) ( 1 − c c * ) (16)

k r o = k r o , r w ( S o − S o r 1 − S w c − S o r ) n (17)

其中 c 为溶剂浓度， m 3 /m 3 ； S o i 和 S o r 分别是初始和残余油饱和度，无量纲； S w c 是原生水饱和度，无量纲； k r o , r w 是最小含水饱和度下的石油相对渗透率，无量纲； n 是Corey系数，无量纲。

非等温条件下的状态方程：

ρ o = ρ o i [ 1 + c f 1 ( p − p 0 ) − β 1 ( T − T 0 ) ] (18)

ρ s = ρ s i [ 1 + c f 2 ( p − p 0 ) − β 2 ( T − T 0 ) ] (19)

得到稠油和溶剂混合流体的密度为：

ρ = ( 1 − c ) ρ o + c ρ s (20)

将(17) (18)带入(19)得：

ρ = ( 1 − c ) ρ o i [ 1 + c f 1 ( p − p 0 ) − β 1 ( T − T 0 ) ] + c ρ s i [ 1 + c f 2 ( p − p 0 ) − β 2 ( T − T 0 ) ] (21)

因此：

ρ ϕ = ϕ ρ = ϕ { ( 1 − c ) ρ o i [ 1 + c f 1 ( p − p 0 ) − β 1 ( T − T 0 ) ] + c ρ s i [ 1 + c f 2 ( p − p 0 ) − β 2 ( T − T 0 ) ] } (22)

带入连续方程等号左侧第一项得：

∂ ( ρ ϕ ) ∂ t = ϕ [ − ∂ c ∂ t ρ o i [ 1 + c f 1 ( p − p 0 ) − β 1 ( T − T 0 ) ] + ( 1 − c ) ρ o i ( ∂ p ∂ t c f 1 − ∂ T ∂ t β 1 ) ] + ϕ [ ∂ c ∂ t ρ s i [ 1 + c f 2 ( p − p 0 ) − β 2 ( T − T 0 ) ] + c ρ s i ( ∂ p ∂ t c f 2 − ∂ T ∂ t β 2 ) ] (23)

由于压缩系数和热膨胀系数很小 [ 1 + c f 1 ( p − p 0 ) − β 1 ( T − T 0 ) ] ≈ 1 ， [ 1 + c f 2 ( p − p 0 ) − β 2 ( T − T 0 ) ] ≈ 1 ，所以上式可简化为：

∂ ( ρ ϕ ) ∂ t = ϕ [ − ∂ c ∂ t ρ o i + ( 1 − c ) ρ o i ( ∂ p ∂ t c f 1 − ∂ T ∂ t β 1 ) ] + ϕ [ ∂ c ∂ t ρ s i + c ρ s i ( ∂ p ∂ t c f 2 − ∂ T ∂ t β 2 ) ] (24)

进一步整理得到：

∂ ( ρ ϕ ) ∂ t = ϕ ( ρ s i − ρ o i ) ∂ c ∂ t + [ ( 1 − c ) c f 1 ρ o i + c f 2 c ρ s i ] ϕ ∂ p ∂ t − [ ( 1 − c ) β 1 ρ o i + β 2 c ρ s i ] ϕ ∂ T ∂ t (25)

同样地，

ρ V = { ( 1 − c ) ρ o i [ 1 + c f 1 ( p − p 0 ) − β 1 ( T − T 0 ) ] + c ρ s i [ 1 + c f 2 ( p − p 0 ) − β 2 ( T − T 0 ) ] } ⋅ [ − k k r o μ ( ∂ p ∂ x ) ] (26)

将上式带入连续方程左侧第二项得到：

∂ ( ρ V ) ∂ x = { ∂ c ∂ x ρ o i [ 1 + c f 1 ( p − p 0 ) − β 1 ( T − T 0 ) ] − ( 1 − c ) ρ o i ( ∂ p ∂ x c f 1 − ∂ T ∂ x β 1 ) } ⋅ k k r o μ ( ∂ p ∂ x ) − { ∂ c ∂ x ρ s i [ 1 + c f 2 ( p − p 0 ) − β 2 ( T − T 0 ) ] + c ρ s i ( ∂ p ∂ x c f 2 − ∂ T ∂ x β 2 ) } ⋅ k k r o μ ( ∂ p ∂ x ) − { ( 1 − c ) ρ o i [ 1 + c f 1 ( p − p 0 ) − β 1 ( T − T 0 ) ] + c ρ s i [ 1 + c f 2 ( p − p 0 ) − β 2 ( T − T 0 ) ] } ⋅ { ∂ ∂ x ( k k r o μ ) ( ∂ p ∂ x ) + k k r o μ ∂ 2 p ∂ x 2 } (27)

简化为：

∂ ( ρ V ) ∂ x = { ∂ c ∂ x ρ o i − ( 1 − c ) ρ o i ( ∂ p ∂ x c f 1 − ∂ T ∂ x β 1 ) } ⋅ k r o k μ ( ∂ p ∂ x ) − { ∂ c ∂ x ρ s i + c ρ s i ( ∂ p ∂ x c f 2 − ∂ T ∂ x β 2 ) } ⋅ k r o k μ ( ∂ p ∂ x ) − [ ( 1 − c ) ρ o i + c ρ s i ] ⋅ { ∂ ∂ x ( k k r o μ ) ( ∂ p ∂ x ) + k k r o μ ∂ 2 p ∂ x 2 } (28)

进一步整理得到：

∂ ( ρ V ) ∂ x = { ∂ c ∂ x ( ρ o i − ρ s i ) − [ ( 1 − c ) ρ o i c f 1 + c ρ s i c f 2 ] ∂ p ∂ x + [ ( 1 − c ) ρ o i β 1 + c ρ s i β 2 ] ∂ T ∂ x } ⋅ k r o k μ ( ∂ p ∂ x ) − [ ( 1 − c ) ρ o i + c ρ s i ] ⋅ { ∂ ∂ x ( k k r o μ ) ( ∂ p ∂ x ) + k k r o μ ∂ 2 p ∂ x 2 } (29)

由于稠油和溶剂的压缩系数以及热膨胀系数均远远小于1，所以上式继续化简为：

∂ ( ρ V ) ∂ x = { ∂ c ∂ x ( ρ o i − ρ s i ) } ⋅ k r o k μ ( ∂ p ∂ x ) − [ ( 1 − c ) ρ o i + c ρ s i ] ⋅ { ∂ ∂ x ( k k r o μ ) ( ∂ p ∂ x ) + k k r o μ ∂ 2 p ∂ x 2 } (30)

因此，连续方程转化为：

ϕ ( ρ s i − ρ o i ) ∂ c ∂ t + [ ( 1 − c ) c f 1 ρ o i + c f 2 c ρ s i ] ϕ ∂ p ∂ t − [ ( 1 − c ) β 1 ρ o i + β 2 c ρ s i ] ϕ ∂ T ∂ t = − { ∂ c ∂ x ( ρ o i − ρ s i ) } ⋅ k r o k μ ( ∂ p ∂ x ) + [ ( 1 − c ) ρ o i + c ρ s i ] ⋅ { ∂ ∂ x ( k k r o μ ) ( ∂ p ∂ x ) + k k r o μ ∂ 2 p ∂ x 2 } (31)

上式即为非等温条件下考虑溶剂影响的压强扩散方程。

3. 对比等温与非等温压强扩散方程

以SAGD技术为例，对比等温与非等温压强扩散方程之间的差异。文献中等温压强扩散方程为[2]：

∂ ∂ x ( k r o μ o ) ⋅ ( ∂ p ∂ x ) + k r o μ o ⋅ ∂ 2 p ∂ x 2 = ϕ C t k ⋅ ∂ p ∂ t (32)

对比方程(13)与方程(32)发现，本文建立的非等温压强扩散方程与文献中等温压强扩散方程相比，两个方程均存在压强对空间的一阶和二阶偏导项以及压强对时间的一阶偏导项，但方程(13)多了温度对时间的偏导项，并且本文稠油密度为温度和压强的函数，而文献中将稠油密度假设为常数。通过模型对比，可以说明本文模型更为全面精准。

由于稠油被蒸汽携带的热量加热导致粘度降低，流动性增强，沿过渡带流至生产井，因此溶剂腔与稠油接触面会不断更新，从而会导致控制方程的边界条件发生变化。为此引入移动坐标系[7]：

{ τ = t      ξ = x − U t (33)

进而得到：

∂ ∂ x = ∂ ∂ ξ ⋅ ∂ ξ ∂ x = ∂ ∂ ξ (34)

∂ 2 ∂ x 2 = ∂ ( ∂ ∂ x ) ∂ x = ∂ ( ∂ ∂ ξ ) ∂ x = ∂ ( ∂ ∂ ξ ) ∂ ξ ⋅ ∂ ξ ∂ x = ∂ 2 ∂ ξ 2 (35)

∂ ∂ t = ∂ ∂ ξ ⋅ ∂ ξ ∂ t + ∂ ∂ τ ⋅ ∂ τ ∂ t = − U ∂ ∂ ξ + ∂ ∂ τ (36)

将(32) (33)和(34)带入方程(12)得：

ϕ ρ o i [ ( − U ∂ p ∂ ξ + ∂ p ∂ τ ) c f − ( − U ∂ T ∂ ξ + ∂ T ∂ τ ) β ] = ρ o i [ ∂ ∂ ξ ( k k r o μ o ) ( ∂ p ∂ ξ ) + k k r o μ o ∂ 2 p ∂ ξ 2 ] (37)

初边值条件为：

p ( 0 , τ ) = 35 0 kpa (38)

p ( L , τ ) = 30 0 kpa (39)

p ( ξ , 0 ) = 30 0 kpa (40)

移动坐标系下的传热方程为：

( 1 − ϕ ) ρ s c p s ( − U ∂ T ∂ ξ + ∂ T ∂ τ ) + ϕ c p o [ ρ o ( − U ∂ T ∂ ξ + ∂ T ∂ τ ) + T ( − U ∂ ρ o ∂ ξ + ∂ ρ o ∂ τ ) ] = K t ∂ 2 T ∂ ξ 2 (41)

温度的初边值条件如下：

T ( 0 , τ ) = T s (42)

T ( L , τ ) = 20 ˚ C (43)

T ( ξ , 0 ) = 20 ˚ C (44)

其中 T s 表示注入温度，℃； ρ s 为固体骨架的密度， kg / m 3 ； ρ o 为稠油密度， kg / m 3 ； c p s 和 c p o 分别为固体骨架和稠油的比热容， J/kg ⋅ K ； K t 为总导热系数， W/ ( m ⋅ ˚ C ) 。

如图3所示，利用有限差分法求解压强扩散和传热方程，得到第30天不同温度下不同模型的压强分布曲线。从图中可以看到，不同温度下，两模型计算得到的压强分布曲线趋势相同，压强都会随着空间距离的增加而下降，但是下降快慢程度有所差别。本文模型得到的压强随着空间距离的增加，下降趋势由快变慢，而等温压强扩散模型得到的压强分布正好相反，压强随着空间距离的增加，下降趋势由慢变快，并且温度越高，两者压强分布的差异越大，说明等温与非等温压强扩散方程计算得到的结果会有差异，而在有温度作用的采油技术中，有必要建立非等温压强扩散模型。

Figure 3. Comparison of pressure distributions calculated using the different pressure diffusion equation at various temperatures

图3. 不同温度下不同压强扩散方程计算所得压强分布对比

此外，将本文非等温压强扩散方程与传热方程耦合求解，还可以得到泄油速度 V ′ ：

V ′ = k k r o μ o ρ o g sin θ (45)

进而可以计算出溶剂腔双侧泄油速率 q ：

q = 2 ϕ W ∫ V d ξ (46)

因此通过数学建模可以预测不同开采条件下的泄油速率，进而辅助和指导实际开采。

4. 总结

本文创新性推导了SAGD技术下的非等温压强扩散方程以及warm VAPEX技术下与浓度有关的非等温压强扩散方程。并以SAGD技术为例，对比等温与非等温压强扩散方程及其计算结果的差异，得出以下主要结论：

(1) 通过以SAGD技术为例，对比等温与非等温压强扩散方程计算结果的差异，发现两者计算得到的压强分布虽变化趋势近似，但是两者之间的差异会随着温度升高而逐渐增大。

(2) 本文研究结果说明，在有温度作用的采油技术建模分析中，有必要建立非等温压强扩散方程，进而结合其他方程计算出泄油速度进而求得泄油速率，方便对不同条件下的采油技术进行预测和研究。

基金项目

中国国家自然科学基金(编号：21878018、22178022)；北京市教育委员会科学研究计划项目资助(KM202210016001)；北京建筑大学青年教师科研能力提升计划资助(X21030)。

NOTES

*通讯作者。

参考文献

[1] 王学慧, 代玉杰, 赵阳. 稠油提高采收率技术现状及发展趋势[J]. 现代化工, 2024, 44(10): 34-38, 43.
[2] Jia, X., Qu, T., Chen, H. and Chen, Z. (2019) Transient Convective Heat Transfer in a Steam-Assisted Gravity Drainage (SAGD) Process. Fuel, 247, 315-323.
https://doi.org/10.1016/j.fuel.2019.03.022
[3] Jia, X., Li, J., Lin, R. and Chen, Z. (2020) Mathematical Modeling of Dynamic Mass Transfer in Cyclic Solvent Injection. Journal of Petroleum Science and Engineering, 184, Article ID: 106573.
https://doi.org/10.1016/j.petrol.2019.106573
[4] Brooks, R.H. and Corey, A.T. (1964) Hydraulic Properties of Porous Media. Hydrology Papers, No. 3, Colorado State University Fort Collins, Colorado.
[5] 孔祥言. 高等渗流力学[M]. 合肥: 中国科学技术大学出版社, 1999.
[6] Corey, A.T. (1954) The Interrelation between Gas and Oil Relative Permeabilities. Prod Monthly, 19, 38-41.
[7] 范杰, 李相方. 蒸汽辅助重力泄油蒸汽腔前缘传热模型研究[J]. 科学技术与工程, 2016, 16(3): 42-47, 65.

Baidu
map