空气流动换热数值计算模型适用性分析

2021-02-10卢瑞博宁可为邱志方余霖赵富龙谭思超

哈尔滨工程大学学报 2021年12期

卢瑞博，宁可为，邱志方，余霖，赵富龙，谭思超

(1.哈尔滨工程大学核科学与技术学院，黑龙江哈尔滨 150001； 2.黑龙江省核动力装置性能与设备重点实验室，黑龙江哈尔滨 150001； 3.中国核动力研究设计院,四川成都 610000)

临近空间距离地面20～100 km，近年来因其在天气预报，电磁通信，高分辨成像等方面具有重要应用前景而引起广泛关注[1]。常规动力飞机和卫星无法在临近空间稳定运行，而采用核动力发动机的飞行器由于具有功率密度高、续航时间长、推力大、覆盖范围广等优点，可作为一种临近空间飞行器。该飞行器以反应堆作为热源，高温高速空气为冷却剂兼推进剂[2]。反应堆处于高热流密度，不均匀功率分布的条件，与可压缩空气间的流动换热特性直接影响核动力发动机的安全和运行性能；现有实验条件较难模拟堆内换热工况，因此选择和构建合理的空气流动换热数值计算模型，对预测复杂条件下空气换热特性具有重要意义。学者对圆管内空气流动换热开展了试验和理论研究，在圆管平均换热系数方面，Humble 等[3]研究了内径10.21 mm，有效加热长度609.6 mm的加热管的流动换热特性，分析了普朗特数对流动换热关系式的影响，实验中雷诺数为1万～25万。Lowdermilk等[4]基于Humble的研究，将雷诺数扩展到50万，内壁面温度到1 139 K，得到的换热关系式较前者偏大。Desmon 等[5]采用了耐高温的铂金试验段，扩展内壁面温度到1 696 K，探究了定性温度对换热关系式的影响，研究表明，定性温度的选择可在壁面温度修正和平均膜温修正之间。Deissler[6]通过理论分析得到了普朗特数为0.73条件下的传热关系式，结果表明，合理选择定性温度，物性对换热系数的影响可以忽略。Humble[7]开展的实验中，内壁面温度达到1 694 K，扩展入口温度到833 K，在高热流密度和入口温度条件下，由于空气温度和速度梯度大、物性差异明显，传统的经验关系式适用性较差，需要引入温度修正系数。文献[3-7]考虑了空气的可压缩性对换热的影响，加热管中的空气密度沿流动方向变化较大。由于等截面加热管道的影响，管内流动均为亚音速，管道出口空气流速接近音速。

针对空气流动换热开展的研究中，传热关系式大多符合D-B公式的形式，关系式系数有所差异；已有研究采用不同的定性温度对空气物性进行处理，但无法就处理方法达成一致；实验中采用了测量试验管道进出口储气罐中空气滞止参数的方法测量温度和压力，与管道进出口空气参数相比有一定差异，需要得到管内空气参数的精细化模型；同时，实验中较难模拟反应堆通道加热条件，且高加热功率工况下的试验成本较大。相比于试验和理论研究方法，在误差允许的范围内，通过合理设置边界条件，数值模拟方法能够得到试验管道全场数据，预测所关心工况的流动换热特性，为进一步的细致化试验提供指导，具有诸多优势。因此，在空气流动换热数值计算中，确定合适的湍流模型、优化全场数据的分析处理方法、开展高还原度的仿真试验、提高多工况数值计算计算效率，以达到较精确模拟空气流动换热过程的目的，对于开展复杂换热设备的初步设计具有重要意义。

本文基于核动力发动机反应堆冷却剂通道，拟以单圆管内变物性空气为研究对象，以Fluent为计算平台，确定合理的高温高速可压空气流动换热数值计算方法。

1 空气流动换热计算模型及方法

1.1 物理及物性模型

核动力发动机反应堆中燃料形式为中空六棱柱(如图1所示)，冷却剂通道直径在5.4～5.8 mm[8-9]。在本数值计算中，根据反应堆通道物理尺寸，选择了直径为4、6和8 mm，长度为0.6 m的管道，加热管长径比均大于60，减少长径比对换热的影响。

