SAGD循环预热注汽参数影响规律数值模拟

2019-04-17李朋张艳玉孙晓飞刘洋谢孟珂王朝陈会娟

中南大学学报（自然科学版） 2019年11期

李朋，张艳玉，孙晓飞，刘洋，谢孟珂，王朝，陈会娟

(1.中国石油大学(华东)石油工程学院，山东青岛，266580；2.中国石油勘探开发研究院，北京，100083；3.上海大学上海市应用数学和力学研究所，上海，200072)

稠油油藏热采技术主要包括蒸汽吞吐[1]、蒸汽驱[2−3]、蒸汽辅助重力驱油(steam assisted gravity drainage，SAGD)[4−6]等，其中SAGD 技术是实现稠油油藏高效开发的一项前沿技术。SAGD过程主要分为循环预热和正式生产[7]2个阶段，其中，循环预热过程井筒沿程地层预热效果对SAGD 转生产后油藏的整体开发有着重要的影响。由于注汽参数的选择与地层预热效果密切联系，因此，研究不同注汽参数条件下井筒沿程关键参数分布规律以及地层预热效果对合理优化注汽参数和提高地层预热均匀程度具有重要意义。目前，与注汽参数影响规律研究相关的模型主要分为解析模型和数值模型。在解析模型研究领域，DONG等[8]研究了不同注汽速度条件下常规单管井筒压力分布规律；WU 等[9]模拟了注汽速度对平行双管水平井井筒内蒸汽干度分布的影响；林日亿等[10]分析了不同蒸汽干度和注汽速度对水平井散热量的影响；SUN 等[11]针对同心双管射孔完井管柱结构特点，预测了不同注汽压力、蒸汽干度条件下井筒内热物理参数分布规律。上述模型仅考虑井筒内流体的流动与传热规律，而割缝筛管管柱结构特点对表皮系数、摩擦力、摩擦力做功和井筒传热的影响以及注汽过程中油藏物性参数随时间变化对井筒内流体流动及传热规律的影响均未考虑。在数值模型研究领域，人们将油藏物性参数随时间变化的影响考虑在内，分别建立了垂直井筒[12−13]、常规单管[14−16]和单管多点[17]数值模型以计算井筒内蒸汽热力参数，但该模型均不能用于模拟双管水平井井筒内流体的流动与传热规律。此外，现有的数值模拟软件如CMG中的灵活井模型虽然可以对双管水平井井筒内流体的流动和传热规律进行模拟，但其水平井井筒压降通过简单的机理模型进行计算[6]，且其模型中的井筒模型仅考虑裸眼完井和射孔完井管柱结构特点，不能用于模拟割缝筛管完井水平井井筒内流体的流动与传热规律。为此，本文作者以实际割缝筛管完井水平井管柱结构特点和SAGD 循环注汽方式为依据，分别建立SAGD循环预热过程中流体在井筒和地层内的流动与传热数学模型，通过编写计算机程序求解井筒内蒸汽热力参数，对不同注汽压力、蒸汽干度和注汽速度进行敏感性分析，优化注汽参数并对预热效果进行评价，以期为合理选择注汽参数、提高地层预热效果、实现稠油油藏高效开发提供参考。

1 数学模型建立

1.1 假设条件

SAGD循环预热阶段流体在井筒内流动示意图如图1所示。基于如下假设建立模型：1)预热阶段地层内为油、水、汽三相流，井筒内为汽液两相流；2)流体黏度是温度的函数，油、水、汽相对渗透率是饱和度的函数；3)长油管内流体流动为单向稳态流动，环空内流动为水平方向稳态流和径向非稳态流；4)在循环预热过程中，长油管跟端A点截面注汽压力、干度以及注汽速度均为常数，且流体在井筒内任意截面处蒸汽参数值相等；5)井筒内为稳态传热，井筒外缘到地层为非稳态传热；6)考虑割缝筛管参数对表皮系数、摩擦力、摩擦力做功以及井筒传热的影响。

图1 SAGD循环预热示意图Fig.1 Diagram of SAGD preheating process

1.2 井筒内流体流动与传热模型

1.2.1 流体流动模型

