引言稠油又称重油，是21世纪一种重要的非常规油气资源，特点是沥青、胶质、石蜡等杂质含量高，因而黏度高、流动性差，导致开采难度大、采收率低。提高稠油采收率并降低生产成本，对保障国家能源战略安全有重要意义。目前，国内稠油开采技术可分为热采和冷采[1]，其中以蒸汽吞吐应用最广，并辅以电加热、循环热流体等井筒伴热方法。重力热管伴热是一项自平衡的“无源”技术，相较其他伴热方法具有清洁高效、能耗低、投资少、运行可靠等优点，最早由李菊香[2]等提出，经过多年的理论研究和试验论证，已被证明为一种实用的稠油开采技术。重力热管应用于稠油热采系统，主要采用重力热管式抽油杆和重力热管式油管两种方式[3-4]。重力热管采油系统的研究多集中在理论简化模型的求解、现场或室内试验，还未开展完整的数值研究，而数值模拟节约成本和时间，便于进行系统的工况匹配与性能优化。稠油井井深约1 000 m，实际应用中多采用串级式重力热管，其单支长度一般在10 m以上，属于超长重力热管，因而引用Wang[5]等提出的压力基相变模型作为重力热管内部的相变模拟方法。鉴于目前的研究现状，选取串级重力热管式稠油热采系统对其进行理论建模和数值研究，分析重力热管的流动传热特性及伴热效果，并对节能效果与应用前景进行评估。1简化传热模型理论分析本文以辽河油田齐108—20—26井的井筒参数为基础进行几何建模，泵挂深度780 m，系统几何参数如表1所示。选取测试的某稠油样品，热物性参数拟合公式如表2所示，典型数值如表3所示。动力黏度为90.30 ℃下的计算值，密度为79.31 ℃下的计算值，导热系数为90.10 ℃下的计算值，比定压热容为30.71 ℃下的计算值。10.3969/j.issn.1004-7948.2021.10.001.T001表1串级重力热管式油管系统几何参数井筒参数数值抽油杆外径d10.020油管内径d20.062外径d30.073热管外管内径d40.098外径d50.114套管内径d60.126外径d70.140水泥环外径d80.240m10.3969/j.issn.1004-7948.2021.10.001.T002表2稠油热物性参数拟合公式热物性参数拟合公式动力黏度μμ=7.238 5×10-12e7 686.520 8T-56.975 3密度ρρ=1 181.1-0.726 7T+1.268 6×10-5T2导热系数λλ=0.330 4-2.898 0×10-4T+2.725 6×10-7T2比定压热容cp—10.3969/j.issn.1004-7948.2021.10.001.T003表3稠油热物性参数典型数值热物性参数数值动力黏度μ/（Pa·s）0.668密度ρ/（kg/m3）927.6导热系数λ/［W/（m·K）］0.261比定压热容cp/［J/（kg·K）］2 454.61.1简化假设（1）稠油井筒传热视为径向一维稳态传热。（2）稠油为连续单相、均质、常物性的牛顿流体，流动视为一维稳态层流，忽略自身轴向导热。（3）忽略抽油杆上下运动（摩擦）引起的温度变化。（4）每段重力热管均正常工作，蒸汽处于饱和状态，忽略热管内部轴向热阻（温差）和连接段长度。（5）热管外管与套管之间的环形空间充满常物性的水，将自然对流折算为纯导热，且忽略热辐射。（6）因套管外的水泥环与地层导热系数相近，将其视为地层导热的一部分[6]。（7）井筒向地层导热为各向同性非稳态导热，物理模型简化为外部边界固定的圆柱体，地温梯度、油层温度、地面环境温度均为常数[7]。1.2热平衡方程单级重力热管式油管的传热模型如图1所示。10.3969/j.issn.1004-7948.2021.10.001.F001图1单级重力热管式油管系统传热模型由此列出能量平衡方程如下：Hin=Hself+Q1 (1)Hself=Q2+Hout (2)Q1+Q2=Qe+Qc (3)dΦ1=-qmcpdtf=k1tf-tvπd2dz (4)dΦ2=-qmcpdtf=k2tf-tvπd2dz (5)式中：Φ1——热管蒸发段单位时间内传热量，W；Φ2——冷凝段单位时间内传热量，W；qm——稠油质量流量，根据油井的日产液量换算得出，kg/s；tf——稠油温度，℃；tv——热管工质的饱和温度，℃；k1——以油管内表面积为基准的热管蒸发段，W/m2⋅K；k2——冷凝段总传热系数，W/m2⋅K。边界条件为：z=0,tf=te' (6)z=Le,tf=te″=tc' (7)z=Le+Lc,tf=tc″ (8)对无热管系统进行能量平衡分析：dΦ3=-qmcpdtf=k3tf-trπd2dz (9)式中：Φ3——单位时间内地层散热损失，W；tr——地层温度，℃。地层温度计算公式如下：tr=td-gTz (10)式中：td——井底原始地层温度，℃；gT——地温梯度，℃/m，通常在0.030~0.035 ℃/m之间。边界条件为：z=0,tf=tf' (11)z=Le+Lc,tf=tf″ (12)式中：tf'——稠油在油管进口的温度，℃；tf″——稠油在油管出口的温度，℃。1.3热阻（集总参数）模型取热管蒸发段井筒某横截面，径向温度分布如图2所示。运用集总参数法建立系统热阻模型，如表4所示。10.3969/j.issn.1004-7948.2021.10.001.F002图2重力热管蒸发段井筒某横截面径向温度分布10.3969/j.issn.1004-7948.2021.10.001.T004表4系统热阻模型热阻表达式R11πh1d2R212πλ1lnd3d2R31πh2,ed3或1πh2,cd3R41πh3,ed4或1πh3,cd4R512πλ2lnd5d4R612πλ3lnd6d5R712πλ4lnd7d6R812πλ5ln4aτd9-0.29注：e——重力热管蒸发段；c——重力热管冷凝段。［K⋅m/W］由传热各环节热阻得出相应的总传热系数：k1=πd2∑i=13Ri,e-1 (13)k2=πd2∑i=13Ri,c-1 (14)k3=πd2R1+R2+R6+R7+R8-1 (15)2CFD数值模拟2.1几何模型与网格划分受计算资源和时间所限，重力热管选取60支，每段长度13 m。同时根据简化假设，每段重力热管的工况相似，某段稠油的出口温度作为下段入口温度，因此本模拟选取井筒底端第一段作为模拟对象。为节省计算量，只计算流体区域，并绘制无热管系统比较。利用ICEM CFD软件进行网格划分，如图3所示。使用结构化网格并在壁面附近设置边界层，网格总数分别为721 603和331 547。10.3969/j.issn.1004-7948.2021.10.001.F003图3几何模型与网格划分2.2数值模型多相流模型选取VOF方法，为兼顾稳定性和收敛速度，将库朗数设置为0.5，体积力选取隐式算法，黏性模型为层流，蒸发—冷凝相变模型选取一种基于压力的相变模型，表达式如表5所示。本模型控制方程组如下：∂αqρq∂t+∂αqρquj∂xj=SM (16)∂∂tρui+∂∂xjρuiuj=ρfi-∂p∂xi+∂∂xj2μsij-2∂3∂xiμskk+SCSF (17)∂∂tρε+∂∂xjρεuj=λ∂2T∂xj2+ρfjuj-p∂uj∂xj+∂∂xj2μsijui-2∂3∂xjμskkuj+SE (18)式中：αq——气相或液相的体积分数；SM——质量源项，kg/m3⋅s；SCSF——表面张力转化的体积力，N/m3；SE——能量源项，J/m3⋅s，等于SM与气化潜热的乘积。表5中psatTmix为混合物温度对应的饱和压力，其拟合公式可由Clausius-Clapeyron方程导出。△pi为形成气化核心所需压差，计算公式如下：△pi=2σLVri1+ρVρL (19)式中：ri——气化核心的半径，μm，一般在1~10 μm之间。上述质量源项、能量源项通过UDF添加到控制方程中。10.3969/j.issn.1004-7948.2021.10.001.T005表5压力基相变模型表达式相变过程判据相质量源项蒸发psatTmix≥p+△pi气相SM,V=βeαLρVpsatTmix-pp液相SM,L=-βeαLρVpsatTmix-ppppsatTmixp+△pi气相SM,V=βeαLρVpsatTmix-pp (0αL1)0 (αL=1)液相SM,L=-βeαLρVpsatTmix-pp (0αL1)0 (αL=1)冷凝psatTmix≤p气相SM,V=-βcαVρVp-psatTmixp液相SM,L=βcαVρVp-psatTmixp注：下标V——气相；下标L——液相；βe——蒸发时间松弛因子，s-1；βc——冷凝时间松弛因子，s-1。重力热管工质和壳体材料分别选取水和内表面经钝化处理的不锈钢，套管材料为普通不锈钢，导热系数均为16.27 W/(m∙K)；地层热阻R9取0.3 (m∙K)/W；环空介质水的当量导热系数取为0.648 W/(m∙K)。定义液相水为第一相，气相水为第二相，稠油为第三相。其中液相水的密度拟合公式和水的气液两相表面张力系数拟合公式分别如下：ρ=919.8+0.923 7T-0.002 2T2 (20)σLV=0.098 058 56-1.845×10-5T+2.3×10-7T2 (21)气相水热物性参数为Fluent默认值。稠油热物性参数如表3所示，其中动力黏度通过UDF设置。本系统边界条件设置如图4所示，显示比例y∶x=1∶10。热管外管设置为第三类边界条件，热管冷凝段假定为Nusselt层流膜状凝结，蒸发段液膜区和液池区均假定为核态沸腾，其表面传热系数通过相应经验或半经验公式计算，均取为6 000 W/m2⋅K，相应总传热系数为9.91 W/m2⋅K；地层温度通过UDF设置，井底原始温度为43.4 ℃，地温梯度0.03 ℃/m。10.3969/j.issn.1004-7948.2021.10.001.F004图4本系统边界条件设置为使蒸发段处于核态沸腾，将水的饱和温度设置70 ℃，汽化潜热相应为2 333 kJ/kg。将前述的压力基相变模型转换为UDF，其中ri为2 μm，βe为0.1，βc为0.5。操作压力设置为饱和压力。为减少计算时间，初始温度设置为饱和温度，初始表压设置为0，重力热管的充液率为20%，即初始液相体积占重力热管总体积的20%。采用SIMPLE算法，基于Fluent软件求以上控制方程，对于本系统，瞬态计算时间步长为0.000 5 s，达到准稳态的时间约30.0 s，对于无热管系统，稳态计算步数为2 000步。2.3结果讨论2.3.1相分布CFD模拟不同时刻的液相水体积分数云图如图5~图7所示。由图5可知，t=0.5 s时，随着稠油自油管入口向上流动，重力热管蒸发段底部被加热形成汽化核心，从而产生气泡。由图6可知，t=2.0~10.0 s时，气泡不断上浮长大并伴随着破碎与合并，液池内部流型经历泡状流、柱塞流、弹状流、块状流的转变，同时新气泡也在不断产生。由图7可知，t=30.0 s时，壁面温度变化趋于平缓，模拟达到准稳态，即重力热管完成启动。图7(a)显示蒸发段呈现核态沸腾；图7(b)显示气液交界面由于气泡浮出而出现波动；图7(c)显示冷凝段中段有较多下落的液滴，同时壁面附有不连续的液膜。10.3969/j.issn.1004-7948.2021.10.001.F005图5t=0.5 s液相水体积分数云图10.3969/j.issn.1004-7948.2021.10.001.F006图6t=2.0~10.0 s液相水体积分数云图10.3969/j.issn.1004-7948.2021.10.001.F007图7t=30.0 s液相水体积分数云图冷凝段中由于气液界面的剪切应力作用使液滴被夹带[8]，且液膜不断受到稠油向外散热，从而蒸发脱离壁面形成液滴，不能连续地流回液池。2.3.2温度分布CFD模拟达到准稳态时的温度云图如图8所示。随着稠油自入口向上流动并散热，截面平均温度不断下降，壁面温度边界层不断发展增厚。10.3969/j.issn.1004-7948.2021.10.001.F008图8准稳态温度云图同时，重力热管蒸发段除底部有较大过热度外基本处于局部饱和温度，而冷凝段内部的温度分布由于下降的液滴与液膜的存在极不均匀，并且其外壁面温度更接近相应地层温度，因此表面传热系数高于蒸发段。从准稳态的体积分数分布和温度分布可以得出，本系统重力热管蒸发段传热方式以核态沸腾为主，冷凝段传热方式包含自然对流、膜状凝结和珠状凝结，而绝热段则几乎不存在。2.3.3速度分布CFD模拟达到准稳态时的速度矢量图如图9所示。图9(a)显示稠油入口段速度边界层不断发展以及蒸发段液池底部气泡不规则的流动；图9(b)显示冷凝段下段脱离壁面的液滴下降过程速度不断增大，并与上升气流进行强烈的动量和能量交换，形成较大的涡旋气流；图9(c)显示冷凝段中段的速度分布较为有序，上升的蒸汽气流受到冷却而回流，但两侧重力热管并不对称分布，由于操作压力点的设置和速度的随机脉动等多方面因素影响。综合以上体积分数、温度和速度分布，表明本系统重力热管的设计参数较为合理，运行较为稳定，而其流动、传热特性与传统重力热管有较大差异。10.3969/j.issn.1004-7948.2021.10.001.F009图9准稳态速度矢量图2.3.4稠油轴向温度分布重力热管式油管系统与无热管系统的稠油轴向温度分布如图10所示，理论计算值通过一维传热模型计算得到，CFD模拟值通过各截面质量加权平均得到。图10中1为重力热管式油管系统，2为无热管系统。10.3969/j.issn.1004-7948.2021.10.001.F010图10稠油轴向温度分布由图10可知，显示单级长度内稠油的轴向温度近似线性分布，串级式重力热管式油管第一段稠油出口温度的理论计算值、CFD模拟值分别为88.82 ℃、89.04 ℃。无热管系统相应的稠油温度的理论计算值、CFD模拟值分别为83.37 ℃、83.98 ℃，即重力热管的应用使得单级油管的稠油温降减缓5.45 ℃（理论计算值）/5.06 ℃（CFD模拟值）。以上结果表明，重力热管与井筒的耦合可减缓稠油举升过程中的温降，从而提高流动性、降低抽油机负荷。此外，流动与散热时间的减少提高开采效率，重力热管的取热作用减少高品位地热能损失，实现热量自平衡采油。文中提出重力热管式油管系统的CFD数值模拟方法对超长重力热管的设计选型与工质筛选、系统推广应用与性能优化等具有一定的指导意义。3结语针对串级式重力热管耦合井筒的稠油热采系统进行理论建模和数值研究，得出以下结论：（1）在合理简化假设基础上，建立本系统一维传热模型，得出稠油轴向温度分布和地层散热损失的表达式。（2）利用Fluent进行CFD数值模拟，得到不同时刻的体积分数云图、准稳态时的温度云图和速度矢量图，并对重力热管的传热特性和流动特性进行分析。（3）将本系统的稠油轴向温度分布与无热管系统比较，结果表明重力热管的伴热效果良好，本系统运行工况合理，具有较好的应用前景。