注蒸汽是稠油开采的重要并且成熟的技术手段^[1-3]。油层非均质性较强或水平段较长时，SAGD传统单管注汽方式很容易造成蒸汽局部突进和油层动用不均等情况，因此SAGD注汽井常采用同心双管注汽，即无接箍油管内和无接箍油管与内油管的环空同时注汽 (以下简称为内油管环空注汽)，以同时提高水平段跟端和趾端的加热效果^[4]。随着稠油开采技术的进一步发展，注过热蒸汽已经在一些矿场实践中取得了较好的开发效果^[5-7]。由于CO₂和N₂对提高稠油油藏采收率有很好的促进作用^[8-9]，因此常伴随蒸汽将高温高压烟道气 (主要组分为CO₂和N₂) 一同注入油层。

目前，关于同心双管注过热型多元热流体井筒热物性参数分布的研究尚处于初级阶段。国外多位学者分别提出了不同的数学模型，来计算环空中汽液两相流的压力降^[10-11]；东晓虎^[12]初步建立了海上同心双管注多元热流体数值模型，但未揭示井筒内流体相互传热特征；在此基础上，Gu Hao等人^[13]对计算内油管环空压力降的中间参数进行了简化，并对当量直径进行了改进；程文龙等人^[14]基于实际气体状态方程建立了普通单管注过热型多元热流体井筒传热数学模型，但不能直接用于同心双管注汽井筒传热的计算。为此，笔者引入R-K-S气体状态方程计算热流体密度等参数，基于质量、动量和能量守恒方程，结合经典地层内非稳态传热模型，建立了完整的同心双管注过热型多元热流体井筒传热模型，可以利用该模型进行井筒内过热型多元热流体传热特征分析、非凝结气含量优选和井底过热度预测。

1 同心双管注多元热流体井筒传热数学模型的建立

基于前人对单管/同心双管注饱和蒸汽的研究^{[8, 13-14]}，做以下基本假设：1) 井口注汽参数为定值；2) 过热型多元热流体向水泥环外缘的传热为稳态传热；3) 水泥环外缘向地层传热为非稳态传热；4) 地层热物性参数不随深度变化；5) 外油管与套管之间的内油管环空内为低压空气。

1.1 井筒内流体流动的数学模型

同心双管注过热型多元热流体井筒结构如图 1所示。

垂直井筒内流体流动过程中无质量交换，其质量守恒方程为：

(1)

式中：w_{ij, an}为无接箍油管或内油管环空中流体的质量速度，kg/s；r_{iji, an}为无接箍油管内半径或内油管内半径，m；ρ_{ij, an}为无接箍油管或内油管环空中流体的密度，kg/m³；v_{ij, an}为无接箍油管或内油管环空中流体的流速，m/s；z为井筒深度，m。

文献[4]未考虑摩擦力做功对热损失的影响，考虑无接箍油管中流体与内油管环空中流体的热传导，无接箍油管中的能量守恒方程为：

(2)

式中：Q_ij为无接箍油管向内油管环空的传热速率，W；W_ij为无接箍油管中摩擦力所做的功，W，采用文献[15]中的方法计算；w_ij为无接箍油管中流体的质量流速，kg/s；v_ij为无接箍油管中流体的流速，m/s；h_ij为无接箍油管中过热型多元热流体的热焓，J/kg；g为重力加速度，m/s²；θ为垂直井筒偏离垂直方向的角度，rad。

考虑到与无接箍油管中流体的热传导和向地层的热损失，内油管环空中的能量守恒方程表示为：

(3)

式中：Q_an为内油管环空中流体向水泥环外缘的传热速率，W；W_an为内油管环空中摩擦力所做的功，W，采用文献[15]中的方法计算；w_an为内油管环空中流体的质量流速，kg/s；v_an为内油管环空中液体的流速，m/s; h_an为内油管环空中过热型多元热流体的热焓，J/kg。

无接箍油管和内油管中过热型多元热流体的流动服从动量守恒方程：