蒸汽在长油管内流动时，无质量流量损失，而环空内的蒸汽则在压差的作用下注入地层，因此，可得长油管和环空内蒸汽流动时的质量守恒方程：

式中：mlt和man分别为长油管和环空内流体质量流速，kg/s；manf为地层吸汽量，kg/(m·s)；l为微元段长度，m。

蒸汽在井筒内流动时压降主要来自摩擦阻力压降和加速度压降，因此，可得长油管和环空内蒸汽流动时的动量守恒方程：

式中：Alt和Aan分别为长油管和割缝筛管有效横截面积，m2；ρlt和ρan分别为长油管和环空内流体密度，流体密度采用Beggs−Brill 算法[18]，kg/m3；vlt和van分别为长油管和环空内流体流速，m/s；plt和pan分别为长油管和环空压力，Pa；τlt和τan分别为长油管和环空内壁摩擦阻力，N。

蒸汽在长油管和环空内流动过程中，蒸汽热损失和摩擦热损失等于流体内能的变化量加上注入地层蒸汽的热损失，根据能量守恒原理，可得长油管和环空内蒸汽流动时的能量守恒方程：

式中：Qlt和Qan分别为长油管和环空内流体的散热量，J/s;hlt和han分别为长油管和环空内流体的焓，J/kg；Wlt和Wan分别为摩擦力在长油管和环空微元段上做的功，J/s；var为蒸汽进入地层的速度，m/s。

1.2.2 流体传热模型

在SAGD 循环预热过程中，井筒内流体的传热计算公式可表示为：

式中：klt和kan分别为长油管和环空导热系数，W/(m·K)；Tlt，Tan和Tres分别为长油管、环空和地层温度，°C。

1.2.3 摩擦力及摩擦力做功计算

当长油管内蒸汽流动时，其摩擦力及摩擦力做功采用常规算法计算[10]。当蒸汽从长油管趾端进入环空时，考虑割缝筛管参数对摩擦力、摩擦力做功的影响，摩擦力及摩擦力做功计算式为[19]：

式中：pgf为微元段割缝排数；τanl(k)为第k排割缝筛管内壁摩擦力，N；vanl(k)为第k排割缝蒸汽质量流速，kg/s；lu为割缝单元长度，m；f为流体与管壁之间的摩阻系数；rani为割缝筛管环空微元段有效半径，m；ΔW(k)为摩擦力在第k排割缝筛管上做的功，J。

如果饱和蒸汽在流动过程中发生相变，冷凝成热水，井筒中的流动为单相流，那么在该情况下管壁摩阻系数与套管壁粗糙度、雷诺数和流体的流动型态相关，其计算可参考文献[20]中的单相管流摩阻系数计算公式。如果饱和蒸汽在流动过程中为汽液两相流，那么摩阻系数则采用经典的Beggs−Brill方法计算[18]。

1.3 地层内流体流动与传热模型

在单位时间内地层微元体的累积质量和能量增量等于单位时间内流入和流出微元体的质量差和能量差之和，地层内流体流动的质量守恒方程和能量守恒方程如下。

油相质量守恒方程：

水相质量守恒方程：

蒸汽相质量守恒方程：

能量守恒方程：

式中：σ为单位转换系数；φ为孔隙度；Pi为压力，i为o，w和g，分别表示油、水、蒸汽三相，Pa；ρi为密度，kg/m3；μi为黏度，Pa·s；mi为单位时间体积地层内增加或减少的流体质量(注入为“＋”，采出为“−”)，kg/(m3·s)；Si为饱和度；kri为相对渗透率，10−3µm2；Hi为热焓，J/kg；Ui为内能，J/kg；qc为单位时间单位体积内冷凝水的质量，kg/(m3·s)；g为重力加速度，m/s2；cr为岩石比热容，J/(kg·°C)；H为深度，m；ρr为岩石密度，kg/m3；Qliq为单位时间体积内产出或注入的流体能量(注入为“＋”，采出为“−”)，J/(m3·s)；λr为地层传热系数，W/(m·°C)；Qloss为单位时间单位体积顶底层吸收的热量，J/(m3·s)；Qwell为单位时间地层通过井筒热传导吸收的热量，J/s；ke为储层渗透率，10−3µm2；τ为时间，s。

