基于雷诺时均模型的高温空气流动传热数值研究

焦守华，王金雨，曾 未，孙 畅，李文强，杨 婷，柴 翔

(1.上海交通大学 核科学与工程学院，上海 200240；2.中国核动力研究设计院 核反应堆系统设计技术重点实验室，四川 成都 610213；3.上海电力大学 自动化工程学院，上海 200090)

持续稳定的动力供应将极大提高飞行器的作战能力和续航能力。使用化石燃料作为动力源，会对飞行器的飞行速度、高度和范围产生一定的限制。核能作为一种新型能源，与传统的化石燃料相比具有能量密度高、体积占比小、清洁无污染等一系列优点。将核动力应用到这些飞行器，将极大提升它们的作战能力和续航能力[1]。20世纪，美苏在核动力发动机方面的技术发展路线上都不约而同地指向了核涡轮喷气发动机。核涡轮发动机反应堆其堆芯通常采用高温空气作为冷却剂，高温空气的流动传热特性直接关系到涡轮发动机的运行性能。因此数值研究反应堆内高温空气流动换热特性，可以为设计高温空气流动传热实验和核涡轮发动机反应堆提供一些参考和指导。

在气体流动传热方面，国内外已经开展了很多研究， Locke J M和Landrum D B[2]对再生冷却通道中的超临界氢气的流动传热特性进行了相关研究，并评估了相关的对流传热关系式对氢气的适用性。Ehsan M M等[3]对处在加热或者冷却条件下的超临界二氧化碳的流动传热和压降特性以及相关的对流换热关系式进行了全面综述。Yangyu Guo和Moran Wan[4]系统地研究了微尺度气体流动和传热的第二定律，阐明了在微观尺度热和流体流动中熵产生的物理性质和数学表达式。卢芳等[5]采用FLUENT软件对板式换热器通道内的气体流动特性进行了数值计算。上述学者多是对非高温空气的流动传热特性研究，涉及高温空气的研究比较少。

本文使用开源CFD软件OpenFOAM,对反应堆冷却剂通道内高温空气的流动传热进行数值模拟，探究了直径为4.24 mm、长度为760 mm的圆形通道内均匀加热和非均匀加热情况下高温空气的流动传热特性。研究结果对于设计高温空气流动传热实验具有一定指导意义。

1 数值方法简介

1.1 控制方程

本文研究的所有流动和传热过程都受三大基本定律支配，即质量守恒定律、动量守恒定律和能量守恒定律。在三维空间连续介质内，牛顿黏性流体遵循的三大控制方程为[6]：

质量方程为：

(1)

动量方程为：

(2)

(3)

(4)

能量方程为：

(5)

式中：t——时间；

x,y,z——空间X、Y、Z方向的位置；

u,v,w——流体流速在空间X、Y、Z方向上的分量；

p——作用在流体微元上面的压力；

Fi——流体微元所受体积力沿空间X、Y、Z方向上的分量；

τij——表面黏性应力沿空间X、Y、Z方向上的分量；

Cp——流体的比定压热容；

T——流体的温度；

λ——流体的导热率；

ST——黏性耗散项。

1.2 空气物性参数

本文选用理想气体模型来研究高温空气在圆形通道内的流动传热。在XNJ-140E核涡轮发动机反应堆中[7]，冷却剂的温度在570～1 250 K范围内，压强在0.1～0.5 MPa范围内。在相同的温度下，压强为0.1 MPa和压强为0.5 MPa对应的空气动力黏度、比定压热容和导热系数的相对误差均在0.1%范围内，因此压强对空气的上述物性影响可以忽略不计，只需考虑温度对上述物性的影响。通过NIST数据库[8]得到不同温度下空气物性参数，然后进行多项式拟合，得到的空气物性随温度变化关系式如下所示：

(6)

λ=-4.676×10-15×T4+

2.257×10-11×T3-4.617×10-8×T2+

