流道结构对超临界二氧化碳传热特性的影响研究

白一凡,吕海财,王 汉

(华北电力大学 核科学与工程学院,北京 102206)

超临界二氧化碳(SCO2)是能源动力系统中最有前景的传热和能量转换工质之一,在工业系统中得到了越来越广泛的应用。在核能发电领域,SCO2气冷快堆是国内外的研究热点,这一系统采用SCO2作为一回路冷却剂吸收燃料元件的裂变热,流出反应堆压力容器后直接进入汽轮机做功[1-2]。相较于传统压水堆,SCO2气冷快堆体积小、系统简单、热效率高,同时易于实现模块化,可大幅降低核电厂建造成本。

堆芯燃料组件内SCO2的流动与传热行为是反应堆热工水力设计中十分关注的问题。燃料组件一般由若干燃料棒按一定的设计方式排列,形成大规模棒束结构。但是在实际研究中,特别是在实验研究方面,研究人员经常针对圆管、方管、环形通道等简单通道开展超临界流体的传热特性实验测量,获得壁面温度和传热系数的变化规律,验证或建立传热关系式。Mori等[3]对超临界压力下氟利昂在圆管、3棒束和7棒束内的传热特性进行了实验研究,发现在相同的实验条件下,圆管出现了传热恶化现象,而在棒束通道中则没有发生。Kim等[4]对热力当量直径为4.5 mm的环形通道内SCO2的传热特性进行了实验测量,并将实验结果与内径为4.4 mm的圆管进行了对比,发现当发生传热恶化时,环形通道的传热系数略微高于圆管。Bae[5]对圆管和环形通道内SCO2的传热系数进行了实验测量,圆管内径为4.57 mm,环形通道加热管直径为8 mm,间隙为2 mm。对比结果表明,在相同的实验条件下,环形通道相比于圆管对传热恶化有一定的抑制作用。Li等[6]在上海交通大学超临界水多用途实验回路上对圆管、环形通道和2×2棒束通道进行了超临界水的传热实验研究,发现2×2棒束的传热性能优于圆管和环形通道。Liu等[7]使用低雷诺数湍流模型对圆管和环形通道中超临界水的传热恶化现象进行了数值模拟研究,发现在高质量流速条件下,圆管的传热恶化现象比内壁加热的环形通道严重,但是当环形通道的内外壁面同时加热时,其传热恶化壁面温度升高的幅度和变化规律与圆管相似。研究还认为,相比于水力当量直径,热力当量直径对高质量流速下传热恶化的发生起重要作用。

目前国内外对超临界流体异常传热现象的实验或数值模拟研究以圆管和环形通道为主,但是由于几何结构的差异,简单通道难以模拟燃料组件内冷却剂的横向交混以及燃料元件周向传热不均匀等流动和传热特点,现阶段缺乏由流道几何结构不同所导致的传热差异的研究。本文采用大涡模拟(LES)方法对SCO2在相同热力当量直径、相同工况条件下的圆管、环形通道和子通道中的流动传热特性进行计算研究,分析流道结构对壁面温度、平均速度分布、湍流统计量等结果的影响,并根据LES结果对典型SCO2传热关系式进行评估。典型的SCO2直接循环反应堆工作压力约为20 MPa,堆芯进口冷却剂温度为400 ℃。在设计工况下,冷却剂的热物性变化较为平缓,不容易发生传热恶化。但是考虑到实际运行中系统参数的瞬态变化,如果发生破口事故,系统参数降至CO2临界点附近时,传热特性与设计工况有所不同,可能会发生传热恶化。因此,本文将计算工况设置在临界点附近来分析此时CO2的传热特性和不同流道对传热的影响。

1 数值模拟设置

1.1 计算参数选取依据

首先,不同几何流道的SCO2入口温度T0相同,保证入口的热物性参数相同;其次,入口位置的流动应满足相似性准则的要求,即3个几何流道的雷诺数相等:

(1)

式中:Re0为入口雷诺数;ρ0为入口流体密度,kg/m3;U0为入口平均流速,m/s;μ0为入口动力黏度,Pa·s;D为热力当量直径,m,其定义如下:

