10 mm垂直单管通道内超临界水传热弱化现象的实验与数值分析

张 戈，张 昊，顾汉洋，杨燕华

(上海交通大学 核能科学与工程学院，上海 200240)

超临界压力下流体物性奇特的变化规律导致其对流传热现象十分独特。Jackson等[1]和Cheng等[2]分别对超临界流体对流传热进行研究，根据传热机理不同，可分为3种传热现象：大流量小热流下的强迫传热、小流量大热流下的混合对流传热和大流量大热流下的强加速流动传热。

当流量较小、壁面热流较大时，壁面附近流体密度随温度的升高而急剧减小，导致浮升力的作用加强，并对传热产生很大影响，形成传热弱化现象。根据浮升力作用方向不同，混合对流可分为浮升力助推流动(浮升力的方向和流动方向相同)和浮升力阻滞流动(浮升力的方向和流动方向相反)。浮升力作用下的传热弱化现象是超临界流体传热的一重要特征，也是人们非常关心的一个现象。在相关文献中，Sharabi等[3-4]、Wilcox[5]、Speziale等[6]用传统的湍流模型模拟了传热弱化现象，其结果能定性预测传热弱化现象的发生与恢复，但在壁面温度的预测上均偏高。Petukhov等[7]总结了浮升力对超临界流动传热的影响，主要归结为两类：第1类是浮升力直接影响平均速度分布，进而间接影响湍流流动，称为外部效应；第2类是浮升力直接影响湍流流动，称为结构效应。

本工作通过实验获取10 mm单管中超临界水流动的传热弱化数据，采用文献[8]改进的k-ε-kt-εt四方程模型对该实验参数下的管内流动传热特征进行数值分析。

1 实验研究

1.1 实验

上海交通大学建成了适用于超临界流体热工水力实验的SWAMUP回路[9]，SWAMUP回路示意图示于图1。该回路运行压力最高30 MPa，出口温度最高可达550 ℃，主泵流量5 t/h。预热段加热功率最大300 kW，试验段加热功率最大900 kW。稳压器通过顶部的氮气空间保持回路整体的压力。试验段为内径10 mm的单管，其加热长度为2.5 m。由布置在壁面上的49根间隔为5 cm、直径为1 mm的热电偶测量壁面温度。实验测量的各参数的不确定度列于表1。

图1 SWAMUP回路示意图

表1 测量参数的不确定度

1.2 实验数据处理

试验段的热平衡实验表明试验段的保温效果很好。因此，加热的电功率可认为全部加入到试验段中。壁面热流密度为均匀分布，因此内侧壁面热流密度为：

q=UI/πDiL

(1)

外壁面由49根热电偶测量温度。将壁面划分为100层，在假设体积发热量相同的情况下，内壁面温度由99次的式(2)迭代计算得到。

/4kw)((Di/2)2-(Di+1/2)2)-

(qVl/2kw)(Di/2)2ln(Di/Di+1)

(2)

其中：tw为壁面温度；i为壁面划分层数的编号；kw为管道热导率；qVl为体积热流密度。qVl的表达式为：

(3)

管道热导率由平均温度算得，平均温度为：

/2

(4)

主流温度由沿流动方向的均匀焓升计算，传热系数由内壁温度和主流温度计算。

实验工况列于表2。其中，1#工况为低流量高热流密度，2#工况为高流量高热流密度。

表2 实验工况

2 数值模拟分析

假设流动为二维轴对称流动，其在柱坐标下的质量、动量方程为：

同样是金枝玉叶的段誉，第一次来燕子坞吃的那些：“茭白虾仁”“龙井茶叶鸡丁”，看看就教人馋涎欲滴。段誉的当时心理评判是这样的：“鱼虾肉食之中混以花瓣鲜果，色彩既美，自别有天然清香。”

ρρ

(5)

(6)

(7)

其中：r为柱坐标的半径；x为柱坐标的长度；ρ为密度；u、v为不同方向的速度；p为压力；g为重力加速度；μe为有效黏性系数，μe=μ+μt，μ为动力黏性系数，μt为湍流黏性系数。

(8)

其中：Cμ为常数；fμ为考虑壁面效应的衰减函数；k为湍动能；ε为湍动能耗散率。

能量方程为：

ρρ

(9)

其中：h为焓；σt为湍流普朗特数，本文取0.9。

k和ε方程为：

μ

ρPk+ρGk-ρε

(10)

(11)

根据文献[9]对本文中的k-ε模型进行改进，改进后的k-ε模型的参数和方程列于表3。

表3 k-ε模型参数和方程

(12)

(13)

(14)

式中：Ret为湍流雷诺数；y+为无量纲参数；Reb为主流雷诺数；εw为ε方程在壁面处的边界条件；y为壁面法线方向；ν为运动黏度。

在k方程中，对于垂直管流动，Gk为：

λ

ρCPPt+ρCPε

(15)

(16)

