重质原油在圆管内由变压力梯度导致的流动及热流固耦合分析

Analysis of Flow and Thermal Fluid Solid Coupling of Heavy Crude Oil in a Circular Tube Caused by Variable Pressure Gradient

摘要: 原油作为重要的能源和战略资源，在日常生活、工业、农业等领域有着广泛的应用，因此原油的需求量持续增加。为了提高石油采收率(EOR)，深入研究重质原油在管道中的流动传热问题具有重要意义。本文探讨了热水开采时抽油泵内重质原油的流动传热问题。根据稠油的流变学实验数据，证实重质原油符合幂律流体本构方程；建立了重质原油在杆式抽油泵内的流动传热模型；首次采用有限差分方法求解了幂律方程的数值解，分析了杆式抽油泵热采温度T_w、抽油泵壁的厚度c、杆式抽油泵冲程次数n₁以及幂律指数n对流量Q的影响，建立了杆式抽油泵热采重质原油的最优方案。

Abstract: As an important energy and strategic resource, crude oil has a wide range of applications in daily life, industry, agriculture and other fields, so the demand for crude oil continues to increase. It is of great significance to conduct in-depth research on the flow and heat transfer of heavy crude oil in pipelines in order to improve oil recovery (EOR). This article explores the flow and heat transfer of heavy crude oil in the oil pump during hot water extraction. According to the rheological experimental data of heavy oil, it is confirmed that heavy crude oil conforms to the power-law fluid constitutive equation; A flow and heat transfer model for heavy crude oil in a rod pump was established; For the first time, the finite difference method was used to solve the numerical solution of the power-law equation. The effects of the thermal recovery temperature T_w of the rod pump, the thickness c of the pump wall, the number of strokes n₁ of the rod pump, and the power-law exponent n on the flow rate Q were analyzed, and the optimal scheme for thermal recovery of heavy crude oil by the rod pump was established.

1. 引言

流变学是研究物质变形与流动的科学，原油的流变学性质取决于原油的组成[1] [2]，即取决于原油中溶解气、液体和固体物质的含量，了解原油的流变特性有利于提高原油的采收率。稠油通常表现出明显的非牛顿流变特性。李一波等人[3]综合研究了温度、压力和剪切速率等情况对原油流变性的影响，发现低温下的粘度剪切速率数据与幂律流体本构方程吻合良好。在此基础上，Wang等人[4]更加系统地探讨了温度对含蜡原油流变特性的影响，通过研究表明，在双板浮顶油罐中，当油温高于蜡沉淀点时，原油表现出牛顿行为；当油温降至蜡沉淀点以下时，蜡晶体开始沉淀，但此时含蜡原油仍表现出牛顿行为。当油温低于异常点时，它开始表现出非牛顿特性，并用幂律流体来描述非牛顿流体特性。此外，Wang等人[5]主要分析了幂律指数对原油产能的影响，通过研究发现，幂律指数n小于0.8时，产量随幂律指数的增加而缓慢增加，当幂律系数n大于0.8时，生产速度会大大增大。随后，Mohammadi等人[6]通过MCR 302流变仪生成原油流量曲线，证明了该曲线符合幂律模型。郭等人[7]将原油及其乳化液视为牛顿和幂律流体，主要探究了垂直管道内压力，流速和管径对流体表观粘度的影响。除此之外，在驱油过程中，水驱动和聚合物驱动也会对原油的流变特性产生影响。Kamyabi和Ramazani [8]将微孔中的驱油剂看成不同幂律指数值的稳态幂律流体来研究对原油的驱动作用，结果表明幂律指数n大于1的驱油剂比n小于1的驱油剂驱油效果好。最后，Zhao和Min [9]为了更好地描述油田开发后期高含水原油的流变特性，建立了具有弹性外边界条件的非牛顿幂律流体渗流模型。更多原油的非牛顿流变特性可见如下文献[10]-[13]。

