庄晓如 徐心海 杨智 赵延兴 余鹏†

1) (南方科技大学力学与航空航天工程系,深圳 518055)

2) (哈尔滨工业大学 (深圳)机电工程与自动化学院,深圳 518055)

3) (广东工业大学材料与能源学院,广州 510006)

4) (中国科学院理化技术研究所,中国科学院低温工程学重点实验室,北京 100190)

研究超临界CO2 在高温吸热管内的传热特性是将其应用于聚光太阳能热发电技术中的基础.本文对此进行了数值模拟研究,分析了流体温度、流动方向、系统压力、质量流率和热流密度对对流传热系数和Nu 数的影响.结果表明: 高温区 (800—1050 K) 的对流传热系数和Nu 数受流动方向和系统压力的影响均很小,但都随着质量流率的增大以及热流密度的减小而明显增大; 而随着流体温度的升高,对流传热系数近似线性增大,Nu 数则近似线性减小.另外,本文研究发现在高温区可忽略浮升力对传热的影响,而由高热流密度引起的流动加速效应会明显恶化传热.最后,选取了八种管内超临界流体传热关联式与模拟结果进行对比,发现使用基于热物性修正的关联式对高温区传热数据预测的结果优于使用基于无量纲数修正的关联式得到的结果,且其中预测效果最优的关联式得到的计算结果与模拟结果之间的平均绝对相对偏差为8.1%.

1 引 言

近年来随着聚光太阳能热发电技术的不断发展,工程上对其传热流体提出了更高的要求.CO2具有来源广泛、化学性质稳定、无毒害等优点,并且在超临界状态 (Tc= 304 K,pc= 7.4 MPa) 下具有优越的传热特性[1].同时,超临界CO2布雷顿循环在773—973 K 温区的循环热效率明显高于蒸汽动力循环在此温区的循环效率[2−4].因此,超临界CO2有潜力作为高温传热流体应用于聚光太阳能热发电站中[5,6].美国能源部在第三代聚光太阳能热发电技术路线的报告中也提出,超临界CO2可作为新型高温传热流体应用于热发电站的吸热器中,并指出开发新型高温传热流体的目标之一是使其在吸热器出口温度超过993 K[7].

深入认识超临界CO2在高温区的管内传热特性是将其应用于聚光太阳能热发电技术中的基础,然而目前关于其管内传热特性的实验和数值模拟研究均大多集中在近临界点温区[8,9].Kim 等[10]对超临界CO2在垂直管内径为4.5 mm、入口温度范围为302—388 K、系统压力范围为7.5—10.3 MPa、质量流率范围为208—874 kg·m–2·s–1和热流密度范围为38—234 kW·m–2的工况下进行了对流传热实验,研究了热流密度和质量流率对管壁面温度分布的影响.在向上流动、热流密度适中以及质量流率较低时,壁面温度有明显的峰值; 在向下流动时,所有工况下壁面温度沿流动方向均单调增加; 浮升力和流动加速效应对传热有明显的影响,并考虑了两者的作用提出了新的传热关联式.Bovard 等[11]模拟研究了超临界CO2在垂直管内径为4.5 和6.3 mm、入口温度为302 K、系统压力范围为6.5—8.3 MPa、质量流率范围为51—236 kg·m–2·s–1和热流密度范围为52—85 kW·m–2的管内传热特性,得到提高质量流率和系统压力可削弱浮升力和流动加速效应对传热的影响,也提出了考虑两种效应的新传热关联式.

然而,将超临界CO2作为传热流体应用在塔式聚光太阳能热发电站中,其吸热器入口温度通常要求高于773 K,目前鲜少有关于高温区超临界CO2管内传热特性研究的公开报道[8].Qiu 等[12]对超临界CO2在槽式太阳能吸热器中的传热特性进行了数值计算,所研究的入口温度范围为315—773.15 K、系统压力为9 和20 MPa,结果表明浮升力作用使圆管横截面出现了明显的二次涡流现象,且近临界点温区的二次涡流速度比高温区大一个数量级.刘占斌等[13]模拟研究了六种泡沫材料填充方式对吸热管内超临界CO2流动与换热性能以及管壁温度分布的影响.所研究的入口温度为813 K、系统压力为20 MPa 和热流密度为500 kW·m–2,结果表明环形填充方式 (沿管内壁填充) 的流动换热综合性能最优、净吸热量最大、管壁最高温度最低且温度分布最均匀.

