不同管径下超临界CO2流动传热的数值分析

李超， 王宗一， 廖敏， 赵钦新， 邵怀爽

(西安交通大学 热流科学与工程教育部重点实验室，陕西 西安 710049)

超临界二氧化碳sCO2布雷顿循环在核电领域具有诸多优势，循环效率能够比蒸汽工质高12%左右，安全性好，并且设备更加紧凑，具有良好的应用前景[1]。然而超临界二氧化碳物性随温度剧烈变化，其流动换热特性与常物性流体迥然不同，充分认识并掌握sCO2流动传热特性对于其在第4代核电站的实际应用具有重要意义。目前针对sCO2的换热特性已有大量的实验和模拟研究，一些学者通过实验研究了加热条件下sCO2的传热恶化现象[2-3]。Wang等[4]通过实验研究了加热条件下sCO2在垂直管中的对流换热特性，并与水平管进行比较，分析了压力、热流密度、质量流速、入口温度以及管径对传热的影响，研究结果表明浮升力会强化垂直向下流动的传热，恶化垂直向上流动的传热。Kim等[5]通过实验分析了sCO2在向上流动的管道中的传热特性，热流密度大于某一值时管道会出现传热恶化现象，且环形通道能够减弱恶化程度。Liu等[6]通过实验研究了浮升力对sCO2在竖直加热管道内流动换热的影响，实验结果表面浮升力在小于某一值时传热恶化消失，且传热恶化随热流密度的增大而向进口处移动。对于水平管，Liao等[7]通过实验研究得出在微细管中水平与垂直向下流动的传热得到强化。Zhang等[8]同样通过实验对不同流动方向的4 mm加热管进行了研究，通过管道壁温变化实验数据说明了水平方向与垂直向上流动出现了传热恶化现象，但水平方向的恶化程度弱于垂直向上流动。

另外，随着计算流体力学的快速发展，数值模拟方法作为一种更加经济的研究方法被广泛用于超临界流体传热特性的研究中[9]。Cai等[10]通过数值模拟研究，表明sCO2在浮升力作用下产生的二次流可以增强管道顶部流体的传热并削弱管底部的传热。闫晨帅等[11]通过数值模拟的方法，解释了加热工况下顶母线和底母线传热出现明显差异的机理，研究了水平管中“类液-类气分层”现象对传热恶化的影响。Zhao等[12]对大直径水平加热管(22.14 mm)内超临界二氧化碳传热恶化现象进行了数值模拟研究，并修正了更加适用于大直径管的浮升力准则数。此外，根据场协同原理[13]，对流换热的增强不仅取决于物理性质和流型，还取决于流体的速度矢量与温度梯度之间的夹角。一些研究者利用场协同原理分析了sCO2管内流动换热[14-15]。Xiang等[16]对水平圆管中sCO2的冷却传热性能进行了数值模拟研究，并利用场协同原理分析了截面上传热性能的差异。Yang等[17]分析了sCO2在带有肋片的竖直管道内的传热，证明了场协同原理可以成功地预测传热恶化或增强的程度。

综合以上研究可知，sCO2管内流动换热规律复杂，现有的实验与数值模拟研究并未完全阐明水平管中传热恶化的物理机制。本文通过数值模拟的方法研究了加热工况下管径对sCO2在水平管中的流动换热特性影响，并结合场协同理论，通过分析管内流场、浮升力参数以及努塞尔数的变化，揭示了管径对于浮升力效应的影响，并量化了新涡流对于顶部的传热恶化程度。

1 水平换热管物理模型与数值方法

1.1 几何模型

本文的管道几何模型如图1所示，包括入口段和加热段2部分，为了减小入口效应对流动换热的影响，设有200 mm长的绝热段作为入口段；加热段为均匀热流密度条件，长度为500 mm；为分析管径对超临界二氧化碳流动换热的影响，分别选用管径为6、8和10 mm的管道。

图1 水平换热管物理模型Fig.1 Physical model of horizontal heat exchange tube