近年来，如何提高石油采收率受到了广泛关注。其中，通过降低原油粘度来提高流动性的热采方法已成为研究热点。上世纪末，Dang [14]发现，忽略管道输送中壁面厚度对传热的影响会使数值计算结果产生较大偏差，于是建立了热–流体–固体耦合模型。Oliveski [15]通过对储罐中的原油进行分层分析，使用有限体积法获得了原油的温度和流场，而Sun等人[16]认为不同的盘管结构对原油温度的影响最大，并使用CFD方法研究了储罐加热中原油的传热和流动的耦合特性。Hao等人[17]考虑了轴承运行期间油和轴承之间的热流固耦合效应，并通过实验验证发现该模型更准确。此外，Monge和Birken [18]建立了水平管道中一维非稳态动量和热流固耦合的能量方程，在时间和空间上选择了不同的数值模拟方法，并验证了算法的有效性。然而，很少有学者使用热流固耦合模型来研究垂直井管道的传热情况对原油温度及速度的影响。

2. 抽油泵内变压力梯度表达式的推导

考虑了一种新的热采方法，将热水注入垂直生产井的同轴圆筒中，通过加热泵壁，降低原油粘度，提高采油率。本文考虑了原油在泵筒内的非稳态径向流动和热流固耦合。有杆泵的物理模型示意图如图1所示，x轴为轴向坐标，r轴垂直于管道。

Figure 1. Physical model of rod pump

图1. 抽油泵物理模型示意图

图1(a)中，抽油泵位于下死点(最低点)，为了防止柱塞运动到最下方时撞击到固定阀阀罩，因此抽油泵的下死点和固定阀阀罩之间有一定的距离L，这个距离就是防冲距。此时，固定阀和游动阀都是关闭的状态。泵筒内气体压强Pa等于抽油泵的出口压强(游动阀上方的压强) P₀。泵筒内残余原油的高度为αL，其中α为原油的含气率[19]。

随着抽油杆带动柱塞向上移动，泵筒内气体压强Pa不断减小，到达某一高度时，固定阀在其上下压强差的推动下恰好开始向上移动，我们称此时的高度为无效抽汲高度h_b，如图1(b)所示。根据理想气体状态方程，图1(b)中泵筒内气体的压强表达式可表示为：

$P_b=\frac{n^{*}{R}^{*}T}{\pi R^2\left(h_b+\left(1-\alpha \right)L\right)}$ (1)

其中，n^*是气体的摩尔质量，R^*表示摩尔气体常数，T为热采的温度，R表示抽油泵内径的最大值。

此时，我们忽略固定阀的重力，结合伯努利方程，得到抽油泵筒内的气体压强P_b，抽油泵出口的压强(固定阀上方的压强) P₁及筒内残余原油的压强之间的函数关系式：

$P_b=P_1-{\rho }_{f}g\alpha L$ (2)

其中ρ_f表示原油的密度。

抽油杆柱塞的移动速度假设为：

$v\left(t^{*}\right)=v_0\sin \left(\frac{\pi }{t_0}{t}^{*}\right),0\le t^{*}\le t_0$ (3)

其中，v₀是抽油杆的最大速度，t₀ = 60/(2n₁)表示抽油杆单次冲程所需要的时间，n₁为冲次(每分钟抽油杆上下往复运动的次数)。

抽油泵的柱塞到达无效抽汲高度h_b所需要的时间为无效抽汲时间t₁，因此有效抽汲时间t₂可以表示为$t_2=t_0-t_1$ 。在有效抽汲时间内，以轴向速度流入抽油泵下泵腔内的原油高度为：

$x\left(t\right)=\frac{1}{R^2}{\displaystyle {\int }_0^t{\displaystyle {\int }_0^{R}ru\left(r,t\right)\text{d}r\text{d}t}},0\le t\le t_2$ (4)

本文在建立动量方程时主要考虑固定阀打开到关闭的这一段时间，即如图1(b)、图1(c)。于是，抽油杆柱塞的移动速度修正为：

$v\left(t\right)=v_0\sin \left(\frac{\pi }{t_0}\left(t+t_1\right)\right),0\le t\le t_2$ (5)

因此，在有效抽汲时间内柱塞移动的位移我们可以表示为：

$x_2\left(t\right)={\displaystyle {\int }_0^t{v}_0\sin \left(\frac{\pi }{t_0}\left(t+t_1\right)\right)\text{d}t}$ (6)

在不同时刻及不同柱塞移动位置处，抽油泵下泵腔内的气体压强可由理想气体状态方程获得，如下所示：