图1 CO2 在p = 7.5 MPa 下的热物性变化 (a) T = 293—1050 K; (b) T = 293—350 K; (c) T = 800—1050 KFig.1.Properties of CO2 at p = 7.5 MPa: (a) T = 293–1050 K; (b) T = 293–350 K; (c) T = 800–1050 K.

为了进一步研究超临界CO2在高温吸热管内的传热特性,本文首先分析CO2在高温区与近临界点温区的热物性区别,再基于ANSYS FLUENT软件建立吸热管内超临界CO2对流换热的三维数值模型,分析流体温度、流动方向、系统压力、质量流率和热流密度等工况条件对对流传热系数和Nu 数的影响规律,揭示不同工况条件下的传热机理.然后,判别浮升力和流动加速效应在高温区对传热的影响情况,并与八种管内超临界流体传热关联式进行对比,找出预测效果较优的关联式.研究结果可为超临界CO2太阳能吸热器的优化设计和安全稳定运行提供理论参考.

2 超临界CO2 的热物理性质

为了了解CO2在近临界点温区 (293—350 K)和高温区 (800—1050 K) 的热物性区别,图1 给出了CO2在压力7.5 MPa、温度范围293—1050 K下的密度、比热、热导率和黏度的变化情况[14].可见,CO2的热物性在近临界点温区随温度剧烈变化,其中密度和黏度在临界点附近急剧下降,比热和热导率则在临界点附近出现极大值.而在高温区,CO2的热物性与近临界点温区相比变化不那么明显 (特别是比热),且热导率和黏度在两个温区呈现相反的变化趋势.这些因素均可能使超临界CO2在两个温区呈现不同的传热特性.另外,图2所示的是在高温区不同临界压力 (7.5—9 MPa)下CO2热物性的变化情况[14].可见,随着压力增大,各热物性均有所提高,但除密度外,压力对热物性的影响均很小.

图2 CO2 在高温区T = 800—1050 K,p = 7.5—9 MPa 下的热物性变化 (a) 密度; (b) 比热; (c) 热导率; (d) 黏度Fig.2.Properties of CO2 at high temperature region of T = 800–1050 K with p = 7.5–9 MPa: (a) Density; (b) specific heat;(c) thermal conductivity; (d) viscosity.

3 传热模型

3.1 数值模型

本文所研究的吸热管,管径为6 mm,管长为500 mm 的圆管,计算过程中忽略圆管壁厚,几何模型如图3 所示.

假设管内流动为稳态,忽略吸热管与外界环境的换热.对于高温区超临界CO2的流动传热数值计算,Qiu 等[12]和刘占斌等[13]的计算结果显示,采用标准k-ε湍流模型[15]可获得准确的计算结果.因此,本文采用标准k-ε湍流模型进行数值计算,模型的控制方程可表示如下.

质量方程:

动量方程:

能量方程:

湍动能k方程:

耗散速率ε方程:

这里

图3 吸热管几何模型Fig.3.Geometric model of the solar receiver tube.

式中,ρ为密度 (kg·m–3);g为重力系数 (m·s–2);u为速度 (m·s–1);δij为克罗内克函数,当i=j时,δij= 1,当i≠j时,δij= 0;p为压力 (Pa);µ为动力黏度 (Pa·s);µt为湍流黏度,µt=ρCµk2/ε;cp为比热 (J·kg–1·K–1);k为湍动能;ε为耗散速率;T为温度 (K);λ为热导率 (W·m–1·K–1);β为热膨胀系数 (K–1);σk和σε分别为湍动能和耗散速率对应的湍流普朗特数;c1,c2,Cµ,Prt均为模型经验常数;湍流模型参数具体设置[12]为σk= 1.0,σε= 1.3,c1= 1.44,c2= 1.92,Cµ= 0.09,Prt= 0.85.