CP2=0.62+1.2exp(-3.2Pr)

(17)

(18)

(19)

根据文献[8]给出了改进后的kt-εt模型，其模型参数列于表4。

表4 kt-εt模型参数和方程

计算中，径向和轴向网格结构根据计算工况进行相应调整，保证最靠近壁面的网格的y+<0.2。为得到黏性支层区域的详细特征，在y+<5.0的区域内布置40个网格点。对流项采用QUICK离散格式，速度场和压力场采用SIMPLEC算法。超临界物性根据APWS-95程序生成，因管道沿程压降相对系统压力很小，物性只是随流体温度变化。

3 结果分析与讨论

3.1 壁面温度的预测

图2示出壁面温度实验测量与计算结果的比较。采用的对比模型是FLUENT中常用的低雷诺数模型，分别为YS模型[10]、AB模型[11]、CH模型[12]、AKN模型[13]、LS模型[14]和常用的两方程SST模型[15]。这几种模型均要求近壁面处具有较多的网格。图2计算了两种工况：工况1，p=23 MPa，G=596.6 kg/(m2·s)，q=772.6 kW/m2；工况2，p=25 MPa，G=2 021 kg/(m2·s)，q=1 385 kW/m2。由图2a可见：YS、AB和CH模型所计算的壁面温度较其他模型计算的均高出很多，其计算结果不能预测传热弱化现象的恢复；SST模型能定性预测到壁面温度的第1个峰值，但定量上还是差了很多。而本文模型能很好预测传热弱化现象的产生及恢复，也能定量预测壁面温度。由图2b可见，传统的低雷诺数模型也高估了壁面温度，本文模型能捕捉到沿流动方向上的壁温逐渐上升趋势，而YS和AB模型的计算结果是下降的。

图2 壁面温度实验测量与计算结果的比较

3.2 传热机理

以工况1为例，为了分析传热弱化现象，对5个不同位置处的物性和湍流特征进行分析，结果示于图3。图3中，x=0对应正常传热区域，x=10R和x=50R分别为传热弱化起始位置和弱化过程区域，x=80R和x=140R分别对应传热弱化最严重区域和后期壁面换热能力恢复区域，R为流道半径。

图3 工况1的近壁面区物性和湍流特征

从图3a可看出，在近壁面区域出现了温度的最大梯度。因超临界流体的物性变化剧烈区域为拟临界点附近区域，将拟临界温度±5 ℃称为大物性变化(LPV)区。可见，在上游x=0位置的LPV区范围十分狭小，而下游x=50R和x=80R位置的LPV区的范围明显扩大。对应不同位置的大比热容变化区的分布特征，从图3b可十分明显地得看出这一点，且LPV区的中心位置开始远离壁面，而x=140R位置的LPV区的中心位置比x=80R更靠近壁面。对于常规流体，因近壁面区域的温度变化最为剧烈，其热物性也在该区域变化最剧烈，但超临界流体的物性变化最剧烈的位置由LPV区的位置决定，这是超临界流体的特殊性。由图3c可见，在x=0和x=10R位置，壁面区域存在很大的密度变化，该区域存在显著的浮升力，并显著改变了流体速度分布，但其大变化区域处仍限制在黏性支层区域(y+<5)，且不可能对湍流特性产生显著影响。在下游x=50R位置，密度剧烈变化区域拓宽到整个近壁面区域，壁面温度的上升促使更下游x=80R和x=140R位置近壁面区域的流体密度进一步降低。近壁面区域的低密度流体与中心区域的高密度流体之间的密度差产生近壁面区域的浮升力作用，显著改变了该区域的流动和湍流特性。图3d示出速度的分布特征。x=0和x=10R的密度仅局限在黏性支层，其速度表现为典型的强制对流湍流速度分布特征。在x=50R位置，浮升力的作用显著增强，近壁面区域的流体速度明显升高，使得整个速度分布变得平坦。而随着浮升力作用的进一步增强，近壁面区域的速度进一步增加，在x=80R区域出现典型的混合对流的“M”型速度分布特征，而x=140R位置因湍流和壁面换热的恢复，“M”型速度分布的近壁面最大值有所减低。

工况1下不同x轴位置的剪切应力分布示于图4。浮升力对湍流的影响分为两部分。

1) 浮升力直接影响剪切应力分布