式(11)～(14)中有7个未知数，故需联立辅助方程进行求解，辅助方程如下。

蒸汽−水相平衡方程[21]：

饱和度方程：

毛管力方程：

式中：θs为蒸汽温度，°C；ps为蒸汽压力，Pa；Pcow为油水两相毛管力，Pa；Pcgo油汽两相毛管力，Pa。

1.4 耦合模型

井筒模型和地层模型可通过注汽量公式进行耦合求解：

式中：mg,i为单位时间内注入地层的蒸汽量，kg/s；Jan,i为水平井指数，m3/(Pa·s)；Pan,i为环空蒸汽压力，Pa；ρan,i为微元段蒸汽密度，kg/m3；Pi为微元段地层压力，Pa；ke为各向异性等价渗透率，10−3µm2；re为泄油半径，m；r为井筒半径，m；Ia,i为微元段吸汽指数；S为表皮系数，考虑割缝筛管参数对表皮系数的影响，采用Furui算法计算[22]。

2 模型求解

通过编制程序求解方程(1)～(20)即可得到SAGD循环预热阶段流体在井筒和地层内的流动及传热规律，具体步骤如下：1)在一定井底流压和注汽速度下，假定每一时间步地层压力、温度及饱和度变化值；2)将长油管跟端A点到趾端B点，环空趾端B点到跟端A点均匀划分，并假定各微元段入口处的蒸汽压力、温度和干度；3)长油管第1个微元段出口处的压力、温度和干度以该微元段入口处的蒸汽参数值为起点，反复迭代计算，直至计算结果与假定的蒸汽参数值在误差允许范围内，则将该微元段出口处的蒸汽参数作为下一微元段入口处的蒸汽参数，依此类推，直到计算至长油管趾端B点；4)环空趾端B点到跟端A点第1个微元段出口处的压力、温度和干度计算以长油管趾端B点处的蒸汽参数为起点，求解思路与长油管内蒸汽参数计算方法一致，直到计算至环空跟端A点；5)采用全隐式有限差分法对地层内流体流动与传热模型进行处理，可得到大型七对角块状稀疏奇异矩阵，将步骤3)和4)中所求得的井筒内蒸汽参数代入进行耦合求解，反复迭代计算，直至求得的地层压力、温度及饱和度变化值与假定的地层各参数变化值在误差允许范围内；6)将步骤5)中所得各地层参数作为下一时间步的地层初始条件，重复步骤1)～5)进行迭代计算，直到运算至给定的模拟时间，计算结束。

3 计算实例

以某油田1对SAGD水平井组为例，模型计算所需基本参数如下：油藏埋深180 m，油层厚为20 m，孔隙度为0.33，横向渗透率为2.7µm2，垂向渗透率为1.89µm2，原始地层压力为0.22 MPa，原始地层温度为9°C，含油饱和度为0.75，油层热扩散系数为4.0×10−3m2/h，岩石导热系数为5.2 W/(m·°C)，岩石比热容为1.138 kJ/(kg·°C)，岩石压缩系数为2.0×10−3MPa−1，原油和水的热膨胀系数分别为4.38×10−4°C−1和1.45×10−4°C−1，原油和水的压缩系数分别为2.44×10−4MPa−1和5.21×10−4MPa−1，原油和水的比热容分别为2.2 kJ/(kg·°C)和4.2 kJ/(kg·°C)。水平井长度为850 m，注汽压力为2.0 MPa，注汽速度为130 t/d，蒸汽干度为95%，筛管外半径和内半径分别为109和105 mm，长油管外半径和内半径分别为44.45和38.95 mm，短油管外半径和内半径分别为36.5和31.0 mm；长油管和短油管的管壁绝对粗糙度为4.37×10−5m；割缝筛管的管壁绝对粗糙度为5.0×10−5m。

3.1 验证模型可靠性

为确保模型准确可靠，将模型蒸汽日注入量预测结果和CMG模拟结果与油田实测结果进行对比，如图2所示。由图2可知：模型蒸汽日注入量预测结果低于CMG模拟结果，且模型蒸汽日注入量预测结果更接近于实测结果，这是因为新建模型考虑了割缝筛管管柱结构特点，同时井筒压降计算考虑了滑脱的影响[18]，能更准确地描述SAGD循环预热过程中流体的流动及传热规律。模型预测结果与实测结果相比，蒸汽日注入量最大相对误差为7.21%，在误差允许范围内，说明模型是可靠的。