(2)

式中:A为通道横截面积,m2;P为通道横截面加热壁面周长,m。

此外,不同几何形状流道的格拉晓夫数与无量纲热流密度对应相等,这两个无量纲数的定义如下:

(3)

式中:Gr0为入口格拉晓夫数;g为重力加速度,m/s2;β0为入口体积膨胀系数,1/K;λ0为入口导热系数,W/(m·K);ν0为入口运动黏度,m2/s;qw为加热壁面热流密度,W/m2。

(4)

式中:q+为无量纲热流密度;cp0为入口定压比热容,J/(kg·K)。

根据以上分析,不同几何流道的U0和qw应选择相同的参数,同时选取适当几何尺寸以保证D相同,这样能够使数值模拟结果的差异只来源于通道几何结构。

1.2 控制方程

本文使用LES方法研究SCO2的湍流传热特性,控制方程采用低马赫数不可压缩假设,忽略热物性随压力的变化,认为SCO2的物性只与温度有关。控制方程包括质量守恒方程、动量守恒方程和能量守恒方程,表示如下:

(5)

(6)

(7)

式中:下标i、j为张量的阶数;ρ为密度;ui、uj为速度矢量;t为时间;xi、xj为位移矢量;μ为动力黏度;p为压力;gi为重力加速度矢量;μSGS为亚格子尺度(SGS)黏度,需要用SGS模型进行计算;h为比焓;α为热扩散系数;αSGS为SGS热扩散系数,可由SGS黏度确定:

(8)

式中,PrSGS为SGS普朗特数,计算中采用OpenFOAM程序的默认常数1.0。

对于LES计算,最常用的SGS模型是Smagorinsky模型[8],SGS黏度的定义如下:

(9)

(10)

本文还选取了壁面自适应局部涡粘(WALE)模型[9]进行SCO2的传热计算,该模型的选择是基于Brockmeyer等[10]以及Busco等[11]的研究。WALE模型的黏度定义如下:

(11)

(12)

(13)

(14)

(15)

1.3 几何结构和边界条件

图1所示为LES选取的3种几何结构,分别为圆管、环形通道和棒束子通道。圆管直径d为5.276 mm;环形通道内管外径为5.276 mm,外管内径为7.461 4 mm;子通道呈四边形布置,加热棒直径为9.5 mm,栅距p为10.5 mm。圆管壁面和子通道棒束壁面均为加热壁面,而环形通道只有内管外表面加热,外管内壁面绝热。根据上述几何结构尺寸的选择,3种通道的D均为5.276 mm。3种通道的流体质量流速相等,D相等,q+相等。圆管与环形通道质量流量相等,但略低于子通道的质量流量。

a——圆管;b——环形通道;c——子通道

以圆管为例,边界条件的设置如图2所示。流动方向设为x方向,横向方向设为y和z方向。3种通道沿流动方向的加热段长度L均为45D。加热段壁面为恒定热流密度边界条件,壁面处速度设置为无滑移边界条件。加热段出口设置为对流边界条件,即∂φ/∂t+Ub·∂φ/∂x=0,其中:φ为标量或矢量;x为位移矢量;Ub为流体平均速度。加热段入口前设置长度L为6D的湍流入口发生器,其壁面绝热,出入口设置为周期性边界条件,以生成充分发展的湍流。在湍流入口发生器横截面采集足够长时间的速度分量,形成入口截面的湍流数据库,然后按设定的时间步长将其添加到加热段的入口,模拟真实且充分发展的入口湍流。在子通道中,4个棒间隙被分为两对,一对在y方向,另一对在z方向,如图1c所示。每对间隙采用周期性边界条件,允许流体通过间隙自由进出子通道区域。

1.4 计算工况

3种通道的LES计算工况均相同,即压力p0=8.0 MPa,U0=0.204 m/s,Reb,0=13 028,T0=301.15 K,qw=36 kW/m2。3种几何通道的流动均为混合对流工况,g的方向与流动方向相反,大小为9.8 m/s2。

1.5 网格划分