(4)

(5)

式中：p_ij和p_an分别为无接箍油管和内油管环空中的流体压力，Pa；ρ_ij和ρ_an分别为无接箍油管和内油管环空中多元热流体的密度，kg/m³；r_iji和r_ai分别为无接箍油管和内油管环空内半径，m；r_ijo为无接箍油管外半径，m；τ_fij和τ_fan分别为无接箍油管和内油管环空中摩擦力，N，采用文献[15]中的计算方法。

式 (2)—式 (5) 即为过热型多元热流体在无接箍油管和内油管环空中流动的控制方程。

1.2 参数计算方法

式 (2) 中，Q_ij的表达式为^[4]：

(6)

其中

(7)

式中：q_ij为无接箍油管向内油管环空单位长度的导热速率，W/m；r_ijo为无接箍油管外半径，m；T_ij和T_an分别为无接箍油管和内油管环空中流体的温度，K；U_ijo为综合传热系数，W/(m²·K)；λ_tub为无接箍油管和双层同心管的导热系数，W/(m·K)；h_fiji和h_fijo分别为无接箍油管内、外管壁强迫对流换热系数，W/(m²·K)。

式 (3) 中，Q_an表示为^[4]：

(8)

其中

(9)

式中：q_an为内油管环空中流体向水泥环外壁单位深度的传热速率，W/m；r_ao为内油管外半径，m；U_ao为内油管环空中流体向水泥环外壁的综合传热系数，W/(m²·K)；T_h为水泥环外壁温度，K；r_bi和r_bo分别为外油管内壁和外壁半径，m；r_ci和r_co分别为套管内壁和外壁半径，m；r_cem为水泥环外壁半径，m；λ_ins, λ_cas和λ_cem分别为隔热层、套管和水泥环的导热系数，W/(m·K)；h_c和h_r分别为对流换热系数和辐射换热系数，W/(m²·K)。

地层内非稳态热损失速率为^[9]：

(10)

式中：T_ei为地层温度，K；f(t) 为注汽时间的函数，采用文献[16]中的方法计算。

根据微元体内质量守恒定律，由式 (8) 和式 (10) 可得：

(11)

式中：T_h为水泥环外壁温度，K。

式 (9) 中, h_c和h_r均为外油管内壁的函数^[9]，是计算q_an的关键参数，需采用迭代法求解，并在达到计算精度要求后获得q_an的精确值，步骤如下^[9]：

1) 估计U_ao初值，一般取0.5；

2) 利用式 (11) 计算T_h；

3) 利用式 (8) 计算q_an；

4) 计算T_bi和T_ci；

(12)

(13)

5) 计算h_c和h_r；

(14)

(15)

其中

(16)

(17)

(18)

式中：λ_a为环空流体导热系数，W/(m·K)；G_r为格拉晓夫数；P_r为普朗特数；σ为斯蒂芬-玻尔兹曼常数，σ=5.67×10^-8；ω_bi和ω_ci分别为隔热层外壁和套管外壁的辐射率；C_a为内油管环空内流体的比热容，J/(kg·K)；μ_a为内油管环空内流体的黏度，Pa·s。

6) 将计算得到的h_c和h_r代入式 (9)，计算U_ao^′；

7) 如果|U_ao^′－U_ao|＜ε，则输出q_an；否则返回步骤2)，重新计算。

式 (2) 和式 (3) 中混合流体的热焓h_ij和h_an为温度和压力的函数。由于实际管流为多组分混合气体，为了消除混合气体中各组分之间的相互影响，采用逸度代替分压力计算可得到精确结果^{[4, 17-18]}。

(19)

(20)

式中：f(N₂)，f(CO₂) 和f(H₂O) 分别为氮气、二氧化碳和过热蒸汽的逸度，Pa；m(N₂)，m(CO₂) 和m(H₂O) 分别为氮气、二氧化碳和过热蒸汽的质量分数；h(N₂)，h(CO₂) 和h(H₂O) 分别为氮气、二氧化碳和过热蒸汽的热焓，J/kg。