3.2 井筒沿程关键参数分布

图2 蒸汽日注汽量预测结果与模拟以及现场结果对比Fig.2 Comparison of predicted result with simulated and measured results of steam injection rate

在油田实际生产参数条件下，井筒沿程关键参数分布规律如图3所示。由图3可知：长油管内压力和温度降幅比环空内的显著，降幅分别为12.77 kPa和0.33°C(见图3(a)和图3(b))，环空内蒸汽干度降幅比长油管内的明显，降幅为0.61(见图3(c))。这是由于长油管内蒸汽质量流量为定值，摩擦损失较大，热损失较小，而当蒸汽经长油管趾端进入环空后，蒸汽不断注入地层，蒸汽质量流量逐渐减少，摩擦损失降低，热损失增加。此外，蒸汽干度损失主要发生在环空跟端，这是因为当蒸汽未返回到跟端时，井筒内蒸汽质量流量较大，蒸汽有充足的热能，对蒸汽干度影响不大。而当蒸汽返回到跟端时，井筒内蒸汽质量流量较少，蒸汽热能明显减小，对蒸汽干度影响较大。

研究表明，在SAGD 循环预热过程中，当环空温度在预热3 d内达到稳定且返回跟端蒸汽干度大于0时，预热效果较好。由于注汽参数对预热效果有重要影响，为探索不同注汽压力、蒸汽干度和注汽速度对SAGD 循环预热阶段井筒沿程关键参数分布和预热效果的影响，选取不同注汽压力(1.5～2.5 MPa)，蒸汽干度(0.55～0.95)和注汽速度(130～170 t/d)进行数值模拟研究。循环预热3 d，不同注汽参数条件下模拟结果见表1。

由表1可知：随注汽压力的增加，井筒蒸汽压力增大，返回跟端的蒸汽干度降低，趾端B点与跟端A点的地层吸热率差值增大(地层吸热率是指单位时间内单位长度井筒沿程地层吸收的热量)；随蒸汽干度和注汽速度的增加，井筒蒸汽压力下降，返回跟端的蒸汽干度升高，趾端B点与跟端A点的地层吸热率差值缩小。

表1 不同注汽参数条件下模拟结果Table 1 Simulation results with different steam injection parameters

3.3 注汽参数影响规律分析

为进一步研究注汽压力、蒸汽干度和注汽速度对井筒沿程关键参数分布以及地层预热效果的影响，以预热3 d返回跟端蒸汽干度大于0为目标，对比不同注汽参数条件下模拟结果并对其进行优化。

3.3.1 注汽压力敏感性分析

图4所示为当注汽速度为130 t/d，蒸汽干度为0.95时，不同注汽压力下井筒沿程关键参数分布规律。由图4可知：随注汽压力增加，井筒内蒸汽压力和地层吸热率均增大(见图4(a)和图4(c))，而返回跟端的蒸汽干度降低(见图4(b))，这是因为注汽压力越大，注入到地层的蒸汽和热量越多，井筒内蒸汽热损失越大。当注汽压力由1.5 MPa 增至2.5 MPa 时，蒸汽压降由21.67 kPa 降至20.99 kPa，减小0.68 kPa(见图4(a))；环空干度降由0.18 增至0.64，增加了0.46(见图4(b))；当环空趾端到跟端蒸汽干度差增大时，受井筒沿程蒸汽干度热焓差的影响，水平段地层从趾端到跟端沿程不同位置吸收的蒸汽热焓差别也就越大，井筒沿程地层吸热率最大值与平均值的比值由1.03 增至1.11(见图4(c))，地层吸热率分布不均匀程度增加，水平段地层均匀预热效果变差。因此，减小注汽压力有利于水平段地层均匀预热，但井筒沿程地层吸热率会相对减小。

图4 不同注汽压力下井筒沿程关键参数分布Fig.4 Distribution of key parameters along wellbore with different steam injection pressures

3.3.2 蒸汽干度敏感性分析