图1 核动力发动机Fig.1 Schematic of a nuclear engine

采用二维轴旋转模型表征三维圆管模型开展数值计算，上述2种模型的计算结果表明：换热系数的相对误差不超过±0.5%[10]，该表征方式可行。二维轴旋转模型网格绘制方案如图2所示，网格形式为四边形长0.6 m，通过改变边界节点的数量，控制网格数，网格节点增长率小于1.1。

图2 流动换热模型局部网格示意Fig.2 Schematic of the local mesh in the model of flow and heat transfer

本文利用网格数为1.8万、3.8万、9.4万和22.4万的模型，开展无关性验证。结果表明网格数在3.8万以上，速度、温度和对流换热系数等参数的相对误差在±5%以内，满足网格无关性的要求，并且最终本文选择网格数为9.4万。对比了网格长径比对计算结果的影响，结果表明：在文中较大长径比下，能够满足计算精度要求。此外，为了适应壁面函数要求，对近壁面网格加密处理。

压力和温度影响空气密度和物性参数，其中温度对物性参数的影响更大。采用理想气体模型，考虑可压缩效应的影响，通过气体状态方程得到其密度。通过调用REFPROP得到不同温度和压力下空气的物性参数，数据结果表明，可将空气物性参数表示为温度的单值函数，并以分段多项式的形式加载到材料物性库中。文中通过交互命令调用NIST物性库，探究了物性对换热特性的影响，结果表明，调用数据库和拟合多项式得到的换热系数偏差在±2%以内，因此在数据后处理过程中，采用分段多形式计算空气物性参数，能够确保数据后处理过程中物性参数的正确性。

1.2 验证工况范围

高空中飞行器的运行速度为3 Ma，运行温度为-56.5 ℃，此时气流进入飞行器增压段。绝热滞止到速度为0时，空气温度可达606 K，即来流空气总温为606 K，实际进入增压段后空气具有一定的流速。若按照绝热滞止到一定的马赫数，计算得到加热区域空气的静温和流速如表1所示，堆芯前的空气马赫数在0.3以下时，其静温达到572 K以上。因此在计算过程中，确定进口温度范围在450～750 K，进口压力在0.3～1.5 MPa，雷诺数范围在104～4×105。加热采用均匀热流或均匀壁温的方式，加热功率范围在2～5 kW，空气出口温度在550～1 200 K。

表1 不同马赫数对应的空气静温和流速Table 1 The static temperature and velocity of air under different Mach number

1.3 计算方法

文中数值计算采用了质量流量入口、压力出口及无滑移壁面的边界条件，壁面为均匀热流密度加热条件；采用二阶迎风格式，残差收敛至10-5。为探究数值计算模型的适用性，综合考虑适用于加热空气管流的湍流模型。Spalart-Allmaras模型适用于绕流，Reynolds Stress模型适用于复杂三维流场、强旋流，因此计算中不考虑上述2种模型；最终选用Standardk-ε(S-KE); RNGk-ε(RNG-KE); Realizablek-ε(R-KE)及Standardk-ω(S-KW); BSLk-ω(BSL-KW); SSTk-ω(SST-KW)等粘性模型对比分析[11]。

不同壁面函数模型对y+的要求有所区别，且文中雷诺数范围较广，壁面平均y+随流速变化，文中绘制网格考虑了该因素的影响，最终输出所有算例的结果表明，在增强壁面函数条件下，88.3%的工况平均y+在5以下，其余在个别高流速工况下y+值在15以下，可基本满足壁面函数对不同工况y+的影响。文献[7]流动换热关系式的计算方法，基于能量守恒方法计算加热管的平均换热系数和通用传热关系式为：

(1)

Nu=hD/λa=f(Rea0.8Pra0.4)

(2)

式中:T为总温；下标o表示出口；下标i表示入口；S为换热面积；m为质量流量；cp,b为根据流体平均温度计算得到的定压比热容；cp为使用不同定性温度得到的空气的定压比热；λ为导热系数；μ为粘度。