圆管、环形通道与子通道横截面的网格划分如图3所示,图中数字表示箭头所示方向的网格划分数量。圆管径向划分87个网格,周向划分136个网格,轴向划分750个网格,网格数为8 007 000。环形通道径向划分70个网格,周向划分180个网格,轴向划分645个网格,网格数为8 127 000。子通道径向划分75个网格,棒束周向划分376个网格,轴向划分720个网格,网格数为1 560万。网格分辨率根据Pandey等[12]的计算方法,代入本文LES工况,计算得到圆管径向、周向和轴向网格分辨率分别为0.24～4.8、18和48;环形通道径向、周向和轴向网格分辨率分别为0.18～3.6、14～20和55;子通道径向、周向和轴向网格分辨率分别为0.225～4.5、12和50。根据Choi等[13]的研究结论,网格分辨率满足LES的要求。

a——圆管;b——环形通道;c——子通道

1.6 计算方法

本文使用开源程序OpenFOAM求解LES的质量、动量和能量守恒方程。SCO2的热物性参数取自NIST数据库[14],制作物性随温度变化的一维查询表并植入到OpenFOAM中。数值模拟设置方面,使用PISO算法求解速度与压力的耦合,压力场采用GAMG求解器,残差为1×10-7;速度场采用smoothSolver求解器,残差为1×10-6。动量和能量方程中的对流项采用二阶中心差分格式离散,时间项采用二阶向后隐式格式离散。无量纲时间步长Δt+=4.7×10-5,其定义为Δt+=Δt/t+,其中无量纲时间t+=D/uτ,并控制时间步长使计算过程中最大库朗数小于0.4。当入口的湍流发展到出口后开始时间平均,采样时间约为7倍的流体流通时间,对应约300个湍流特性时间(D/U0)。

1.7 验证

本文从两方面对LES计算结果的准确性进行验证。图4所示为子通道入口湍流发生器充分发展的近壁速度与Moser等[15]的DNS数据对比。近壁速度的采样线从棒壁面到子通道中心,横坐标y+=uτy/ν为网格中心到壁面的无量纲距离,纵坐标U+=U/uτ为无量纲速度。由图4可见,LES计算的速度与DNS结果吻合较好,说明本文使用的OpenFOAM能够准确预测绝热条件下的流体流动特征,同时网格分辨率也满足LES的计算要求。

图4 近壁速度与DNS结果对比

另一方面,选取李志辉等[16]的圆管SCO2实验数据验证LES计算结果。圆管垂直布置,直径为2 mm,长度为290 mm。CO2入口温度为298.15 K,压力为8.8 MPa,热流密度为20 263 W/m2。图5所示为Smagorinsky模型和WALE模型计算的壁面温度和流体温度与实验数据的对比。计算发现,所选取的两个模型均能准确预测流体温度沿轴向的变化,但是对壁面温度的预测存在一定偏差。Smagorinsky模型计算得到的壁面温度在进出口段与实验数据吻合较好,但是无法预测传热恶化时的壁面温度峰值。WALE模型计算得到的壁面温度虽然整体偏高于实验数据,但壁面温度分布趋势与实验完全一致,壁面温度峰值与谷值位置也相同。因此,本文认为OpenFOAM中的WALE模型能够较为合理地预测超临界流体的传热规律,选取WALE模型进行后续计算研究,其所预测的壁面温度偏差可能来源于恒定湍流普朗特数的假设。

图5 LES结果与实验数据对比

2 计算结果分析

2.1 壁面温度分布

图6给出了圆管、环形通道和子通道的壁面温度沿轴向的分布。可以看出,3种通道壁面温度分布趋势大致相似,在入口附近存在温度峰值,随着流动向下游发展,峰值逐渐消失,随后壁面温度沿流动方向逐渐升高。环形通道壁面温度峰值略高于圆管,但远高于子通道。SCO2在圆管和环形通道内上升流动时,由于受到壁面的加热作用,靠近壁面的流体温度高于截面中心处的流体温度,且近壁流体温度高于拟临界温度,中心流体温度低于拟临界温度。此时,由于SCO2的物性剧烈变化,近壁流体密度减小,受到向上的浮升力作用,流速增大,速度梯度减小,从而降低了流体的切应力;边界层厚度逐渐增加,进而抑制了湍流的扩散,最终导致了传热恶化的发生,即壁面温度出现峰值。由于子通道存在横向交混,横截面冷热流体混合均匀,对壁面发热量的输运能力更强,因此壁面温度峰值低于圆管和环形通道。在热充分发展区,圆管的壁面温度最高,其次是环形通道,子通道壁面温度最低。从壁面温度沿轴向分布可以看出,子通道对流换热效果优于圆管和环形通道。