9.25×10-5×T+0.002 17

(7)

(8)

2 计算过程

2.1 几何模型

本文以XNJ-140E核涡轮发动机反应堆[7]的单个冷却剂通道作为计算对象。该对象是直径D为4.24 mm，长度L为760 mm的圆柱管。为了在加热段进口得到充分发展的流动，在入口处增加了一长为200 mm的辅助进口段；为了降低出口边界条件对加热段结果的影响，在加热段出口增加了一长为40 mm的延伸出口段，因此整个几何长度为1 000 mm。本文研究的对象几何结构如图1所示，x方向为流体的流动方向，y方向为径向。为了减少计算量，将网格创建为二维轴对称。

图1 几何模型Fig.1 Geometric model

2.2 边界条件

计算模型的边界条件设置如下：

(1)进口边界：给定质量流量和温度;

(2)出口边界：给定压力;

(3)非加热壁面边界：无滑移和绝热;

(4)加热壁面边界：给定热流密度。

根据边界条件的不同，本文的计算工况如表1所示，表1中的热流密度为管道的平均热流密度。其中，工况1、工况2、工况3用来研究均匀加热下热流密度的影响，工况2、工况4、工况5用来研究均匀加热下进口温度的影响，工况2、工况6、工况7用来研究均匀加热下进口雷诺数的影响。工况8-13是非均匀加热工况，工况8、工况9、工况10是余弦功率分布，归一化功率最大值为1.57，工况11、工况12、工况13的功率沿轴向先不变后线性降低，这和美国研究人员建成的Tory试验堆的轴向功率分布是一致的，本文将这种功率分布称之为Tory型功率分布。

表1 计算工况Table 1 Case for simulation

2.3 网格敏感性分析

为了在有限的计算资源下，更高效地得到合理精确的计算结果，本文首先使用三套网格进行预计算。预计算设置相同的边界条件，仅在网格数量上有所不同。并取加热段出口的速度及温度分布情况和加热壁面温度分布情况作为网格敏感性分析的判据。具体的网格划分方案如表2所示。本文计算工况的Y+平均值为0.773，符合湍流模型适用范围。

表2 不同的网格划分方案Table 2 Different meshing schemes

不同网格数目下，工况2计算结果如图2至图4所示。由图可知，三种网格数量计算的温度和速度分布相差很小，相对误差在1%范围内，因此选用网格3的网格数目作为最终计算的网格标准。

图2 加热段出口温度分布图Fig.2 Temperature distribution of outlet

图3 加热段出口速度分布图Fig.3 Velocity distribution of outlet

图4 加热壁面温度分布图Fig.4 Temperature distribution in heated wall

2.4 湍流模型选取

在各湍流数值模型中，k-ε模型能较好地模拟充分发展的湍流流动，k-ω模型则能较好地模拟各压力梯度下近壁面处的流动情况。Menter提出k-ωSST模型，不仅在近壁面处保留原始的k-ω模型，也对充分发展的湍流流动状况运用k-ε模型[9]。k方程和ω方程分别为：

(9)

(10)

(11)

(12)

(13)

式中：t——时间；

ρ——流体密度；

Gk——湍动能k的产生项；

Gω——比耗散率ω的产生项；

U——流体速度；

μ——动力黏度；

x——三维坐标系；

μt——流体的比定压热容,μt=ρkT;

T——湍流时间尺度；

S——流体平均应力张量的模；

y——距壁面的距离；

ν——流体的运动黏度。

σk、σω、β和β*可通过混合函数求解得到；γ为混合系数；CDkw与流体的横向扩散项相关。

针对圆管内均匀加热的气体流动传热， Taylor等[10]总结氢气、氦气、氮气和空气在均匀加热条件下的实验数据，提出的均匀加热条件下的气体换热关系式，如公式(14)所示，关系式计算值与实验值偏差在25%以内。本文选用k-ωSST湍流模型对工况2进行预计算，计算得到沿流动方向上不同位置处的努塞尔数，并将模拟结果与Taylor关系式计算结果进行对比，结果如图5所示。