根据流动换热通用传热关系式 (2)，采用不同的定性温度如进出口平均温度Tb和平均膜温Tf，可得到不同的传热关系式，总结不同学者提出的空气流动换热关系式为：

Nu=c×Re0.8Pr0.4=c(md/Aμ)0.8(cpμ/λ)0.4

(3)

Tb=(Ti+To)/2，Tf=(Ts+Tb)/2

(4)

式中：c为努塞尔数表达式的系数，壁面与流体温度修正及长径比修正项可加入其中；m为空气质量流量；d为加热管径；A为加热管流通面积。

在计算开展之前，对数值模拟算法和求解器适用性进行了对比。压力基求解器已经过优化，能够求解可压流。在保证计算精度的前提下，计算可压缩空气问题，压力基求解器拥有更小的收敛步数，和更好的收敛性，因此文中选用了该求解器。对比不同算法的求解精度和收敛步数结果表明，压力基求解器中不同算法对加热空气管流问题求解精度误差在±1.5% 以内，因此可忽略算法对计算结果的影响。由于耦合算法能够同时求解流体3大方程，在计算高速可压流动时更有优势，因此文中较少网格数量的算例采用该算法收敛步数较少，求解过程所需要的时长较短。综合考虑算法计算精度和效率，文中计算采用Coupled耦合求解算法。

2 模型适用性分析

2.1 湍流模型适用性分析

针对选用的不同直径管道，参考实验中空气流量调节方法，在确定边界条件合理性的基础上，逐渐增大空气质量流量，以满足不同雷诺数的要求；改变加热管壁面热流密度，以满足不同加热功率的要求；改变空气进口温度参数，以满足加热管不同入口温度条件的要求。数值计算模型工况的选择满足单一变量准则，不同湍流模型采用相同的工况。通过数值计算方法，得到了不同加热功率、入口温度和压力条件下加热管流的空气参数，为获取管道的平均传热关系式，依据式 (1)～(3) 的数据处理方法，得到了加热管Nu、Pr和Re的关系趋势如图3所示。

图3 不同模型传热关系式对比Fig.3 The comparation of heat transfer correlations with different models

为了验证模型的适用性，文中经计算得到不同模型的传热关系式，并将其与较为成熟的空气换热公式[7](系数c为0.023，适用范围：Re在39万以下，入口温度在297～833 K，壁面温度在1 694 K以下) 和精度较高的Gnielinski公式 (适用范围为：Re为2 300～100万，Pr为0.6～10万，为准确度较高的换热关系式，90%的数据与关系式偏差在±20%以内，大部分在±10%以内) 进行对比，以确定湍流模型的适用性，对比不同模型和经验关系式中的Nu与S-KE中Nu的偏差如图4所示。通过数据后处理可知，不同模型的相同工况下，管内气体的流动参数相近，进出口温差的相对偏差多数在±0.50%，少数在1.7%以内。若仅关心管内空气流动参数，上述模型均能较好描述流体流动状态参数。

图4 不同模型Nu与S-KE模型的偏差Fig.4 The deviation of Nu from S-KE model

由图3可知，上述6种模型下的Nu、Pr和Re的关系均满足传统传热关系式的一般形式，即式 (3)。不同模型能够定性描述空气加热管流换热效果，但不同模型对流体和固体间换热的描述存在差异，固体壁面温度求解不确定。因此，基于传统试验关系式，选择适合的模型准确模拟壁面温度，是解决数值模拟流固耦合问题的关键。将不同模型的Nu与S-KE模型Nu对比如图4所示，Nu主要受到固体壁面与流体的传热温差的影响。几种模型对比结果可知，RNG-KE模型对壁面温度的预测值偏低，而R-KE对壁面温度的预测结果更高，在文中工况范围内，不同模型对壁面温度的计算偏差在±15.0%左右。半数湍流模型中的Nu相对偏高，S-KW、SST-KW和R-KE模型得到的Nu和2个经验关系式偏差均相对较小，初步满足计算精度要求。