对于实际气体，高压下各组分之间的相互作用力不可忽略，此时道尔顿分压定律和基于理想气体提出的混合法则不再适用，可采用R-K-S状态方程^[19]求解式 (19) 和式 (20) 中的逸度f(N₂)，f(CO₂) 和f(H₂O)。实际气体R-K-S状态方程表示为^{[4, 20-21]}：

(21)

其中

(22)

(23)

(24)

(25)

R-K-S状态方程的逸度系数表示为：

(26)

过热型多元热流体各组分逸度的计算公式为：

(27)

式中：Z_m为过热型多元热流体压缩系数；A，B，a和b均为混合汽常数^[14]；R_m为摩尔气体常数，J/(mol·K)；y_i为过热型多元热流体中各组分的摩尔分数；a_i和b_i为过热型多元热流体中各组分汽/气体的常数^[14]；φ_i为过热型多元热流体中各组分的逸度系数。

1.3 模型求解方法

采用内迭代求函数零点，外迭代求过热型多元热流体压力和温度的方法对模型进行求解^[4]。由于求解过程中井筒内共有p_ij，p_an，T_ij和T_an等4个未知数，因此内迭代过程中需要4个独立方程才能封闭求解。分别对动量守恒方程和能量守恒方程进行差分变形得：

(28)

(29)

(30)

(31)

具体计算方法为：1) 将无接箍油管和内油管环空等分为m段，输入井口处注汽参数；2) 估计微元段出口处温度，利用式 (28) 和式 (29) 对出口端压力进行封闭求解；3) 将步骤2) 的出口端压力计算结果代入式 (28) 和式 (29)，得到一组新的出口端温度，判断估计值与计算值是否满足精度要求，若不满足，返回步骤2) 重新计算，若满足则进行下一步；4) 判断热流体是否到达油层，若未到达油层，则将微元段出口端压力、温度值作为下一微元段入口值，返回步骤2) 继续计算，若到达油层，则输出井筒内压力和温度的沿程分布。

2 实例分析

为了指导现场注汽井注汽参数的优选，首先对同心双管注多元热流体井筒传热数学模型进行了验证，分析了井筒内典型传热曲线，然后将不同注汽参数输入模型，得到不同注汽条件下井筒内热物性参数分布，并分析注汽参数对热物性参数分布的影响规律。

2.1 注汽井井筒内典型传热曲线分析

对于水平段较长的水平井，水平段跟端、趾端处的油层物性可能存在较大差异^[22]，在对无接箍油管和内油管环空分别注汽时，通过合理选择无接箍油管和内油管环空注汽参数，可以达到高效动用跟端、趾端油层的目的^[13]。W井采用同心双管注过热型多元热流体的方式开发，无接箍油管和内油管环空注汽压力、温度和日注汽量分别为4.1 MPa、650 K、205 t/d和3.6 MPa、600 K、120 t/d。无接箍油管内、外半径分别为19.0和24.2 mm，内油管内、外半径分别为38.0和44.0 mm，外油管内、外半径分别为50.9和57.2 mm，套管内、外半径分别为80.7和88.9 mm，水泥环外半径为123.6 mm。

将上述注汽参数及井筒、地层基本参数代入数学模型，计算得到W井井筒内热物性参数分布预测结果，并与实测值进行了对比，结果如图 2所示。