$p_{air}=\frac{n^{*}{R}^{*}T}{\pi R^2\left(h_b+\left(1-\alpha \right)L+x_2-x\right)}$ (7)

所以，通过计算可推理出下泵腔内原油的压强梯度表达式：

$\frac{\partial P}{\partial x}=-\frac{n^{*}{R}^{*}T}{\pi R^2}\frac{1}{{\left(h_b+\left(1-\alpha \right)L+\frac{H}{2}\left(\text{cos}\left(\frac{\pi }{t_0}{t}_1\right)-\cos \left(\frac{\pi }{t_0}\left(t+t_1\right)\right)\right)-\frac{1}{R^2}{\displaystyle {\int }_0^t{\displaystyle {\int }_0^{R}ru\left(r,t\right)\text{d}r\text{d}t}}\right)}^2}$ (8)

表1给出了抽油泵最大内半径R、防冲距L、重力加速度g和原油含气量α相关参数的取值[1]。

Table 1. Values of relevant parameters for oil pump

表1. 抽油泵相关参数取值

	H (cm)	T (K)	R (cm)	L (cm)	g (cm/s^−2)	α
Value	200	363	1.6	30	980	50%

3. 重质原油的流变学特性

稠油相对于轻油而言更需要热驱动，在山东胜利油田取适量样品，借用流变仪对其流变性进行研究，从图2中可以发现稠油符合幂律流体本构，并且幂律指数n < 1，具有剪切变稀流变特性。

同理，根据不同温度下重质原油的粘度值见表2，拟合出重质原油粘度–温度之间的函数关系式：

$\mu =10.7942{\text{e}}^{-0.097\left(T_f-293\right)}$ (9)

根据幂律流体本构方程和重质原油粘度–温度函数表达式，可以得出：

$\tau =10.7942{\text{e}}^{-0.097\left(T_f-293\right)}{\left|\frac{\partial u}{\partial r}\right|}^{n-1}\dot{\epsilon }$ (10)

Figure 2. Shear stress-shear rate curve of heavy crude oil

图2. 重质原油的剪切速率–剪切应力曲线

Table 2. Dynamic viscosity of heavy crude oil at different temperatures

表2. 重质原油在不同温度下的动态粘度

温度(K)	303	323	333	343
粘度(Pa∙s)	4.0923	0.5635	0.2396	0.1346

4. 重油流动与热流固耦合模型的建立

重油在抽油泵内的流动环境与轻油保持一致，在此基础上，建立重质原油在抽油泵内的非定常径向流动和热流固耦合模型：

$\begin{array}{c}\frac{\partial u}{\partial t}=\frac{n^{*}{R}^{*}T}{\pi R^2}\frac{1}{{\left(h_b+\left(1-\alpha \right)L+\frac{H}{2}\left(\text{cos}\left(\frac{\pi }{t_0}{t}_1\right)-\cos \left(\frac{\pi }{t_0}\left(t+t_1\right)\right)\right)-\frac{1}{R^2}{\displaystyle {\int }_0^t{\displaystyle {\int }_0^{R}ru\left(r,t\right)\text{d}r\text{d}t}}\right)}^2}\\ \text{+}\frac{\mu \left(T_f\right)}{{\rho }_f}{\left|\frac{\partial u}{\partial r}\right|}^{n-1}\left(\frac{{\partial }^{2}u}{\partial r^2}+\frac{1}{r}\frac{\partial u}{\partial r}\right)-g\end{array}$ (11)

$\frac{\partial T_f}{\partial t}={\alpha }_f{\left|\frac{\partial T_f}{\partial r}\right|}^{n-1}\left(\frac{{\partial }^2{T}_f}{\partial r^2}+\frac{1}{r}\frac{\partial T_f}{\partial r}\right)$ (12)

$\frac{\partial T_s}{\partial t}={\alpha }_s\left(\frac{{\partial }^2{T}_s}{\partial r^2}+\frac{1}{r}\frac{\partial T_s}{\partial r}\right)$ (13)

边界条件为：

$T\left(r,0\right)=T_0,u\left(r,0\right)=0,0\le r\le R$ (14)

$\frac{\partial T_f\left(0,t\right)}{\partial r}=0,\frac{\partial u\left(0,t\right)}{\partial r}=0$ (15)