图5 工况2对应的Nu随x/D变化Fig.5 Nu vs x/D for case2

(14)

式中：Tb——沿流动方向某一位置处的流体主流温度；

Ts——与Tb相对应的壁面温度。

由图5可知，模拟计算值与由Taylor关系式计算的值在x<20D吻合很好，误差极低，在x>20D趋势一样，误差略微增大，但误差总体还是在10%以内。因此本文选用k-ωSST湍流模型进行均匀和非均匀加热下的数值计算。

3 计算结果及分析

3.1 均匀加热情况下结果分析

图6 均匀加热情况下模拟值与Taylor关系式值对比Fig.6 Comparison of Nu obtained by simulation against Nu calculated by Taylor correlation under uniform heating

选取工况1、工况2和工况3分析热流密度对对流换热的影响，得到的结果如图7所示。由图7可知，沿流动方向，努塞尔数一直在减小，即对流换热能力越来越弱。此外，随着热流密度的增加，努塞尔数在减小，但幅度并不大。

图7 不同热流密度下努塞尔数沿轴向分布Fig.7 Profiles of Nu along the axial direction at different heat flux

选取工况2、工况4和工况5分析进口温度对对流换热的影响，计算结果如图8和图9所示。由图8可知，沿流动方向，努塞尔数一直在减小，对流换热能力越来越弱。在进口雷诺数一定的情况下，随着进口温度增加，努塞尔数会略微降低，温度升高40 K，平均努塞尔数仅降低2%。

图8 不同进口温度下努塞尔数沿轴向分布Fig.8 Profiles of Nu along the axial direction at different inlet temperature

选取工况2、工况6和工况7来探究进口雷诺数对对流换热的影响。计算结果如图9所示。由图9可知，沿流动方向，努塞尔数一直在减小，进口速度越大，对应的努塞尔数越高，进口雷诺数每增加7 500，平均努塞尔数会增加30%左右。

图9 不同进口雷诺数下努塞尔数沿轴向分布Fig.9 Profiles of Nu along the axial direction at different Re of inlet

3.2 非均匀加热情况下结果分析

选取工况1、工况2、工况3和工况8、工况9、工况10来研究均匀加热和余弦加热对对流换热的影响，计算结果如图10和图11所示。由图10可知，在余弦功率加热下，沿流动方向，由于进口效应，对流换热系数会先下降，之后由于热流密度增加，对流换热系数会上升。由图11可知，与均匀加热一样，沿流动方向，努塞尔数也是一直在下降。在x<60D，均匀加热下的努塞尔数要低于余弦加热的，之后则相反。冷却剂在入口段后形成了充分发展的流动。均匀加热下进口段处加热功率较高，近壁面处冷却剂的密度和流速受温度的变化而剧烈变化；在余弦加热下，进口段处的加热功率较低，湍流边界层随着加热功率的增加而逐渐演变。

图10 对流换热系数沿轴向分布Fig.10 Profiles of convective heat transfer coefficient along the axial direction

图11 努塞尔数沿轴向分布Fig.11 Profiles of Nu along the axial direction

将所有余弦加热工况的模拟值与Taylor关系式计算值相比较，结果如图12所示。由图12可知，模拟计算值与Taylor关系式计算值偏差在10

图12 余弦加热情况下模拟值与Taylor关系式值对比Fig.12 Comparison of Nu obtained by simulation against Nu calculated by Taylor correlation under cos-type heating

