竖直上升管内超临界CO2异常传热机理研究

胡娟娟，李增光，高一民，杨建，刘浪

中国舰船研究设计中心 上海分部，上海 201108

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]，其中SSTk-ω模型对超临界流体管内流动的预测效果较好[14-15]。同时，该模型结合了k-ε和k-ω模型的优点，即k-ε模型适于预测远离壁面的区域，k-ω模型适于预测近壁面区域[16-18]。

鉴于此，本文拟对d=2 mm的竖直上升管内超临界CO2流动传热特性进行数值模拟，选取多组典型数据，对模型进行验证；然后采用SSTk-ω模型研究质量流速、热流密度、压力等因素对传热的影响规律；最后选取典型工况，分析不同轴向截面处各物性的径向分布情况，揭示传热强化与传热恶化的发生过程以及其内在机理。

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%，均为湍流流动。

图 1 物理模型与边界条件设置Fig.1 Physical model and the setting of boundary condition

采用三维结构化网格，利用ICEM软件对网格进行划分，圆形截面采用O型剖分，近壁面网格加密处理，径向网格尺寸以1.15的比例增长；沿管长方向采用均匀网格，网格划分结果如图2所示。为精确捕捉层流底层中超临界流体的流动，近壁面第1个网格的无量纲壁面距离y+必须小于1。本文工况中，网格的最大无量纲壁面距离y+=0.06，满足计算要求。

图 2 网格划分Fig.2 Mesh generation

利用Fluent软件对控制方程进行离散求解。采用控制体积法求解连续方程、动量方程与能量方程，当上述三大方程与湍流参数方程残差小于10-4时，即认为达到收敛要求。设置SIMPLEC算法进行耦合求解；离散格式均采用二阶迎风格式。考虑到超临界CO2物性对温度变化较为敏感，且对传热的影响较为显著，为准确计算各离散点处的物性，调用了制冷剂物性查询软件NIST REFPROP对物性进行即时求解。

1.2 网格无关性验证与模型验证

1.2.1 网格无关性验证

网格密度会直接影响计算准确度，因此在数值模拟前期需要进行网格无关性验证。本文在运行压力P=8 MPa、质量流速G=100 kg/(m2·s)、热流密度q=30 kW/m2的工况下，采用SSTk-ω模型进行网格无关性验证，逐步增加网格数量，直到计算得到的壁温-焓值变化曲线随网格数目变化不明显，即可认为此时计算结果与网格数目无关。

图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网格数量进行计算。

图 3 网格无关性验证Fig.3 Mesh independence validation

1.2.2 模型验证

湍流模型的选择对超临界流体计算的影响很大。为选取合适的湍流模型，本文选择了超临界流体数值模拟常用的SSTk-ω模型、RNGk-ε模型和Launder-Sharma低Re-k-ε模型，分别对Nathan[19]与Song[20]的实验结果进行了计算，其中包括1组传热强化工况与2组传热恶化工况。3个模型的计算结果如图4所示。

图 4 模型验证Fig.4 Model validation

由图4(a)可以看出，在强化工况下，SSTk-ω模型结果与实验结果吻合较好，最大误差为2.8%；而RNGk-ε模型和Launder-Sharma模型的计算值则均偏离实验值，预测精度较低，最大误差分别为12%与19.6%。由图4(b)和图4(c)可以看出，在恶化工况下，实验数据有明显的壁温峰值出现，此时，SSTk-ω模型计算得到的壁温曲线中峰值位置与实验结果基本一致，峰值高低有所出入，最大误差为18.9%；峰值外区域壁温与实验结果吻合较好，最大误差为8.7%；而RNGk-ε模型和Launder-Sharma模型计算得到的壁温曲线均无峰值出现，与实验得到的传热规律不符，误差高达66.7%。综上，为保证计算结果的准确性，本文选择SSTk-ω模型进行数值计算。

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)时超临界CO2的换热特性随热流密度的变化Fig.5 The variation of heat transfer characteristics of S-CO2 with heat flux under the condition of G=500 kg/(m2·s)