${\lambda }_f\frac{\partial T_f\left(R,t\right)}{\partial r}={\lambda }_s\frac{\partial T_s\left(R,t\right)}{\partial r},T_f\left(R,t\right)=T_s\left(R,t\right),u\left(R,t\right)=0$ (16)

$T_s\left(R+c,t\right)=T_w$ (17)

同样根据有限差分方法，利用$r_i=i\Delta r\left(i=0,1,\cdots ,M\right)$ , $t_k=k\Delta t\left(k=0,1,\cdots ,K\right)$ , $\Delta r=R/M$ , $\Delta t=t_2/K$ 将计算区域分割成网格点，对模型进行离散，结果如下：

$\begin{array}{l}\left[\frac{1}{\Delta t}-\frac{\mu \left(T_f{}_i^k\right)}{{\rho }_f{r}_i\Delta r}B\left(i,k-1\right)+\frac{2\mu \left(T_f{}_i^k\right)}{{\rho }_f{\left(\Delta r\right)}^2}B\left(i,k-1\right)\right]u_i^k\\ =-g+A\left(t\right)+\frac{\mu \left(T_f{}_i^k\right)}{{\rho }_f{\left(\Delta r\right)}^2}B\left(i,k-1\right)u_{i+1}^k+\left[\frac{\mu \left(T_f{}_i^k\right)}{{\rho }_f{\left(\Delta r\right)}^2}-\frac{\mu \left(T_f{}_i^k\right)}{{\rho }_f{r}_i\Delta r}\right]B\left(i,k-1\right)u_{i-1}^k+\frac{1}{\Delta t}{u}_i^{k-1}\end{array}$ (18)

$\begin{array}{l}\left[\frac{1}{\Delta t}+\frac{2{\alpha }_f}{{\left(\Delta r\right)}^2}C\left(i,k-1\right)-\frac{{\alpha }_f}{r_i\Delta r}C\left(i,k-1\right)\right]T_f{}_i^k\\ =\frac{{\alpha }_f}{{\left(\Delta r\right)}^2}\left(i,k-1\right)T_f{}_{i+1}^k+\frac{1}{\Delta t}{T}_f{}_i^{k-1}+{\alpha }_{f}C\left(i,k-1\right)\left[\frac{1}{{\left(\Delta r\right)}^2}{T}_f{}_{i-1}^k-\frac{1}{r_i\Delta r}\right]T_f{}_{i-1}^k\end{array}$ (19)

$\left[\frac{1}{\Delta t}+\frac{2{\alpha }_s}{{\left(\Delta r\right)}^2}-\frac{{\alpha }_s}{r_i\Delta r}\right]T_s{}_i^k=\frac{{\alpha }_s}{{\left(\Delta r\right)}^2}{T}_s{}_{i+1}^k+{\alpha }_s\left[\frac{1}{{\left(\Delta r\right)}^2}{T}_s{}_{i-1}^k-\frac{1}{r_i\Delta r}\right]T_s{}_{i-1}^k+\frac{1}{\Delta t}{T}_s{}_i^{k-1}$ (20)

$Q=\frac{\Delta t}{2}\frac{\Delta r}{2}\left(u_0^0+u_0^K+2{\displaystyle \sum _{j=1}^{K-1}{u}_0^j}+u_M^0+u_M^K+2{\displaystyle \sum _{j=1}^{K-1}{u}_M^j}+2{\displaystyle \sum _{i=1}^{M-1}{u}_i^0}+2{\displaystyle \sum _{i=1}^{M-1}{u}_i^K}+4{\displaystyle \sum _{j=1}^{K-1}{\displaystyle \sum _{i=1}^{M-1}{u}_i^j}}\right)$ (21)

其中，

$A\left(t\right)=\frac{n^{*}{R}^{*}{T}_f{}_i^k}{\pi {\rho }_f{R}^2{\left\{\begin{array}{l}{h}_b+\left(1-\alpha \right)L+\frac{H}{2}\left(\text{cos}\left(\frac{\pi }{t_0}{t}_1\right)-\cos \left(\frac{\pi }{t_0}\left(t_k+t_1\right)\right)\right)\\ -\frac{\Delta t{\left(\Delta r\right)}^2}{R^2}\left[2{\displaystyle \sum _{m=2}^{k-2}{\displaystyle \sum _{p=2}^M\left(p-1\right)u\left(p,m\right)}}+{\displaystyle \sum _{p=2}^M\left(p-1\right)u\left(p,k-1\right)}\right]\end{array}\right\}}^2}$ (22)