3.2 边界条件和求解设置

模型的边界条件设置如下.

1) 入口: 采用质量流量入口边界条件,速度垂直于入口截面,流量研究范围为G= 200—500 kg·m–2·s–1.

2) 出口: 采用压力出口边界条件,流动出口的绝对压力研究范围为p= 7.5—9 MPa.

3) 壁面: 设置为静态无滑移边界条件,给定恒定热流密度,研究范围为q= 100—800 kW·m–2.

使用ANSYS FLUENT 软件进行吸热管内超临界CO2传热特性的数值计算,采用SIMPEC 算法求解速度与压力的耦合方程,动量和能量方程采用QUICK 格式进行离散,湍动能和耗散速率方程采用二阶迎风格式.模型中的CO2热物性数据来自美国国家标准与技术研究所 (NIST) 的数据库,通过在ANSYS FLUENT 中激活nist-real-gas-model功能进行直接调用.当各控制方程的残差小于10–6,且出口流体温度和壁面温度没有明显变化,则认为计算收敛.

3.3 网格无关性验证和模型验证

图4 吸热管三维模型网格划分 (网格数量: 841000)Fig.4.Mesh generation of the three-dimensional solar receiver tube (grid quantity: 841000).

采用ANSYS ICEM 软件对图3 吸热管三维模型进行网格划分,网格生成效果如图4 所示.由于近壁面位置的温度梯度和速度梯度较大,所以对近壁面处的网格进行了加密,保证近壁面处的无量纲壁面距离 (y+) 小于1.为了对计算网格进行无关性验证,分别对比了网格数量为2652250,1332250,841000,462250 和196000 时 (网格质量均大于0.7),模拟所得的出口截面壁面温度 (Tw,o) 和对流传热系数 (ho),如表1 所列.可见,算例1 和算例3 的Tw,o和ho相对偏差分别仅为0.07%和0.14%.因此,在同时考虑计算精度和运算时间的情况下,最终选择了网格数量为841000 的网格进行后续的数值计算.其中,对流传热系数h和对应的努塞尔数Nu可由下式计算得到:

式中,q为壁面热流密度 (kW·m–2);D为吸热管直径 (m);Tw为壁面温度;Tb为流体温度,可由下式计算得到

其中A为吸热管横截面流体面积 (m2).

表1 网格无关性验证结果Table 1.Verification for grid independence.

另外,为了验证模型的准确性,将三种工况条件下的计算结果与文献[16]中相同工况条件下的实验数据进行了对比,如图5 所示.可见,三种工况条件的计算结果与实验数据均符合较好,且壁面温度和对流传热系数的平均相对偏差分别为0.85%和3.35%,从而验证本文的数值计算模型是可靠的.

图5 模型计算结果与文献[16]实验数据的对比 (a) 壁面温度; (b) 对流传热系数; 其中,工况1,G = 868 kg·m–2·s–1,q = 231.0 kW·m–2,p = 9.22 MPa; 工况2,G = 873 kg·m–2·s–1,q = 216.2 kW·m–2,p = 9.09 MPa; 工况3,G = 874 kg·m–2·s–1,q = 191.8 kW·m–2,p = 8.71 MPaFig.5.Comparisons between numerical results and experimental data of Ref.[16]: (a) Wall temperature; (b) convective heat transfer coefficient.Case 1: G = 868 kg·m–2·s–1,q =231.0 kW·m–2,p = 9.22 MPa.Case 2: G = 873 kg·m–2·s–1,q =216.2 kW·m–2,p = 9.09 MPa.Case 3: G = 874 kg·m–2·s–1,q =191.8 kW·m–2,p = 8.71 MPa.

4 结果与分析

4.1 流体温度和流动方向的影响