图 6 G=750 kg/(m2·s)时超临界CO2的换热特性随热流密度的变化Fig.6 The variation of heat transfer characteristics of S-CO2 with heat flux under the condition of G=750 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)

由图5可以看出，在质量流速G=500 kg/(m2·s)工况下，热流密度较低（q=50 kW/m2，q/G=0.1）时，壁温最低，变化较平缓，随着流体焓值的提高单调增加，无峰值出现，相应地，其换热系数最高，且在拟临界温度点附近达到峰值，为传热强化工况；随着热流密度的升高（q=75 kW/m2，q/G=0.15），拟临界点前开始出现较为平缓的壁温峰值，换热系数随之降低，达到谷值，然后，随着流体焓值的增加，换热迅速恢复，在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时，壁温随流体焓值的变化单调上升，无峰值出现。换热系数的趋势较为一致，均在较低焓值处出现换热系数低谷，不同压力下低谷处对应的流体焓值基本一致，然后换热逐渐恢复，出现换热系数峰值。随着压力的升高，换热系数显著升高，传热恶化程度减弱。

图 8 G=1 000 kg/(m2·s)，q=100 kW/m2时超临界CO2的换热特性随压力的变化Fig.8 The variation of heat transfer characteristics of S-CO2 with pressure under the condition of G=1 000 kg/(m2·s) and q=100 kW/m2

图 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 传热强化工况下的径向物性分布Fig.12 Radial physical properties distribution under the condition of heat transfer enhancement

图12给出了管径d=2 mm时管内流体在传热强化工况下的物性径向分布。图中：u为轴向速度；( ρb- ρ)g· ρ-1为浮升力，其中ρb为截面平均密度，g为重力加速度；cp为定压比热容；k为湍动能；λ为导热系数。由图12(a)可以看出，随着y+的增大，加热段流体温度逐渐降低，并在轴线处达到最低，随着温度的变化，流体物性也发生变化；绝热段温度不变，相应地，物性也不变。由图12(b)可知，在近壁面处，流体密度较低，为保持质量守恒，流体迅速向上流动，近壁面处流速提高，如图12(c)所示；同时，由于近壁面处流体密度较低也导致此处浮升力较高，如图12(d)所示，较高的浮升力也使得近壁面处流体速度有所提高。由于湍动能与雷诺切应力正相关，雷诺切应力又随速度梯度的增大而增大，因此随速度梯度先增大后减小，湍动能也呈现先增大后减小的趋势，如图12(e)所示。同时由图12(e)可知，随着流体焓值的增加，湍动能不断增强，湍动能峰值越来越高。

当ib=275～350 kJ/kg时，截面流体会达到或跨越拟临界温度，此时，大比热容区形成。如图12(f)所示，随着ib的增加，大比热容区逐渐向远离壁面的方向扩散，大比热容区流体占据的份额越来越大，当ib=350 kJ/kg时，大比热容区开始扩散到过渡区，大比热容区流体份额达到最大，此时换热系数达到峰值。ib＞350 kJ·kg-1后，截面流体温度均已跨越拟临界温度，此时，比热容cp单调变化，峰值消失，大比热容区流体所占份额减少，换热系数有所降低。

为定量计算各截面流体平均比热容水平，引入截面平均比热容cp,ave，其定义如下：

式中：r为截面半径；uc为工质流速；A为截面面积。

由式(2)，求得6个截面处的截面平均比热容如表1所示。可以看出，随着流体焓值的升高，截面平均比热容先迅速升高，在ib=350 kJ/kg时达到最高，然后随着截面流体焓值的增加，平均比热容有所降低。这与上文的定性分析结果一致。

表 1 各截面的平均比热容Table 1 Average specific heat in each section

3.2 传热恶化机理分析