$B\left(i,k-1\right)={\left|\frac{u_i^{k-1}-u_{i-1}^{k-1}}{\Delta r}\right|}^{n-1}$ (23)

$C\left(i,k-1\right)={\left|\frac{T_f{}_i^{k-1}-T_f{}_{i-1}^{k-1}}{\Delta r}\right|}^{n-1}$ (24)

5. 结果与讨论

当P₀ = 8 MPa，P₁ = 2 MPa时，主要讨论冲次n₁，壁面厚度c，热采温度T_w，幂律指数n对稠油流量Q的影响。图3表明冲次n₁对稠油流动速度u的影响，具体指在不同时刻，轴心处的速度在不同冲次下的流动曲线。从图中可以看出，冲次不但会影响重质原油的最大流动速度，还会影响原油产生倒流的时间。在初始时刻，原油流入抽油泵的速度大于抽油杆带动柱塞的移动速度，因此下泵腔内原油体积不断增加，空气体积减小，与进口压强之间的差值减弱，随后原油流速开始不断下降，甚至还会出现倒流。

Figure 3. Impact of stroke times n₁ on speed u of heavy crude oil

图3. 冲次n₁对稠油流速u的影响

图4揭示了幂律指数n对原油流动速度的影响，根据实验验证了此类原油的幂律指数0.9 < n < 1，具有剪切变稀的特性。从图中看出，随着幂律指数的增大，原油流动速度逐渐减弱，在n = 0.9时，流速流量同时取得最大值。

Figure 4. Impact of power-law index n on speed u of heavy crude oil

图4. 幂律指数n对稠油流速u的影响

从图5(a)中可以得出结论，热采温度T_w对原油的温度T_f有显著影响，通过加热抽油泵的外壁，使得管内原油温度也随之升高。图5(b)揭示了在轴心处不同时刻热采温度对抽油泵壁面和原油温度的影响，从图中可以发现，在不同时刻，起初原油温度受热采温度影响较小，1 s后热采温度对原油温度的影响逐渐增大。

(a) (b)

Figure 5. Influence of thermal recovery temperature T_w on wall and fluid temperature T_f

图5. 热采温度T_w对抽油泵壁面和流体温度T_f的影响

如图6(a)主要表示了热采温度对原油流动速度的影响，从图中可以看出，原油的流动速度会随着温度的增大而增大，主要原因是温度降低了原油的粘度，从而增大了原油在管道中的流动速度。图6(b)可以得出363 K是热水开采的最佳温度，此时，原油流速最大。

(a) (b)

Figure 6. Influence of thermal recovery temperature T_w on heavy crude oil velocity u

图6. 热采温度T_w对重油速度u的影响

抽油泵壁根据厚度不同分为普通壁和薄壁。图7显示了壁厚对原油流动的影响。结果表明，壁厚对泵内原油的流动有明显的阻碍作用。壁厚通过影响固体壁的导热性来影响原油的温度，从而影响原油粘度和速度。因此，抽油泵壁面的厚度在技术参数允许的范围内可以尽可能地减小。

6. 结论

本文基于稠油的流变特性和推导的压强梯度表达式，建立了重质原油在圆筒内的非稳态动量和热流固耦合模型，首次使用有限差分方法求解了符合幂律行为的数值解，主要结论如下：

1) 实验结果表明，东营胜利油田的重质原油属于剪切变稀的幂律流体，且幂律指数的取值范围为：0.9 < n < 1。

(a) (b)

Figure 7. Influence of wall thickness c on crude oil velocity u and flow Q

图7. 壁厚c对原油速度u和流量Q的影响

2) 对于符合幂律本构关系的原油，其在柱面坐标系下的动量和能量方程可以采用有限差分方法求解，结论与牛顿流体保持一致。

3) 原油流速随着幂律指数的增大而减小。壁面厚度会阻碍热量的传递，因此导致速度减小，然而对于热采温度则呈现相反的结果，热采温度会促进传热过程，加速流体流动。

基金项目

国家自然科学基金(Nos. 22478028, 22178022)。

NOTES

^*通讯作者。