实际工程应用中,超临界CO2在吸热管内的流动方向存在多种情况,本文主要研究其在管内进行垂直向上、垂直向下和水平方向流动时的传热特性变化情况.图6 为在G= 300 kg·m–2·s–1,q=500 kW·m–2和p= 7.5 MPa 的工况下,三种流动方向在高温区对管内超临界CO2传热特性的影响.可见,三种流动方向的对流传热系数均随流体温度的升高而近似线性的增大,其中垂直向下流动的对流传热系数最高,而垂直向上流动的对流传热系数最低,但它们的数值偏差不大,其绝对相对偏差小于0.68%.这是因为CO2的密度在高温区会随着流体温度的升高而减小,在定质量流率 (G=ρu)下则使得流速增大,从而增强了湍流强度,导致对流传热系数随流体温度的升高而增大.另外,根据这三种极端流动方向对传热的影响结果,可推断其他倾斜角度的流动方向对高温区超临界CO2吸热管内传热特性的影响应该也不大.从图6 还可得到Nu 数随流体温度的升高是减小的,与对流传热系数的趋势相反,这是因为CO2的热导率在高温区是随流体温度的升高而增大的,从而说明在高温区Nu 数受CO2热导率的影响比对流传热更为明显.

图6 管内流体温度和流动方向在高温区对超临界CO2传热特性的影响Fig.6.Effects of fluid temperature and flow direction on heat transfer of supercritical CO2 at high temperature region.

图7 z = 250 mm 管截面平面上三种流动方向的流体 (a) 温度云图; (b) 径向速度矢量图; (c) 湍动能云图Fig.7.(a) Temperature contours,(b) radial velocity fields,(c) turbulence kinetic energy contours on the plane of z = 250 mm for three flow directions.

图7给出了在z = 250 mm 横截面上的管内流体温度和湍动能云图以及径向速度矢量图.对于垂直向上和垂直向下流动,从图7(a)可见,流体的温度分布均呈径向对称,最低流体温度位于管截面的中心,且向上流动的流体径向温度梯度略高于向下流动的流体径向温度梯度.从图7(b)可见,垂直向上和垂直向下流动的流体径向速度分布均呈现M 形.这是由于径向温度梯度的存在使得管截面流体产生径向密度梯度,管中心流体由于温度较低密度增大,导致其流速降低,而近壁面流体由于温度较高密度减小,反而促使其流速增大,从而形成“M”形径向速度分布.另外,对于吸热管内向上流动,浮升力效应会降低边界层剪切力,减弱湍流涡流的产生和扩散[17],从而削弱对流换热.而对于向下流动,浮升力效应的影响则相反,所以向下流动的对流传热系数高于向上流动的对流传热系数.但从图7(c)可见,它们的湍动能分布差别不大,所以它们的对流传热系数也偏差不大.而对于水平流动,从图7(a)可见,其流体温度分布并非呈径向对称,最低流体温度位置是略微的向管底部偏移.从图7(b)可见,在水平流动的流体径向速度分布中,管截面两边的流体是沿管壁面向上运动,然后再沿管中心下降,从而在管截面两边形成二次涡流.而二次涡流的生成,使得最低流体温度位置向管底部偏移.另外,水平流动的流体径向速度比垂直向上和向下流动的径向速度大一个数量级,且水平流动的管截面两侧流体径向速度明显高于其他管截面区域.从图7(c)可见,虽然水平流动的流体径向速度较高,但其湍动能反而比垂直向下流动略低.这可能是因为水平流动中剪切力大多用于加速管截面两侧流体速度,消耗了其湍动能.综上所述,三种流动方向的湍流度和传热效率均相差不大.

4.2 系统压力的影响

图8 为在G = 300 kg·m–2·s–1,q = 500 kW·m–2和垂直向上流动的工况下,系统压力在高温区对管内超临界CO2传热特性的影响.可见,系统压力在高温区对对流传热系数的影响极小,随着系统压力的增大,对流传热系数最多仅减小0.17%; 而Nu 数随系统压力的增大也仅有略微的减小,最大减小量为0.51%,且流体温度越高系统压力的影响越小.由图2 可见: 在高温区CO2的密度会随着系统压力的增大而增大,在定质量流率下,会使得流体流速减小,从而削弱传热; 但另一方面CO2的热导率也会随着系统压力的增大而增大,从而强化传热.两者综合作用,使得对流传热系数的变化不大.然而,图2 也显示在高温区系统压力对CO2密度的影响明显强于热导率,但它们对对流传热系数的影响却相当,说明在高温区和本文所研究的系统压力范围 (7.5—9 MPa),CO2热导率对其传热特性起着重要的作用,这也与4.1 节中Nu 数随流体温度增大而减小的原因一致.