图6 不同通道壁面温度沿轴向分布

图7 x=44D截面不同通道壁面温度沿周向分布

2.2 平均流动特性

a——圆管;b——环形通道;c——子通道

图9 x=40D横截面上沿径向平均速度分布

2.3 湍流统计量

图10 x=40D横截面上沿径向湍动能分布

图11 x=40D横截面上轴向湍流热流密度沿径向分布

2.4 传热关系式评估

图12所示为计算得到的圆管、环形通道和子通道努塞尔数(Nu)沿轴向的分布,横坐标为无量纲化的加热段长度x/D。从图12中可以看到,在通道前半段,环形通道的Nu与圆管相近,而在后半段环形通道的Nu分布在圆管和子通道之间。子通道Nu始终高于圆管和环形通道,说明子通道换热能力强于相同D的圆管和环形通道。

图12 Nu沿轴向分布

本文选取6个SCO2传热关系式[17-22](表1),将圆管、环形通道和子通道LES计算得到的Nu与传热关系式进行比较,如图13所示。Dittus-Boelter(D-B)公式预测的Nu与LES计算结果以及其他传热关系式的偏差很大,不能准确预测SCO2的传热规律。这是由于D-B公式仅由Re和普朗特数组成,因为SCO2的热物性在拟临界温度附近变化剧烈,需要对传热关系式进行额外修正,以反映浮升力、流动加速以及传热性质变化的影响。表2总结了6种传热关系式的预测准确度,其中参数δ和σ分别为平均相对误差和标准差,定义为:

表1 选择的典型传热关系式

表2 传热关系式预测偏差统计

a——圆管;b——环形通道;c——子通道

(16)

(17)

式中:n为LES计算的数据点个数;δi为第i个相对误差,其计算公式为:

(18)

式中:Nucorr,i为传热关系式计算得到的Nu;NuLES,i为LES计算得到的Nu。

从图13和表2中可以看到,在选择的6种传热关系式中,Kim等[20]提出的传热关系式对3种不同通道的传热预测都具有良好的准确度,平均相对误差仅为5%左右,最大标准差约为11%。Kim等[20]使用q+表征流动加速效应,其关系式相比其他关系式中对流动加速效应做出了修正,且本研究的q+范围与Kim关系式范围接近,这可能是Kim关系式对LES计算结果预测准确度高的原因。Kim关系式是基于竖直圆管实验数据拟合得到的,其对本文3个几何结构的LES计算结果都具有较好的预测效果,这可能与3个几何结构的D相同有关。因此,Kim关系式并不能体现流道几何结构对传热的影响。

3 结论

本文采用LES方法对SCO2在圆管、环形通道和子通道中的流动与传热特性进行了研究,对比了不同流道结构的传热差异,并根据计算结果对典型传热关系式进行了评估,主要结论如下。

1) 在相同热力当量直径、相同工况条件下,子通道的传热性能最好,其次是环形通道和圆管;因此,选取圆管作为研究对象开展SCO2的传热特性研究具有保守性。

2) 子通道的传热性能优于圆管和环形通道,但其近壁面湍动能和湍流热流密度最低,说明不同几何通道湍流变量的强弱与传热性能优劣并无直接对应关系。

3) D-B公式不能准确预测SCO2的传热强弱,而Kim关系式对LES计算得到的圆管、环形通道与子通道的努塞尔数都具有较高的预测精度。

需要指出的是,以上结论是通过对比分析一个混合对流传热工况所得,不一定适用于SCO2的其他传热模式,通道几何结构对SCO2传热特性的影响仍需更全面的深入研究。