选取数值模拟工况P=8 MPa，G=500 kg/(m2·s)，q=100 kW/m2作为传热恶化典型工况进行分析。图13给出了该工况下壁温与对流换热系数随流体焓值的变化曲线。由图13可以看出，进入加热段后，壁温急剧升高，传热恶化现象逐渐发生，并在ib=300 kJ/kg附近达到峰值，然后逐渐降低，变化趋于平缓，传热逐渐恢复；与此对应，进入加热段后换热系数急剧降低，在ib=300 kJ/kg附近达到谷值，然后缓慢升高，换热有所恢复。可以看出，ib=300 kJ/kg为该工况下传热恶化最为剧烈的地方，因此在接下来的机理分析中将重点分析此处的物性变化特征。

图 13 传热恶化工况下壁温与对流换热系数随流体焓值的变化曲线Fig.13 The variation of wall temperature and convective heat transfer coefficient with fluid enthalpy under the condition of heat transfer deterioration

为研究传热恶化机理，分别截取绝热段ib=232 kJ/kg，加热段ib=250，275，300，350，400 kJ/kg这6个截面，得到各截面处参数的径向分布情况如图14所示。由图14(a)可见，各截面处壁面温度均高于拟临界温度，与传热强化工况相比，在传热恶化工况下，径向温度变化十分剧烈，在ib由235 kJ/kg增大到300 kJ/kg的过程中，壁温越来越高，径向温度梯度也越来越大。受到温度分布的影响，近壁面低密度区形成，并在ib=300 kJ/kg处达到最低值，如图14(b)所示。在低密度区域的影响下，近壁面浮升力增大使得流体加速，如图14(c)与图14(d)所示。随着速度梯度的增大，黏性底层区的湍动能有所增加。随着流动发展到过渡区，流体密度的急剧上升使浮升力迅速下降，速度沿径向的变化逐渐变得平缓，此时，雷诺切应力减小，因此湍流强度逐渐减弱。由图14(c)可以看出，ib=300 kJ/kg时，出现了速度拐点，拐点处速度梯度降为0，由于管内流场为轴对称分布，此时管内呈现“M”型速度曲线；与之对应，湍动能也达到最低，此时传热传质的剧烈程度降低，造成热量积聚，发生传热恶化。在ib＞ 300 kJ/kg后，“M”型速度曲线消失，速度梯度增大使湍动能升高，换热性能得到改善，壁温有所降低。

4 结 论

本文利用经验证的数值模拟方法，对超临界CO2换热特性随不同实验参数的变化规律进行了研究，并选取典型工况深入分析了异常传热机理。在本文研究的参数范围内，得到如下结论：

图 14 传热恶化工况下的径向物性分布Fig.14 Radial physical properties distribution under the condition of heat transfer deterioration

1） 热流密度较低（q/G＜0.13）时，壁温沿管长方向单调增长，换热系数在拟临界点附近达到峰值；热流密度较高（q/G≥0.13）时，会出现壁温峰值，且随着热流密度的增加，传热恢复区换热系数峰值逐渐降低，并逐步向流体焓值较大的区域推移，直至消失。

2） 随着压力远离临界压力，超临界流体物性变化剧烈程度减弱，由此引起的传热强化与传热恶化程度也随之削弱。

3） 传热强化工况下，径向湍动能随着流体焓值的增大而增大，截面流体平均比热容随流体焓值的增大先增大后减小，在ib=350 kJ/kg时达到最大，此时换热达到最强，因此大比热容区流体份额显著增加是传热强化工况下换热性能变好的主要原因。

4） 传热恶化工况下，受流体径向密度的影响，近壁面浮升力有所增大，使得流体加速流动，随着流动发展到过渡区，流体密度急剧上升使浮升力迅速下降，速度也随之下降。在ib=300 kJ/kg时，出现了“M”型速度曲线，在拐点处速度梯度降为0，湍动能也达到最低，发生传热恶化，此时，速度梯度降低是发生传热恶化的主要原因。