图8 系统压力在高温区对超临界CO2 传热特性的影响Fig.8.Effect of pressure on heat transfer of supercritical CO2 at high temperature region.

4.3 质量流率的影响

图9 质量流率在高温区对超临界CO2 传热特性的影响Fig.9.Effects of mass flux on heat transfer of supercritical CO2 at high temperature region.

图9为在p = 7.5 MPa,q = 500 kW·m–2和垂直向上流动的工况下,流体质量流率在高温区对管内超临界CO2传热特性的影响.可见,高温区的对流传热系数和Nu 数均随着质量流率的增大而明显增大.这是因为质量流率越大,流体流速也越大,使得流体的湍流扩散速率增大,从而增强了对流换热效应,导致对流传热系数增大.另外,管壁面边界层厚度也会随着流体流速增大而减小,导致热阻变小,同样强化了传热.

4.4 热流密度的影响

图10 为在G = 300 kg·m–2·s–1,p = 7.5 MPa和垂直向上流动的工况下,壁面热流密度在高温区对管内超临界CO2传热特性的影响.可见,在高温区随着热流密度的增大,对流传热系数和Nu 数均减小,且热流密度越大其影响程度越明显.由于随着壁面热流密度增大,近壁面的流体温度会急剧升高,而流体温度升高会造成其体积膨胀,引发流动加速效应[17].流动加速效应所形成的速度梯度方向与流动方向是相反的[18],这会导致近壁面处的湍流度减小.当近壁面处的低密度流体层达到足够厚,会使得湍流边界层转化为层流边界层,从而导致对流传热被削弱.

4.5 浮升力和流动加速效应的影响

根据现有文献[19−22]可知,浮升力和流动加速效应在近临界点温区对管内超临界CO2传热特性具有重要的作用.然而,目前不同文献对两种效应的定义各有不同,本文采用Liu 等[17]提出的判断标准,分析在高温区浮升力和流动加速效应对管内超临界CO2传热特性的影响情况.本文所使用的浮升力因子 (Bu) 和流动加速因子 (Ac) 定义如下:

图10 热流密度在高温区对超临界CO2 传热特性的影响Fig.10.Effect of heat flux on heat transfer of supercritical CO2 at high temperature region.

式中,

Liu 等[17]指出当Bu< 1.3×10–5和Ac<3.3×10–6时,浮升力和流动加速效应对传热的影响小于3%,可以忽略不计.本文计算了不同工况条件下的Bu和Ac,如图11 所示.如图11(aI) 和图11(aII) 所示,在高温区随着流体温度的升高,Bu和Ac均减小,流动方向对它们几乎没有影响,这与流动方向对对流传热系数的影响规律一致(图6).如图11(bI) 和图11(bII) 所示,随着系统压力增大,Bu会增大,而Ac几乎不变.但所研究系统压力范围的Bu值均远低于其阈值1.3×10–5,说明此时浮升力效应对传热的影响不大,因此即使Bu值随着系统压力的增大而增大,对流传热系数受系统压力的影响仍很小 (图8).如图11(cI) 和图11(cII) 所示,Bu和Ac均随着质量流率的增大而减小,这是因为Re数会随着质量流率的增大而增大,由 (11) 式和 (12) 式可见,Re数的增大会减小Bu和Ac.而对于垂直向上流动,浮升力和流动加速效应会恶化传热,Bu和Ac值的减小,使得流体流速增大强化传热的程度高于浮升力和流动加速效应恶化传热的程度,从而促使对流传热系数随着质量流率的增大而增大 (图9).如图11(dI) 和图11(dII) 所示,随着热流密度的增大,Bu和Ac均增大,这是因为壁面和流体的温差会随着热流密度的增大而增大,使得管内流体密度梯度增大,从而导致浮升力和流动加速效应增强,削弱了对流传热,这与热流密度对对流传热系数的影响一致 (图10).