图5所示为当注汽压力为2.0 MPa，注汽速度为170 t/d 时，不同蒸汽干度下井筒沿程关键参数分布规律。由图5可知：随着蒸汽干度增加，井筒内蒸汽压力降低(见图5(a))，这是因为蒸汽干度越大，其密度越小，流速越大，井筒摩擦损失也就越大；而蒸汽干度和地层吸热率均随蒸汽干度增加而上升(见图5(b)和图5(c))，这是因为蒸汽干度越高，其焓越大，地层吸收的热量越多。当蒸汽干度由0.55 增至0.95 时，井筒压降由19.44 kPa 升至21.25 kPa，增加了1.81 kPa(见图5(a))；环空干度降由0.28 降至0.26，减小了0.02(见图5(b))；当环空趾端到跟端蒸汽干度差减小时，水平段地层从趾端到跟端沿程不同位置吸收的蒸汽热焓差别也就越小，井筒沿程地层吸热率最大值与平均值的比值由1.10 降至1.05(见图5(c))，地层吸热率分布更加均匀，水平段地层均匀预热效果增强。因此，增大蒸汽干度有利于水平段地层均匀预热，且井筒沿程地层吸热率会相对增加。

3.3.3 注汽速度敏感性分析

图6所示为当注汽压力为2.0 MPa，蒸汽干度为0.95 时，不同注汽速度下井筒沿程关键参数分布规律。由图6可知：随注汽速度增加，井筒内蒸汽压力降低(见图6(a))，这是因为蒸汽质量流量增加，井筒摩擦损失增大；而蒸汽干度和地层吸热率均随注汽速度增大而上升(见图6(b)和图6(c))，这是因为注汽速度越大，热损失越少，蒸汽干度上升使得地层吸热率增大；当注汽速度由130 t/d增至170 t/d 时，蒸汽压降由12.83 kPa 升至21.24 kPa，增加了8.41 kPa(见图6(a))；环空干度降由0.61 降至0.25，减小了0.36(见图6(b))；当环空趾端到跟端蒸汽干度差减小时，水平段地层从趾端到跟端沿程不同位置吸收的蒸汽热焓差别也就越小，井筒沿程地层吸热率最大值与平均值的比值由1.11 降至1.06(见图6(c))，地层吸热率分布更加均匀，水平段地层均匀预热效果增强。因此，增大注汽速度有利于水平段地层均匀预热，且井筒沿程地层吸热率会相对增加。

3.4 注汽参数优化及预热效果分析

图5 不同蒸汽干度下井筒沿程关键参数分布Fig.5 Distribution of key parameters along wellbore with different steam qualities

上述研究表明低注汽压力、高蒸汽干度和高注汽速度有利于地层均匀预热，为分析注汽参数对SAGD 循环预热效果的影响，将未优化注汽参数(注汽压力为2.0 MPa，蒸汽干度为0.95，注汽速度为130 t/d)与优化注汽参数(注汽压力为1.5 MPa，蒸汽干度为0.95，注汽速度为170 t/d)模拟所得第3 天和第140 天的水平井筒沿程地层温度分布进行对比，如图7所示。

图6 不同注汽速度下井筒沿程关键参数分布Fig.6 Distribution of key parameters along the wellbore with different steam injection rates

为分析优化前后地层受热均匀性，引入地层温度非均匀系数αT：

图7 水平井筒沿程地层温度分布对比Fig.7 Comparison of reservoir temperature distribution along horizontal wellbore

为表征优化参数对地层预热均匀性的改善效果，定义改善系数β：

式中：αT,油田参数和αT,优化参数分别为油田在实际注汽参数条件下和注汽参数优化后对应的地层温度非均匀系数；Ti为水平井井筒沿程各微元段地层温度；-T为水平井井筒沿程平均地层温度；β为地层预热效果改善系数。

地层预热效果评价系数计算结果如表2所示，非均匀系数越小，说明地层受热越均匀；改善系数越大，说明地层预热效果改善越明显。由表2可知，优化注汽参数后，第3 天和第140 天的地层温度非均匀系数均小于未优化注汽参数时的地层温度非均匀系数，且地层预热效果改善系数分别为83.84%和27.64%。由此可知，优化注汽参数，可提高地层受热均匀程度，改善地层预热效果。