1.2 控制方程

稳态湍流对流换热的控制方程组在直角坐标系中为[18]：

连续性方程：

(1)

动量方程：

(2)

能量方程：

(3)

式中μeff是有效粘性系数，其计算公式为：

(4)

式中μt为湍流动力粘度，其计算公式为：

(5)

(6)

(7)

式中：ρ为密度，kg/m3；u为速度，m/s；T为温度，K；p为压力，Pa；μ为动力粘度，Pa·s；y为到最近壁面的距离；a1为常数，值为0.31；k为湍动能，m2/s2。

本文基于Fluent商业CFD软件，选择SSTk-ω湍流模型作为数值计算中的湍流模型[19]。该模型通过使用混合函数，在近壁面处采用k-ω模型，其他区域采用k-ε模型。与标准k-ω模型相比，SSTk-ω湍流模型中使用了不同的湍流常数项，引入了横向耗散导数项，并且在湍流粘度的定义中考虑了湍流剪应力的输运过程，使得SSTk-ω模型适用性更高。Wang等[20]通过对比SSTk-ω模型与其他湍流模型得出SSTk-ω湍流模型能够对传热系数、壁面温度和主流温度进行更好的预测，除此之外，也有许多研究采用了该湍流模型进行研究并得到了准确的结果[16,21-22]。对于SSTk-ω湍流模型，湍动能k的运输方程为：

(8)

比耗散率ω的运输方程为：

Gω-Yω+Dω+Sω

(9)

式中：Gk和Gω为湍动能和比耗散率的产生项；Гk和Гω为湍动能和比耗散率的有效扩散系数；Yk和Yω是湍动能和比耗散率的耗散项；Dω为交叉扩散项，Sk和Sω为湍动能k和比耗散率ω的自定义源项。

1.3 网格划分及计算方法

本文通过ANSYS ICEM软件划分模型网格，截面通过O型切分生成结构化六面体网格，由于超临界二氧化碳的物理性质随温度发生剧烈变化，因此对近壁面处的网格进行局部加密，由于使用了SSTk-ω湍流模型，因此需保持第1层网格y+<1，网格划分如图2所示。为确定网格数量对于计算结果的影响，需进行网格无关性验证。根据网格疏密程度不同，生成了5套网格，网格数量分别为659 610、1 074 920、1 790 370、3 015 360、4 149 610。不同网格下管道总压降与流体出口温度的计算结果如图3所示，管道压降随网格数的增加而减小，出口温度随网格数的增加而增加。从图中数据走势可以看出，当网格数量增大到一定程度后，压降与出口温度随网格数量的进一步增加基本趋于稳定。综合考虑计算资源和计算精度，本文选用第3套网格生成方案执行相应的计算任务。

图2 水平管网格划分Fig.2 Meshing diagram of horizontal tube

考虑到超临界二氧化碳物性变化剧烈，选用FLUENT中的NIST 实际气体模型以反映流体物性的剧烈变化给传热带来的影响。模拟入口设置为质量流速入口，出口为压力出口，入口段设置为绝热条件，加热段为均匀热流密度加热。数值模拟计算采用压力基、稳态求解器进行求解，求解方法采用基于压力-速度耦合的COUPLED算法。压力标准差分、连续性方程、动量方程、能量方程、湍动能以及比耗散率项均采用二阶迎风格式进行离散。

图3 水平管道模型网格独立性验证Fig.3 Mesh independence verification of horizontal tube model

2 数值模拟计算结果及分析

2.1 数据处理方法

超临界二氧化碳流体在管内流动的局部对流换热系数h为：

(10)

(11)

格拉晓夫数Gr和主流区雷诺数Reb分别表示为：

(12)

(13)

努塞尔数表示为：

(14)

式中：q为壁面热流密度，W/m2;T为温度，K；A为管道换热面积，m2；ρ为密度，kg/m3;u为速度，m/s；d为管道直径，m；μ为动力粘度，Pa·s；λ为导热系数，W/(m·K)；下标b表示主流区参数，下标w表示壁面，下标in表示流体入口，下标out表示流体出口；Tw为加热壁面平均温度。