图12 为在本文所研究的工况条件下,Bu和Ac的分布图.可见,所有工况的Bu值均远低于其阈值1.3×10–5,而除了最低热流密度工况 (G=300 kg·m–2·s–1,q= 100 kW·m–2,p= 7.5 MPa和垂直向上流动) 外,其余工况的Ac值均高于其阈值3.3×10–6.说明在本文所研究的工况条件(高温高热流密度) 下可忽略浮升力对传热的影响,而由高热流密度引起的流动加速效应会明显恶化传热.

4.6 与现有文献传热关联式计算结果的对比

现有文献中关于管内超临界CO2的传热关联式大多是基于近临界点温区的传热数据提出的,其对高温区传热数据的预测还有待评估.表2 为目前常用的八种管内超临界流体传热关联式,其可分为基于热物性修正和基于无量纲数修正两种类型.将关联式计算结果与模拟结果进行对比,并采用平均绝对相对偏差 (mean absolute relative deviation,MARD) 和平均相对偏差在±20%以内的数据点占总数据量的比例 (η) 作为关联式预测能力的评估标准,MARD 具体定义如下:

图11 (a) 流动方向、(b) 系统压力、(c) 质量流率、(d) 热流密度在高温区对 (I) Bu 和 (II) Ac 的影响Fig.11.Effects of (a) flow direction,(b) pressure,(c) mass flux,(d) heat flux on (I) Bu and (II) Ac at high temperature region.

图12 Bu 和Ac 的分布图Fig.12.Distributions of Bu and Ac .

式中Nucal,Nunum为传热关联式计算、模拟结果.

图13 为其中预测效果较好的五种关联式计算结果与模拟数据对比的散点图.从表2 和图13 可得,基于热物性修正的传热关联式对本文高温区数据预测的结果优于基于无量纲数修正的关联式得到的结果,且模拟数据均大于计算结果.由前面的分析可得,高温区的浮升力效应对传热的影响可忽略,而虽然流动加速效应对传热有一定的影响,但Nu数受导热传热的影响比对流传热的影响要大,因此高温区的管内超临界CO2传热数据的预测使用基于热物性修正的传热关联式更合适.其中,Krasnoshchekov 和Protopopov[24]的关联式对本文传热数据预测效果最好,其平均绝对相对偏差为8.1%,且100%的数据点在±20%的平均相对偏差范围内.

图13 传热关联式计算结果与模拟数据的对比Fig.13.Comparisons of the calculated heat transfer results by the correlations with the numerical results.

表2 传热关联式计算结果与模拟数据的对比Table 2.Comparisons of the calculated heat transfer results by the correlations with the numerical results.

表2 （续） 传热关联式计算结果与模拟数据的对比Table 2 (continued).Comparisons of the calculated heat transfer results by the correlations with the numerical results.

5 结 论

本文对超临界CO2在高温吸热管内的传热特性进行了数值模拟研究,分析了流体温度、流动方向、系统压力、质量流率和热流密度对对流传热系数和Nu数的影响,揭示了不同工况条件下的传热机理,探究了浮升力和流动加速效应在高温区对传热的影响情况,并选取了八种常用的管内超临界流体传热关联式与模拟结果进行了对比,主要得到以下结论:

1) 在高温区 (800—1050 K),随着流体温度的升高,对流传热系数近似线性的增大,而Nu数近似线性的减小;

2) 垂直向上、垂直向下和水平三种流动方向以及系统压力对高温区的对流传热系数和Nu数影响均很小;

3) 高温区的对流传热系数和Nu数均随着质量流率的增大和热流密度的减小而明显的增大,且热流密度越大其影响程度越明显;

4) 在本文所研究的工况范围内,浮升力效应对传热的影响可忽略,而由高热流密度引起的流动加速效应会明显恶化传热;

5) 高温区的管内超临界CO2传热数据的预测使用基于热物性修正的传热关联式优于基于无量纲数修正的关联式,且Krasnoshchekov 和Protopopov[24]的关联式预测效果较优,对模拟结果预测的平均绝对相对偏差为8.1%.