本文以工况2和工况9为例分析了在靠近进出口处和中间部分的流体径向温度分布情况，结果如下所示。由图13可知，在x=10D处，均匀加热下，流体径向温度分布很不均匀，壁面温度是主流温度的1.52倍，数值模拟值与Taylor关系式计算值偏差0.53%；余弦加热下，该处由于热流密度较低，流体径向温度分布相对比较均匀，壁面温度是主流温度的1.125倍，数值模拟值与Taylor关系式计算值偏差7.54%。均匀加热和余弦功率加热两种加热模式在x=90D处的热流密度都较高，分别是277 155 W/m2和433 133 W/m2。由图14可知，由于较高的热流密度，壁面温度分别是主流温度的1.33和1.55倍。努塞尔数的数值模拟值与Taylor关系式计算值偏差分别是4.57%和2.86%。由图15可知，在x=170D处，均匀加热下，流体径向温度分布很不均匀，壁面温度是主流温度的1.24倍，数值模拟值与Taylor关系式计算值偏差5.12%；余弦加热下，热流密度较低，仅为均匀加热下的1/4，流体径向温度分布相对比较均匀，壁面温度是主流温度的1.07倍，数值模拟值与Taylor关系式计算值偏差27.27%。由于Taylor关系式是针对均匀加热的实验数据总结而出，根据主流温度确定的物性校正出来的，实验中由于热流密度较高，流体径向温度分布不均匀情况比较明显，因此该关系式较适用于高热流密度的加热区域。但是在余弦加热情况下，在进出口处，加热功率较低，流体径向温度分布相对比较均匀，因此数值模拟值与Taylor关系式计算值偏差比较大。

图13 x=10 D处温度分布图Fig.13 Temperature distribution at x=10 D

图14 x=90 D处温度分布图Fig.14 Temperature distribution at x=90 D

图15 x=170 D处温度分布图Fig.15 Temperature distribution at x=170 D

将所有Tory型加热工况的模拟值与Taylor关系式值相比较，结果如图16和图17所示。Tory型加热情况下，功率沿轴向先不变后线性降低，这和美国研究人员建成的Tory试验堆的轴向功率分布是一致的。

图16 Tory型加热情况下模拟值与Taylor关系式值对比Fig.16 Comparison of Nu obtained by simulation against Nu calculated by Taylor correlation under Tory-type heating

图17 Tory型加热情况下模拟值与Taylor关系式值偏差Fig.17 Error of Nu obtained by simulation against Nu calculated by Taylor correlation under Tory-type heating

由图16和图17可知，在189区域，功率线性降低，属于非均匀加热情况，模拟计算值与Taylor关系式计算值的偏差沿流动方向逐渐增大，在出口处达到了最大值，最大值为15%。

3 结论

本文基于开源CFD软件OpenFOAM，对直径为4.24 mm、长度为760 mm的圆形通道内均匀加热和非均匀加热情况下高温空气的流动传热进行了研究，可以得到如下结论。

(1)均匀加热情况下，k-ωSST湍流模型数值模拟结果与Taylor关系式计算结果吻合较好，趋势一致，k-ωSST湍流模型适用于均匀加热下的高温空气流动传热模拟。

(2)均匀加热情况下，沿流动方向，努塞尔数一直在降低，对流换热能力越来越弱。进口雷诺数对努塞尔数影响比较大，进口雷诺数每增加7 500，努塞尔数会增加30%。热流密度和进口温度对其影响较小。

(3)余弦加热情况下，k-ωSST湍流模型数值模拟结果与Taylor关系式计算结果在进出口偏差比较大。

(4)由于Taylor关系式是针对均匀加热的实验数据总结而出，根据主流温度确定的物性校正出来的，实验中由于热流密度较高，流体径向温度分布不均匀情况比较明显，因此该关系式较适用于高热流密度的加热区域。但是在余弦功率加热情况下，在进出口处，加热功率较低，流体径向温度分布相对比较均匀，因此数值模拟值与Taylor关系式计算值偏差比较大。

(5)Tory型加热情况下，沿流动方向，k-ωSST湍流模型模拟结果与Taylor关系式计算结果在前半段偏差比较小，在后半段偏差逐渐增大，并在出口处达到了最大值。