Chinese Journal of Ship Research
竖直上升管内超临界CO2异常传热机理研究
胡娟娟*，李增光，高一民，杨建，刘浪中国舰船研究设计中心上海分部，上海 201108
摘 要:［目的］小通道内的换热已被广泛应用于核反应堆与航天发动机冷却等方面。为保证该系统内热力装置安全、高效运行，需系统、深入研究超临界CO2的传热规律与异常传热机理。［方法］通过CFD 数值模拟，对不同质量流速下内径为2 mm的竖直上升管中超临界CO2的换热随热流密度与压力的变化规律展开研究，并对传热强化与传热恶化工况下不同截面处径向物性参数的变化规律进行研究，深入分析异常传热机理。［结果］研究结果表明，当质量流速为500~1 000 kg/(m2·s) 时，随着热流密度的升高，超临界CO2 由传热强化转变为传热恶化，且在传热恶化工况下，随着热流密度的升高，超临界CO2传热恢复区的换热系数峰值有所降低，直至消失；随着运行压力的升高，超临界CO2的异常传热程度均被削弱。通过机理分析发现，在传热强化工况下，超临界CO2的换热在拟临界点附近达到最强， 其主要原因在于管截面上大比热容区流体的份额显著增加；在传热恶化工况下，壁温峰值对应的截面处速度呈“M”型分布，拐点处速度梯度降为0，湍动能也降到最低，导致传热能力变弱。［结论］研究结果对于改善小通道内超临界CO2的传热性能具有一定的指导意义。关键词：超临界CO2；竖直上升管；异常传热；机理中图分类号: U664.15文献标志码:A DOI：10.19693/j.issn.1673-3185.01921
0 引 言
超临界 CO2 （S-CO2 ）布雷顿循环因其优异的性能，在船舶余热利用、舰艇动力推进、能源发电、分布式能源等领域均已获得应用[1- 2]，并展现出广阔的应用前景。与水蒸汽相比，超临界CO2有诸多优点：临界条件易实现、稳定性好；临界点附近具有高密度，可降低压缩功；功率密度高，有利于减小循环系统的尺寸[3-4]。循环工质的换热性能与系统稳健运行密切相关。超临界压力下，CO2的物理性质十分特殊：密度与液体相近，动力黏度与气体相近，扩散系数介于气、液之间[5]。此外，在拟临界点附近，超临界 CO2的定压比热容的变化非常敏感，例如当压力为 7.5 MPa时，比热容峰值附近的温度变化0.2 ℃时对应的比热容变化可高达460%。物性的特殊性使超临界CO2的流动换热规律十分特殊，而在工程实际中，超临界机组工况随运行参数的变化极其复杂。此外，随着科技的进步，热交换器与制冷器正在不断向微小化、紧凑化发展，在核反应堆、空调与热泵以及航天卫星等领域均有较大的应用前景。因此，有必要对小通道内超临界CO2在加热管内的流动和换热规律进行深入研究，并对传热强化、传热恶化机制进行探究，用以从根本上减轻传热恶化带来的危害，保障超临界机组稳健运行。早在 20 世纪 50 年代，就有学者研究了超临界 CO2。Kim 等[6-8] 对超临界CO2 在内径 d=4.5 mm上升管与下降管内的对流换热特性进行了研究，结果显示热流密度和质量流速是影响加热管管壁温度分布的主要因素，并基于实验数据提出了换热关联式；Zahlan 等[9] 在亚临界、近临界和超临界压力下，分别对内径 d=8，22 mm竖直上升管中超临界 CO2的流动与传热特性进行了实验研究，建立了超临界CO2流动换热数据库，并验证了其有效性；张丽娜等[10] 对 d=4 mm水平圆管内超临界CO2的换热特性进行了数值模拟，研究了热流密度、质量流速、压力等参数的影响。由此可见，现有研究多针对较大通道内超临界CO2 的传热规律，对小通道内超临界CO2异常传热机理的研究不足。同时有研究表明，管径越小，轴向速度的偏心程度越大，会对换热产生影响[11]，常规通道内的研究结果并不一定适用于小通道。因此，对小通道内超临界CO2的异常传热机理展开研究具有重要意义。近年来，数值模拟方法越来越广泛地被应用
于超临界流体流动换热研究。直接数值模拟（DNS）可准确对传热进行预测，但计算量巨大，计算时间长，目前仅用于对雷诺数较低的情况进行预测[12]，不适用于工程应用。雷诺平均法（RANS）在对超临界流体流动换热的数值模拟中得到了广泛应用[13]，其中 SST k-ω模型对超临界流体管内流动的预测效果较好[14-15]。同时，该模型结合了k-ε和 k-ω 模型的优点，即 k-ε模型适于预测远离壁面的区域，k-ω模型适于预测近壁面区域[16-18]。鉴于此，本文拟对d=2 mm的竖直上升管内超临界 CO2流动传热特性进行数值模拟，选取多组典型数据，对模型进行验证；然后采用SST k–ω模型研究质量流速、热流密度、压力等因素对传热的影响规律；最后选取典型工况，分析不同轴向截面处各物性的径向分布情况，揭示传热强化与传热恶化的发生过程以及其内在机理。
1数学物理模型
1.1物理模型
为保证加热段流动充分发展，需在加热段前设置绝热段。如图1 所示，超临界CO2 先流经内径 d=2 mm、长度 Liso=280 mm的绝热流动段，在绝热段出口达到流动充分发展的湍流后进入与绝热段直径相等、长度 Lh=440 mm的加热段。绝热段入口设置为质量流入口边界，加热段出口设置为压力出口边界，加热段壁面设置为无厚度、无滑移边界，加热方式采用均匀热流法。在本文工况下，入口的流体焓值取为iin=235 kJ/kg，特征雷诺数 11 878≤Re≤63 694；湍流度 4.01%≤I ≤4.95%，均为湍流流动。
采用三维结构化网格，利用ICEM 软件对网格进行划分，圆形截面采用O型剖分，近壁面网格加密处理，径向网格尺寸以1.15 的比例增长；沿管长方向采用均匀网格，网格划分结果如图2所示。为精确捕捉层流底层中超临界流体的流动，近壁面第1个网格的无量纲壁面距离y+必须小于1。本文工况中，网格的最大无量纲壁面距离 y+=0.06，满足计算要求。利用 Fluent 软件对控制方程进行离散求解。
图2网格划分Fig. 2 Mesh generation
采用控制体积法求解连续方程、动量方程与能量方程，当上述三大方程与湍流参数方程残差小于10−4时，即认为达到收敛要求。设置SIMPLEC 算法进行耦合求解；离散格式均采用二阶迎风格式。考虑到超临界CO2 物性对温度变化较为敏感，且对传热的影响较为显著，为准确计算各离散点处的物性，调用了制冷剂物性查询软件NIST REFPROP对物性进行即时求解。
1.2网格无关性验证与模型验证
1.2.1 网格无关性验证网格密度会直接影响计算准确度，因此在数值模拟前期需要进行网格无关性验证。本文在运行压力 P=8 MPa、质量流速 G=100 kg/(m2 ·s)、热流密度 q=30 kW/m2的工况下，采用SST k-ω 模型进行网格无关性验证，逐步增加网格数量，直到计算得到的壁温−焓值变化曲线随网格数目变化不明显，即可认为此时计算结果与网格数目无关。图 3 给出了网格数量为105×104 ，280×104， 540×104时的计算结果。由图3可以看出，网格数量（105×104 ）较少时，在截面流体焓值 ib<290 kJ/kg时，壁温 Tw 单调上升，但当 ib>290 kJ//kg 后，模拟得到的壁温迅速上升，并出现壁温峰值；当网格数量增加到 280×104 后，在整个焓值区间内，不再出现壁温峰值；进一步增加网格数量到540×104后，与网格数为 280×104 时相比，计算结果无明显变化。为提高计算效率，本文选用280×104 网格数量进行计算。1.2.2 模型验证湍流模型的选择对超临界流体计算的影响很大。为选取合适的湍流模型，本文选择了超临界流体数值模拟常用的SST k-ω 模型、RNG k-ε模型和 Launder-Sharma 低 Re-k-ε 模型，分别对Nathan[19] 与 Song[20] 的实验结果进行了计算，其中包括1 组传热强化工况与2 组传热恶化工况。3个模型的计算结果如图4所示。
由图 4(a) 可以看出，在强化工况下，SST k-ω模型结果与实验结果吻合较好，最大误差为2.8%；而 RNG k-ε 模型和 Launder-Sharma 模型的计算值则均偏离实验值，预测精度较低，最大误差分别为 12% 与 19.6%。由图 4(b) 和图 4(c) 可以看出，在恶化工况下，实验数据有明显的壁温峰值出
现，此时，SST k-ω模型计算得到的壁温曲线中峰值位置与实验结果基本一致，峰值高低有所出入，最大误差为18.9%；峰值外区域壁温与实验结果吻合较好，最大误差为8.7%；而 RNG k-ε 模型和 Launder-Sharma 模型计算得到的壁温曲线均无峰值出现，与实验得到的传热规律不符，误差高达66.7%。综上，为保证计算结果的准确性，本文选择 SST k-ω模型进行数值计算。
2 运行参数对超临界CO2 在竖直上升管内传热特性的影响
2.1热流密度的影响规律
图 5~图 7 分别给出了质量流速G=500 ，750， 1 000 kg/(m2·s)，P =8 MPa时壁温Tw与换热系数h随热流密度的变化规律。图中，ipc 为流体拟临界温度处对应的焓值，Tb为流体温度。由图5可以看出，在质量流速G =500 kg/(m2·s)工况下，热流密度较低（q=50 kW/m2，q/G=0.1）时，壁温最低，变化较平缓，随着流体焓值的提高单调增加，无峰值出现，相应地，其换热系数最高，且在拟临界温度点附近达到峰值，为传热强化工况；随着热流密度的升高（q=75 kW/m2，q/G=0.15），
图 5 G=500 kg/(m2·s) 时超临界 CO2 的换热特性随热流密度的变化
Fig. 5 The variation of heat transfer characteristics of S-CO2 with heat flux under the condition of G=500 kg/(m2·s)图7 G=1 000 kg/(m2·s) 时超临界CO2 的换热特性随热流密度的变化Fig. 7 The variation of heat transfer characteristics of S-CO2 with heat flux under the condition of G=1 000 kg/(m2·s)
拟临界点前开始出现较为平缓的壁温峰值，换热系数随之降低，达到谷值，然后，随着流体焓值的增加，换热迅速恢复，在 ib≈376 kJ/kg 处达到峰值；随着热流密度的进一步升高（q=100 kW/m2，q/G= 0.2），壁温在拟临界点前出现较尖锐的峰值，换热系数也随之迅速下降，随后，随着流体焓值的增加，换热系数逐渐恢复，在 ib≈433 kJ/kg 处达到最大；热流密度较高（q=150， 200 kg/m2 时，q/G=0.3, 0.4），壁温整体水平与壁温峰值显著升高，换热系数急剧降低，且由图5(b) 可以看出，在较高热流密度下，换热系数随流体焓值的变化曲线几乎重合，都处于较低水平，且在恢复区也未见换热系数峰值。由图6可以看出，当质量流速G=750 kg/(m2·s)， q=50 kW/m2（q/G=0.07）时，换热性能最好，换热系数出现了较为明显的尖峰；当q=100 kW/m2（q/G= 0.13）时，换热在低焓值区发生较为微弱的恶化，随后换热及时恢复，并出现换热系数峰值；q=200， 300 kW/m2 （q/G=0.27，0.4）时，壁温出现显著的峰值，换热系数出现谷值，且恢复区换热系数仍处于较低水平。质量流速G=1 000 kg/(m2·s) 的壁温与换热系数随焓值的变化趋势如图7 所示，其换热规律与 G=500， 750 kg/(m2 ·s)时较为一致，此处不再赘述。
2.2运行压力的影响规律
为研究超临界CO2换热特性随运行压力的变化规律，选取1个传热强化工况（G=1 000 kg/(m2·s)， q=100 kW/m2）与 2个传热恶化工况（G=1 000 kg/(m2·s)， q=200 kW/m2；G =500 kg/(m2·s)，q =100 kW/m2 ），在运行压力分别为 P=8，9，10 MPa时对壁温与换热系数随流体焓值的变化规律进行了对比，结果如图 8~图 10 所示。由图8可以看出，在传热强化工况下，不同压力下壁温的变化趋势基本一致，均随着流体焓值的升高呈平稳上升的趋势，且随着压力的升高，壁温略有升高；不同压力下的换热系数在拟临界点附近出现了峰值，随着压力的升高，换热系数峰值有所降低，位置基本不变，传热强化强度减弱。图 9给出了传热恶化工况G=1 000 kg/(m2·s)， q=200 kW/m2下换热特性随压力的变化规律。从中可以看出，在该工况下，不同压力下的壁温曲线趋势呈现出不同的规律：压力P=8 MPa 时，壁温在低焓值区出现了较为平缓的凸起；P=9，10 MPa时，壁温随流体焓值的变化单调上升，无峰值出现。换热系数的趋势较为一致，均在较低焓值处
图9 G=1 000 kg/(m2·s)，q=200 kW/m2 时超临界 CO2 的换热特性随压力的变化Fig. 9 The variation of heat transfer characteristics of S-CO2 with pressure under the condition of G =1 000 kg/(m2·s) and q=200 kW/m2
图 10 G=500 kg/(m2·s)，q=100 kW/m2 时超临界CO2 的换热特性随压力的变化Fig. 10 The variation of heat transfer characteristics of S-CO2 with pressure under the condition of G =500 kg/(m2·s) and q=100 kW/m2
出现换热系数低谷，不同压力下低谷处对应的流体焓值基本一致，然后换热逐渐恢复，出现换热
系数峰值。随着压力的升高，换热系数显著升
高，传热恶化程度减弱。图 10 给出了传热恶化工况G =500 kg/(m2·s)， q=100 kW/m2下换热特性随压力的变化规律。由图 10可以看出，不同压力下的壁温曲线趋势呈现出不同的规律：P=8 MPa时，壁温在低焓值区出现了非常尖锐的峰值；P=9， 10 MPa时，壁温随流体焓值的变化单调上升，无峰值出现。当P=8 MPa时，换热系数在加热段入口处急剧降低，并达到最低值，随后缓慢上升；当P=9， 10 MPa 时，换热系数显著升高，传热恶化程度减弱。
综上可以得出如下结论：在传热强化工况下，
随着压力的升高，传热强化效果有所削弱；在传
热恶化工况下，随着压力的升高，传热恶化情况
也有所改善。这是由于随着压力的升高，压力离临界压力越来越远，定压比热容、密度、动力黏度等物性对温度的敏感程度降低，变化由剧烈变得
平缓，因此由物性变化引起的传热强化与恶化现象均有所减弱。
3 竖直上升管内超临界CO2 的异常传热机理分析
3.1传热强化机理分析
选取传热强化工况P =8 MPa，q =50 kW/m2， G=500 kg(m2·s) ，进行机理分析。图11 给出了上述工况下壁温与对流换热系数随流体焓值的变化曲线。图 11传热强化工况下壁温与对流换热系数随流体焓值的变化曲线Fig. 11 The variation of wall temperature and convective heat transfer coefficient with fluid enthalpy under the condition of heat transfer enhancement
由图 11 可以看出，在该工况下，壁温随流体焓值的升高单调增加，无明显的峰值或低谷出现。随着流体焓值的升高，管壁与流体的温度差逐渐减小，在拟临界焓值附近达到最小，然后又逐渐增大。对流换热系数的变化与之呼应，也呈先增大后减小的趋势，在拟临界焓值附近达到最佳换热效果，即发生传热强化。为研究传热强化机理，分别截取绝热段ib= 232 kJ/kg ，加热段 ib=275 ，300 ，325 ，350 ， 375 kJ/kg截面，得到各截面处参数的径向分布情况如图12所示。其中，横坐标为无量纲壁面距离y+，其定义式如下：式中：τw为壁面剪切应力，Pa，可通过数值模拟软件输出; ρ 为密度； µ为动力黏度。一般认为， y+<5 时工质处于黏性底层区，5≤y+<30 时处于过渡区，30≤y+<60 时处于对数律层，y+≥60 时处于湍流核心区[21]。图 12 给出了管径 d =2 mm时管内流体在传热强化工况下的物性径向分布。图中：u为轴向速度；( ρb− ρ)g· ρ−1 为浮升力，其中ρb 为截面平均密度， g 为重力加速度； cp 为定压比热容；k为湍