由图 2可知，压力和温度计算结果与实测值的最大相对误差均小于6%，满足工程计算精度要求，证明了模型的可靠性。由图 2(a)可知：1) 无接箍油管和内油管环空中多元热流体的压力均不断下降；2) 无接箍油管中的压力梯度大于内油管环空中的压力梯度，当过热型多元热流体到达油层时，无接箍油管中压力下降12.6%，内油管中压力下降3.7%。这是由于W井无接箍油管内流体质量流速比内油管环空中质量流速高72.5%，因此无接箍油管中摩擦力做功增加，压力梯度增大。由图 2(b)可知：1) 在0~25 m井段，无接箍油管中流体温度迅速下降，在井深50 m后温度梯度趋于稳定，并与内油管环空中温度梯度相近；2) 在0~25 m井段，内油管环空中过热型多元热流体温度有所升高，在井深25 m后温度开始下降，在井深50 m后，其温度梯度与无接箍油管中温度梯度相近。分析认为，在0~25 m井段，无接箍油管和内油管环空中温差较大，由于无接箍油管为金属材料，导热性能良好 (模型中U_ijo值为1 450 W/(m²·K))，其流体热损失速率很大，导致流体温度剧烈变化。

2.2 非凝结气含量分析

在实际生产过程中，过热型多元热流体中的非凝结气是由柴油和空气以1.0:14.9的质量比经过充分燃烧得到的，按照元素质量分数关系，非凝结气中N₂和CO₂质量比约为4:1^[22-33]。因此，在其他注汽参数不变的条件下，设定模型中氮气与二氧化碳的质量分数比分别为4%:1%，20%:5%和40%:10%，基于R-K-S方程的模型计算井筒内的热物性参数分布，结果如图 3所示。

由图 3(a)、图 3(b)可知，随着过热型多元热流体中非凝结气含量增加，无接箍油管和内油管环空中流体温度均下降; 随着非凝结气质量分数增大，无接箍油管中流体向内油管环空传热速率减小，内油管中流体净热损失速率增大，温度降低。由图 3(c)、图 3(d)可知，随着非凝结气质量分数增大，井筒内流体过热度均下降。实际上，过热型多元热流体中的非凝结气主要用于改善油藏内渗流特征，而非用于携带热能。因此，多元热流体中过热蒸汽比例过低会降低油层加热效果，应结合实际油藏与非凝结气的作用特征，对过热蒸汽和非凝结气的质量比例进行优选。

2.3 注汽温度优化

在井口注汽压力和注汽速度不变的条件下，可以利用该模型对注汽温度进行优化。以无接箍油管注汽温度为例，取注汽温度分别为550，570，590，610，630和650 K，基于R-K-S方程的模型计算井筒内热物性参数分布，结果如图 4所示。

由图 4(a)可知，随着注汽温度升高，无接箍油管中压力梯度增大。这是因为温度升高，流体密度减小，流速增大，摩擦力增大，压力梯度增大。

由图 4(b)可知，无接箍油管注汽温度对内油管环空压力梯度影响不大。

由图 4(c)、图 4(d)可知：1) 当无接箍油管注汽温度低于590 K时，以无接箍油管注汽温度550 K为例，在0~40 m井段无接箍油管中过热度迅速升高，此时无接箍油管中流体从内油管环空中吸收热能，内油管环空中过热度下降；2) 当无接箍油管注汽温度高于610 K时，在0~40 m井段，热量由无接箍油管传向内油管环空，而内油管环空开始吸收热能，此时内油管环空中过热度增加；3) 在井深约40 m处至井底，无接箍油管中过热度不断下降，此时无接箍油管与内油管环空间的相互传热量很小，对过热度影响很小。需要指出的是，受锅炉最高输出蒸汽温度的限制，现场应结合实际情况合理选择注汽温度，以获得较高的井底过热度。

3 结论

1) 模型计算结果表明，在近井口处，无接箍油管和内油管环空之间的温差较小，导致流体热物性参数剧烈变化，但温度梯度快速趋于一致；非凝结气含量增大，井底过热度降低；随着无接箍油管注汽温度升高，无接箍油管和内油管环空井底过热度均增加。

2) 现场应用时应结合实际油藏与非凝结气的作用特征对过热蒸汽和非凝结气质量比例进行优选；结合蒸汽发生装置性能，适当提高注汽温度，以获得较高的井底过热度。

3) 需要进一步研究油管与套管发生接触时的井筒内热参数分布特征和注汽速度随时间变化时的井筒内热参数分布特征等问题。