综上，数值计算得到的传热关系式与Gnielinski公式和文献[7]中公式的对比表明，在本文工况范围内，采用S-KW、SST-KW和R-KE均能较好描述空气流动换热特性，其中由于R-KE模型偏差在两经验关系式之间，同时该模型可以较好解决粘性底层低雷诺数区域的计算[12]，其适用性更强。基于本节定性分析结果，下文将采用R-KE模型计算空气流动换热工况，壁面处理采用增强壁面函数 (y+值满足要求)，对数值计算模型展开进一步分析。

2.2 定性温度选择

水的物性参数随温度变化较小，在加热管足够长、温差范围小的情况下，充分发展段壁面和流体的温差几乎不随长度变化，可取充分发展段水的物性参数计算管道平均传热系数。空气物性参数随温度变化较大，且在文中高加热功率工况下，加热管道进出口物性参数变化较大，考虑了加热管为普通结构(Type-A)、200 mm非加热入口(Type-B)和圆柱形非加热入口(Type-C)条件下空气的局部换热系数如图5所示。入口形状不能消除入口段效应对空气换热强化的影响，在充分发展段，由于空气加热后膨胀，密度变小，换热能力降低，局部对流换热系数随流向距离增加而减小，与文献[13]的结果相一致。若选取充分发展段局部换热系数作为平均换热系数，会引入误差。因此，对于加热空气管流，需选用合适的定性温度计算平均换热系数，以保证数值计算的后处理精度。本文依据进出口平均温度、平均膜温和壁面温度，此3种定性温度描述数值计算中空气流动换热关系式，其中，为了保证数据覆盖较多工况范围，每组关系式中的数据由120组工况构成。

图5 不同入口段条件下局部换热系数Fig.5 Local heat transfer coefficient under different inlet section conditions

将仿真关系式仿真Nua与文献 [7]中的关系式和Gnielinski公式进行对比如图6所示。其中，基于Gnielinski公式数据，用±10%的误差带表示误差范围。经过数据处理可知，3种定性温度中平均壁面温度数值最高、进出口平均温度最低，平均膜温同时考虑了流体和壁面温度，因此数值在两者之间。由于空气热导率、黏度和比热容，在文中定性温度范围内均为单调递增。空气物性参数影响传热关系式，根据式 (2)、 (3) 可知：在相同换热量的条件下，以壁面温度作为定性温度计算得到的Nu和Re数值均偏小。定性温度对Nu与Re的综合影响效果表现为基于壁面温度计算的传热关系式偏高，基于进出口平均温度的传热关系式偏低，如图6 (d) 所示。在误差方面，对此不同定性温度下的传热关系式与经典关系式，在使用平均壁温作为定性温度计算空气物性时，数值模拟结果仿真Nus与文献[7]和Gnielinski 公式的偏差大于±10%；平均膜温下仿真Nuf偏差略小，部分数据误差超过±10%；而进出口平均温度下仿真Nub与Gnielinski公式误差在±10%以内。因此，在空气流动换热数值模拟中，可优先选用进出口平均温度作为定性温度。

图6 不同定性温度下的传热关系式Fig.6 Conceptual schematic of a nuclear engine

2.3 粘性加热效应影响

流体粘度很大或流动速度很高的情况下需要考虑粘性加热的影响[14]，文中雷诺数范围较大，空气流速范围广，因此在相同的边界条件下，计算了考虑粘性加热和不考虑条件下加热管流的Nu如图7所示。在Re小于20万时，不考虑粘性加热效应条件下的Nu偏大，但较多数据与Gnielinski公式的相对误差仍在10%的范围内；Re大于20万时，由于空气流速较快粘性加热效应明显，两者相对误差超过10%。对于高速空气加热管流，属于热壁面问题，粘性加热作用类似于内热源，表现为降低传热效果[14]。在文中开式堆工况范围内，不同雷诺数条件下，空气流速高，粘性加热效应明显，均须考虑该效应对换热的影响。