由图4a可见，当浮升力作用较小时，总剪切应力基本呈线性分布。当有浮升力作用时，剪切应力发生改变。因壁面热流的输入，近壁面处和主流区间存在温差，产生如图3c所示的密度差，导致近壁面处的速度变大，相对主流处的速度形成趋于平坦的层流化现象(图3d)。同时，密度差还导致近壁面处的剪切应力减小(图4b、c)。剪切应力和平均速度梯度对湍流传热的影响很大。平均速度梯度的存在导致较大涡的存在，大涡通过剪切应力增强湍流的交混作用，提高热量传递的效率。在浮升力作用下，剪切应力和速度梯度变小，使得湍流能力变弱，并导致传热能力变弱，进而使近壁面处和主流处温差进一步加大，使得流体密度差进一步加大，导致剪切应力和速度梯度进一步减小，形成一正反馈模式，导致传热能力急剧弱化。随着传热能力的弱化，近壁面处和主流处的密度差进一步增大，导致负的剪切应力的出现(图4d、e)，剪切应力的绝对值增加。另一方面，近壁面处的速度变大，出现“M”型速度分布，近壁面处与主流处的速度梯度再次增加。这使得湍流能力加剧，传热能力开始恢复。

图4 工况1不同x轴位置剪切应力分布

2) 浮升力直接影响湍流结构

湍流产生项和浮升力项是k-ε方程中的重要项。图5示出工况1近壁面处的湍流特性。由图5可见，湍流产生项在x=80R处达到最小值，而此处恰好是传热恶化最剧烈的点。当流动状态由湍流向层流转变时，湍流产生项急剧变小，当流动状态由层流向湍流转变时，湍流产生项又急剧增加。浮升力项在x=0位置很小，随着流动状态的变化有一定增加。但因为湍流产生项的减小，导致浮升力项的影响不能被忽略，成为湍动能的一主要源项。

图5 工况1近壁面处的湍流特性

4 结论

1) 改进的低雷诺数模型能准确预测超临界流体流动传热的传热恶化现象发生的起始、峰值和结束位置，而传统模型均高估了壁面温度。

2) LPV区域的位置影响浮升力效应。当LPV区在过渡区时，浮升力明显影响了平均速度和湍流剪切应力，导致了传热弱化现象的产生与恢复。

3) 除浮升力的外部效应影响，其结构效应也会影响湍流产生项和浮升力项，而这两项均是传热弱化情况下湍流方程中重要的项。

参考文献：

[1] JACKSON J D, HALL W B. Forced convection heat transfer to fluids at supercritical pressure: Turbulence forced convection in channels and bundles[M]. New York: Hemisphere Publishing Corporation, 1978: 563-611.

[2] CHENG X, SCHULENBERG T. Heat transfer at supercritical pressures, literature review and application to an HPLWR[R]. Kalsruhe, Germany: FZKA, 2001.

[3] SHARABI M, AMBROSINI W, FORGIONE N, et al. Prediction of experimental data on heat transfer to supercritical water with two-equation turbulence models[C]∥3rd International Symposium on SCWR-Design and Technology. Shanghai, China: [s. n.], 2007.

[4] SHARABI M, AMBROSINI W. Discussion of heat transfer phenomena in fluids at super critical pressure with the aid of CFD models[J]. Annals of Nuclear Energy, 2009, 36(1): 60-71.

[5] WILCOX D C. Turbulence modeling for CFD[M]. California, USA: DCW Industries Inc., 2000.

[6] SPEZIALE C G, ABID R, ANDERSON E C. Critical evaluation of two-equation models for near wall turbulence[J]. AIAA Journal, 1992, 30(2): 324-331.

[7] PETUKHOV B S, POLYAKOVA H. Heat transfer in turbulent mixed convection[M]. New York: Hemisphere Publishing Corporation,1988.

[8] ZHANG H, XIE Z R. Numerical study on supercritical fluids flow and heat transfer buoyancy[C]∥The 8th International Topical Meeting on Nuclear Thermal-hydraulics, Operation and Safety (NUTHOS-8). Shanghai, China: [s. n.], 2010.

[9] GU H Y, ZHANG G, WEN Q L, et al. Supercritical water heat transfer test in a vertical tube[C]∥The 8th International Topical Meeting on Nuclear Thermal-hydraulics, Operation and Safety (NUTHOS-8). Shanghai, China: [s. n.], 2010.

[10] YANG Z, SHIH T H. New time scale basedk-εmodel for near-wall turbulence[J]. AIAA Journal, 1993, 31(7): 1 191-1 198.

[11] ABID R. Evaluation of two-equation turbulence models for predicting transitional flows[J]. Int Eng Sci, 1993, 31(6): 831-840.

[12] CHIEN K Y. Predictions of channel and boundary layer flows with a low Reynolds number turbulence model[J]. AIAA Journal, 1982, 20(1): 33-38.

[13] ABE K, KONDOH T, NAGANO N. A new turbulence model for predicting fluid flow and heat transfer in separating and reattaching flows, Ⅰ: Flow field calculation[J]. Int J Heat Mass Transfer, 1994, 37(1): 139-151.

[14] LAUNDER B E, SHARMA B I. Application of the energy dissipation model of turbulence to the calculation of flow near a spinning disk[J]. Lett Heat Mass Transfer, 1974, 1(2): 131-138.

[15] MENTER F R. Two equation eddy-viscosity turbulence models for engineering application[J]. AIAA Journal, 1994, 32(8): 1 598-1 605.