2.2 可靠性验证

为验证数值模型的准确性，将数值模拟结果与Dang和Hihara的实验结果[23-24]进行对比实验工况进口温度为350 K, 质量流速为200 kg/(m2·s)，壁面热流密度为12 kW/m2，压力为8 MPa，实验测量结果与数值模拟计算结果如图4所示。可以看出，模拟结果和实验结果的吻合较好，具有相同的变化趋势，最大误差发生在跨临界区，误差为18.84%，最小误差为0.23%，平均误差为7.49%，说明本文采用的数值方法准确可行。

图4 数值方法验证Fig.4 Numerical method validation

2.3 管径对换热效果的影响

为研究加热工况下管道直径大小对sCO2管内流动换热特性的影响，分别对3种不同管径的管道进行数值模拟分析。管道直径取为6、8和10 mm，入口质量流量为0.006 kg/s，加热段的壁面热流密度为50 kW/m2，流体入口压力为7.5 MPa，入口温度为310 K。表1为3种管径对应的模拟结果，在相同的入口质量流率下，管径越大，管内流速越小，而对流传热系数受到管内流速的影响，因此3种管径的平均传热系数相差较大。为了研究管径对传热效果的影响，本文通过改变重力加速度的大小来分析不同管径中浮升力在加热情况下对整体传热性能的影响。浮升力是由于截面内密度分布不均匀，在重力的影响下，密度较大的流体向下运动，而密度较小的流体则向上运动，在截面内形成二次环流，在重力加速度为0时截面内则不会出现二次环流现象。从表中数据可以看出，重力加速度为0.5g时，管径为6、8和10 mm的管内平均传热系数比0g时分别增加了3.45%、12.86%和23.52%；重力加速度为g时，管径为6、8和10 mm的管内平均传热系数比0g时分别增加了4.80%、15.73%和25.84%。模拟结果表明，在本文选定实验工况范围内，浮升力效应会强化管道整体的换热效果；并且在等质量流率条件下，管径增大会让浮升力对于整个管道的流动换热特性影响增大；同时也可以看出重力加速度增大后，浮升力对于整个管道的流动换热特性影响也随之增大。

表1 3种管径的传热效果数值模拟结果

2.4 管径对流场的影响

图5为不同管径对应的密度分布云图和比热容分布云图。受重力影响，sCO2在水平管道内流动时呈现出上下分层的流动形态，顶部流体密度较小，底部流体密度较大，其流体物性分别与气态和液态sCO2相近，也可称为“类气”和“类液”sCO2[11]。由比热容云图分布可以看出，管道顶部“类气”流态的sCO2不仅密度较小，其比热容同样较低，随管径增大管道流通截面积增加，相同质量流量下sCO2工质的平均流速降低，管道顶部积聚的低密度“类气”流体相应增多，超临界流体上下分层愈加明显，其聚集在管道顶部会使得壁面局部传热恶化现象更加严重。

图5 不同管径密度和比热容分布Fig.5 Density and specific heat capacity distribution of different pipe diameters

2.5 浮升力对传热性能的影响

(15)

图6 浮升力参数沿流向变化Fig.6 Variation of buoyancy parameters along the flow direction

为研究不同管道直径下浮升力效应的影响程度，进一步分析了壁面局部传热系数随流动方向的变化，如图7所示。有重力时的平均传热系数和底母线上局部传热系数均大于无重力时的传热系数，但顶母线的传热系数小于无重力时的传热系数，说明在加热情况下传热强化与传热恶化同时存在，在重力作用下产生的浮升力，强化了整个管道和底母线的传热效果，但同时也恶化了管道壁面顶母线附近的传热；通过对比不同管径的传热情况可发现，随着管径的增大，顶母线和底母线的传热系数相差越来越大，表明管径增加引起的浮升力增强，会同时加剧传热强化和传热恶化这2种现象。

图7 不同管径表面传热系数沿流向变化Fig.7 Variation of heat transfer coefficient of different pipe diameters along the flow direction

为分析管径对管道顶部传热恶化的影响，本文取Z=375 mm截面处的流线及密度分布云图进行分析，如图8所示。随着浮升力的增大，截面内流体的密度分布更加不均匀，并且在管径为8和10 mm时，管道顶部出现了新的涡流，而这部分低密度流体导热能力较差且不参与主流区的流动，聚集在管道顶部使得传热恶化进一步加剧。

图9为顶母线Nutop和周向平均Nuave与无重力工况下Nuzero的比值。在加热时，顶母线的传热由于浮升力的出现而发生恶化，而无重力工况时没有浮升力，本文用努塞尔数的比值来评估顶母线附近流体传热的恶化程度，该比值越小，顶母线的传热能力越差，说明传热恶化越剧烈。由顶母线Nu比值可以看出，管径为8和10 mm的比值要远小于管径为6 mm时的Nu比值，这与图8中管道顶部新涡流的出现有关。管径为6 mm时，Nu的平均值为0.55，管径为8和10 mm时平均值分别为0.35和0.33，说明了新涡流的出现会使管道顶部传热恶化进一步加剧，Nu的比值降低了约0.2。从周向平均Nu比值可以看出，管径越大，浮升力对管内传热影响越大，这与上文所得结论一致，在加热时管内传热恶化和传热强化现象同时存在，管道整体传热效果得到强化，而顶部则会发生剧烈的传热恶化现象。

图8 不同管径下Z=375 mm截面处的流线及密度分布云图Fig.8 Streamline diagram and density distribution at Z=375 mm under different pipe diameters

2.6 场协同分析

场协同原理[13]从流体速度场和温度场相互作用的角度研究了对流传热机理，认为对流传热的强化或弱化不仅和物性、流动状态有关，还与速度矢量和温度梯度场的协同性有关。局部场协同角θ为：

(16)

式中：V为流体的速度矢量；▽T为流体的温度梯度。场协同角表示速度场和温度场之间的协同程度，场协同角为0°时，场协同程度最好。图10为不同管径对应的局部场协同角分布图。

图9 顶母线、周向平均与无重力工况的Nu比值Fig.9 The Nu ratio of top bus to zero gravity and circumferential average to zero gravity

图10 局部场协同角分布Fig.10 Local field synergy angle distribution

本文沿流向分别取Z=125,250,375 mm 3个截面进行分析，具有较小场协同角的流体主要聚集在管道的底部，说明管道底部的传热得到了强化，同时结合图7、8进行分析，超临界二氧化碳在加热时传热强化与恶化现象同时出现在管道底部和顶部；随着管径的增加，场协同角分布变得更不均匀，具有较小场协同角的流体更加靠近壁面；并且可以看出在管径增大后，管道下游的局部场协同角增大，接近90°，即速度场和温度梯度的夹角趋近于90°，这是由于管径增大后，x方向与y方向的速度分量减小，即截面内主涡流的强度减小，因此在管道底部流体的局部场协同角随之增大。

3 结论

1) 浮升力效应会强化管道整体的换热效果，并且在等质量流率条件下，管径越大，浮升力对于整个管道的流动换热特性影响越大。

2) 随着管径的增大，管道顶部积聚的低密度“类气”流体相应增多，超临界流体上下分层愈加明显，顶母线和底母线的传热系数相差越来越大，表明管径增加引起的浮升力增强，会同时加剧传热强化和传热恶化这2种现象。并且管径较大时，“类气”流体会在管道顶部形成新的涡流，本文通过努塞尔数的比值得出新涡流的出现会使管道顶部传热恶化进一步加剧。

3) 通过场协同原理分析了管内速度场与温度场的协同程度，管径增大后会使主涡流强度减小，并且场协同角分布更加不均匀，具有较小场协同较的流体更加靠近壁面。