圆管内氦氙混合气体与超临界二氧化碳换热特性对比分析

宁可为, 刘凯, 孙汝雷, 赵富龙, 游尔胜, 余霖, 谭思超

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

陆上长时供能是小型核反应堆未来的应用方向之一。出于任务场景对体积、重量的限制,传统压水堆已无法满足要求,新型冷却剂工质得到了广泛的关注,以氦氙混合气体为代表的气体冷却剂与超临界二氧化碳为代表的超临界工质具有压缩性能好、热物性稳定、传热能力优秀等特点,适合与能量转换系统结合,构成运行安全、轻小灵活的核动力系统。

冷却剂换热特性对反应堆运行性能具有重要影响。近年来,诸多学者以空间应用为背景开展了氦氙混合气体物性及流动换热的实验及数值模拟研究,Tournier等[1]基于查普曼-恩斯科格(Chapman-Enskog)输运理论及赫希菲尔德(Hirschfelder)气体混合方法提出了用于预测二元惰性混合气体物性的半经验公式;Vitovsky等[2]针对空间反应堆紧密排布的棒束堆芯结构开展了圆形及准三角通道流动换热实验,研究得到了氦氙混合气体的流动结构,并证明适合采用质量流量平均温度作为定性温度;Dragunov等[3]对不同通道形式下的层流、湍流及变物性条件下的流动换热关系式进行了对比,其推荐采用别祖霍夫(Petukhov)关系式用于氦氙混合气体流动换热计算,该结论与秦浩模拟计算所得相同,同时,Qin等[4]进一步验证发现,未经修正的原始Petukhov关系式计算结果预测精度较高。

超临界二氧化碳作为性能优异的新型冷却剂工质,其跨临界及超临界态下的流动换热特性得到国内外学者的广泛研究。Liu等[5]以TID(tube in duct)型燃料组件为基础,完成了300 MWt超临界二氧化碳反应堆概念设计;刘新新[6]采取实验测量和数值模拟方法研究了超临界二氧化碳在加热直管和螺旋管的换热特性和换热恶化机理;Song等[7]实验对比了超临界二氧化碳在不同管径竖直管内的换热特性,研究发现超临界二氧化碳具有优秀的热力学性质,但在不同的流动条件下,其换热规律差异较大。

综上所述,目前的研究对象基本局限于单一工质,对于不同冷却剂工质换热特性的对比分析不够充分。本文将采用数值计算方法对堆芯圆形冷却剂通道的流动换热特性进行分析,研究不同参数下氦氙混合气体及超临界二氧化碳对堆芯流动换热特性的影响机理,讨论二者作为工质时的反应堆安全特性。

1 几何物理模型及网格划分

1.1 几何模型

陆上多用途的小型核反应堆要求体积高度紧凑,因此,提高燃料比例、组件紧密排布的堆芯设计方案应运而生,六棱柱燃料元件即解决思路之一。中空六棱柱燃料元件以燃料作为基体,圆形冷却剂孔道在基体内部呈三角排布,该设计能够增加燃料体积分数,提供更高的热惯性,有助于抵御异常运行工况下的功率激增。用空气直接冷却的Tory-ⅡC型冲压发动机已对该燃料型式进行了成功实验;Pope[8]及Liu 等[5]采用该型燃料元件进行了超临界二氧化碳反应堆概念设计,结果表明堆芯物理及流动换热特性均满足反应堆运行要求。本文计算采用7冷却剂孔道中空六棱柱燃料元件,截面六边形对边距14.2 mm,冷却剂孔道直径4 mm,燃料元件长0.69 m,具体型式如图1所示。

图1 计算模型及简化Fig.1 Calculation model and simplification

对于该燃料组件型式,由于关注重点在于流体域,因此对模型进行适当简化。Lu等[9]通过计算发现,中空六棱柱型的燃料元件与圆环型的燃料元件在流体域部分的流动参数特性误差最大在1.45%左右,因此对流体域单独建模,得到直径4 mm的圆形冷却剂通道。简化后的模型可以增大网格密度,增加网格的精细度,以提高计算的精确度。以该模型为基础,开展氦氙混合气体与超临界二氧化碳流动换热对比分析。

1.2 物理模型

本文基于计算流体力学(computational fluid dynamics, CFD)软件进行堆芯流动换热数值模拟计算。其中,分离求解器是基于压力的求解器,按顺序逐一的求解每个方程,主要用于不可压流动和低速微可压流动,对于冷却剂最高速度超过0.3 Ma的情况,应当视为不可压流动,采用耦合流模型;而当冷却剂最高速度不超过0.3 Ma时,可视为低速微可压流动,采用分离流模型。

对于氦氙混合气体,以美“普罗米修斯”反应堆为参考,其冷却剂最高速度约40 m/s,因此计算采用分离流模型。Meng[10]以实验值为参考,对多种氦氙混合气体流动的计算模型进行了对比,结果发现相比于SST(shear stress transport) 模型和V2F模型,RKE模型符合程度较好,计算精度较高。k-ε模型通常适用于热传递的工业型应用,并可有效平衡稳定性、计算量和精度。氦氙混合气体在堆芯中的流动同时存在径向与纵向热交换,因此选择k-ε模型。其中标准k-ε模型和可实现k-ε模型均为高雷诺数模型,可实现的k-ε两方程模型可以为全y+壁面处理带来更多灵活性,因此本文计算基于可实现k-ε两方程模型开展。

对于超临界二氧化碳,大量学者已开展直管、螺旋管内的传热特性实验与数值模拟研究,He[11]以实验值和直接数值模拟(DNS)计算结果为参考,提出k-ε或V2F模型具有较高的预测精确性。由于标准的k-ε方程只适用于离开壁面一定距离的高Re数湍流区域,壁面附近低Re数区域须由壁面函数法来修正相关的参数[12]。本文研究中,工质流动的Re基本均大于104,对比其结论得出的工况,选择k-ε模型、全y+壁面处理完成计算。

1.3 网格划分及无关性验证

为确保计算结果精度不受网格影响,且节省计算时间,在计算前首先完成两模型的网格无关性验证。为了将近壁面流体计算更精确,将流体域进行了近壁面网格加密处理,采用了棱柱层网格在流体域近壁面处绘制边界层网格。网格无关性验证结果见图2。由图2可知,对于氦氙混合气体模型,不同网格尺寸对于温度参数的影响微乎其微,但对于对流换热系数的影响较大,低网格数量下,入口段效应体现并不明显,这将对工质换热对比产生不利影响。网格尺寸在2 mm处与0.3 mm处相对误差为1.12%,在可接受范围内。因此,为提高计算效率,网格尺寸确定为2 mm,y+=0.96,网格数量为1 218 434。

表1 网格无关性验证参数及结果Table 1 Grid independence verification parameters and results

图2 计算模型网格无关性验证Fig.2 Verification of grid independence of computational model

2 计算结果及对比分析

2.1 工质物性对比

作为二元混合气体,氦氙混合气体物性与混合比例直接相关。混合气体的热物性采用基于维里系数的稀有气体修正公式计算得出[1]。不同混合比例对堆内通道流动换热具有显著影响,余霖[15]通过研究发现,单冷却剂通道活性区域内,30～40 g/mol为气体流动换热最佳的摩尔质量比例区间,在2 MPa下,31.5 g/mol和40 g/mol氦氙混合气体物性随温度变化如图3所示。可以看出,不同混合比例并不影响物性的变化趋势,其中,较小的混合比例具有较小的密度、较大的定压比热与导热能力。对于动力粘度而言,不同混合比例的影响微乎其微。

图3 氦氙混合气体及超临界二氧化碳物性随温度变化Fig.3 Variation of physical properties of He-Xe mixture and S-CO2 over temperature

作为一元工质,超临界二氧化碳物性仅受温度、压力的影响,其超临界状态使得物性变化趋势与常规流体存在较大差异,该特点主要表现为物性参数在拟临界区附近存在剧烈变化。以10、17 MPa为例,定压比热在临界点附近存在一局部极大值,而在亚临界向超临界跨越的过程中,密度、动力粘度出现陡降,在拟临界区之外,超临界二氧化碳物性几乎不存在显著变化。压力对物性的影响同样出现在临界点附近,跨过拟临界区后,压力提升对物性参数的影响较小。

通过图3中2种冷却剂工质物性对比可以看出,在温度大于400 K后,超临界二氧化碳的密度高于氦氙混合气体,且其具有2倍于氦氙混合气体的定压比热,这使得超临界二氧化碳具有较强的载热能力;氦氙混合气体导热系数高于超临界二氧化碳,但就数值而言,只有约0.15～0.25 W/(m·K),二者导热能力均较差。

2.2 流动换热特性对比

为对比二者流动换热特性,明确不同冷却剂的热工安全性能,针对飞行器闭式反应堆堆芯开展不同运行工况下的反应堆热工安全特性对比。以局部对流换热系数与平均对流换热系数作为主要对比参数,从热流密度、温度、速度等方面对氦氙混合气体与超临界二氧化碳进行对比计算。其中,为了获得相对较优的换热能力,氦氙混合气体相对分子质量选择31.5 g/mol,该值与美“普罗米修斯”方案保持一致;超临界二氧化碳压力的选择需考虑能量转换系统特性[16-17],以17 MPa作为总压进行计算。

2.2.1 加热功率对换热系数的影响

为了研究不同加热功率对换热系数的影响,在入口温度、流速相同的条件下,改变加热功率,得到2种冷却剂管内流动局部对流换热系数。入口温度及流速分别为1 200 K、10 m/s,对加热功率为5、20、80、140 kW/m2下的工况展开计算。其中,氦氙混合气体出口压力2 MPa,超临界二氧化碳取17 MPa,加热方式均取余弦加热。

取点时,沿轴向等距取10个点,导出每点处对应截面的流体平均温度与壁面上的流体温度,二者之差即为该位置处的温度差。计算得到热流密度数据,根据牛顿冷却公式得到传热系数。图4(a)、(b)为不同加热功率下,2种冷却剂工质各自的对流换热系数沿轴向的分布。

图4 不同加热功率下的氦氙混合气体及超临界二氧化碳对流换热系数沿轴向分布Fig.4 Axial distribution of convective heat transfer coefficients of He-Xe mixture and S-CO2 under different heating power

由图4可以看出,对于氦氙混合气体和超临界二氧化碳,局部对流换热系数均呈现入口段先下降后上升、湍流充分发展段基本不变、出口位置下降的趋势,这主要受入口段效应和出口冷却剂物性变化的影响。入口段热边界层存在从零开始生长直至充分发展的过程,后由于边界层中出现湍流,其扰动与混合作用使对流换热系数出现波动,并最终趋于稳定;出口段区域,冷却剂温度升高、压力降低,两者共同作用导致密度加剧降低,冷却剂载热能力相对降低,此时,变物性效应导致冷却剂与固体壁面温差相对增加,对流换热系数降低。

整体而言,加热功率增大会导致对流传热系数先增大后维持,低功率下加热量增加对于对流换热系数的影响较大;当处于较高功率水平后,功率再次提升所带来的影响主要反映在入口段。这主要受加热壁面与流体温差的影响。低功率下的冷却剂携带热量的能力充足,因此,功率增加时,固体仍能够有效地加热流体;功率到达一定程度时,冷却剂无法充分带走壁面热量,导致气体和固体间的温差增大,此时对流传热系数基本维持不变,甚至出现下降趋势。当功率增大到一定值时,由于冷却剂携带热量有限,入口段效应带来的气体温度上升幅度减小,固体温度与流体温差减小,入口段带来的影响效果减弱。

由图4可以看出,对于2种工质,在中间段湍流充分发展区域,对流换热系数均保持基本稳定,不出现较大波动。因此,在该区段取点并作平均,将其作为整段管道的平均对流换热系数。本文取0.2、0.28、0.35、0.4、0.5 m处的局部对流换热系数,取5个点的平均值得到整管平均对流换热系数。

通过将沿轴向位置处的传热系数平均,得到不同加热功率下的平均对流换热系数如图5所示。

图5 氦氙混合气体与超临界二氧化碳对流换热系数随加热功率的变化Fig.5 Variation of convective heat transfer coefficient of He-Xe mixture and S-CO2 with heating power

通过对比可以发现,在5～100 kW/m2的加热功率区间内,超临界二氧化碳对流换热系数始终高于氦氙混合气体,该现象反映了二者换热能力的区别。

在该工况区间下,超临界二氧化碳的换热能力显著高于氦氙混合气体,这主要由流体导热系数、密度、普朗特数等物性参数所决定。

2种工质的对流换热系数对加热功率均不十分敏感,在中低功率水平时,超临界二氧化碳对流换热系数变化较低,随加热功率增加,对流换热系数小幅上升;与之对比,氦氙混合气体的对流换热系数在小加热功率下的增加幅度较大,而在大功率下仍维持在较高的水平。该现象主要反映了2种不同工质载热能力的差别。功率水平较高时,超临界二氧化碳载热的能力不足以带走壁面加热量,故其换热温差相对增加值大于氦氙混合气体。

2.2.2 入口速度对换热系数的影响

为了研究不同入口速度对换热系数的影响,在入口温度、加热功率相同的条件下,改变入口流速。控制壁面恒热流密度为80 kW/m2,定入口温度为1 200 K,入口流速分别为1、5、10、15、20 m/s。其中,考虑到模型在低雷诺数下存在较大的计算误差及堆内实际应用条件,舍去氦氙混合气体1 m/s、超临界二氧化碳20 m/s的工况。

图6(a)、(b)为不同入口速度下,2种冷却剂工质各自的对流换热系数沿轴向的分布。由图6可知,当入口速度变化时,局部对流换热系数分布趋势基本保持一致,即存在入口段、充分发展段和出口段。入口段对流换热系数先下降后上升,充分发展段基本保持不变,而出口段对流换热系数大幅降低。2种冷却剂局部对流换热系数沿轴向分布趋势的差异主要表现在充分发展段,氦氙混合气体在轴向距离中间至后部的位置对流换热系数基本不变或稍有提升,而超临界二氧化碳则随轴向距离增加而始终保持下降趋势。

图6 不同入口速度下的氦氙混合气体及超临界二氧化碳对流换热系数沿轴向分布Fig.6 Axial distribution of convective heat transfer coefficients of He-Xe mixture and S-CO2 under different inlet velocity

不同入口速度下的平均对流换热系数如图7所示。由图7可知,随着入口速度增加,2种冷却剂工质对流换热系数均增大,对流换热效果增强。而在入口速度大于1 m/s后,超临界二氧化碳对流换热系数与其随流速增长的斜率均显著高于氦氙混合气体。从数值上看,当流速从5 m/s增长至20 m/s时,氦氙混合气体对流换热系数仅从620 W/(m2·K)增长至1 013 W/(m2·K),且在15 m/s后,对流换热系数随入口速度的增幅放缓,通过增大流速以改善换热的效果有限;而超临界二氧化碳的对流换热系数随入口流速近似线性增加,在15 m/s时达到3 060 W/(m2·K)。因此,在反应堆设计时,氦氙混合气体应尽量选择较大流速以实现较优的换热效果;而超临界二氧化碳可适当放宽对流速的要求以满足系统运行的要求。

图7 氦氙混合气体与超临界二氧化碳对流换热系数随入口速度的变化Fig.7 Variation of convective heat transfer coefficient of He-Xe mixture and S-CO2 with inlet velocity

2.2.3 入口温度对换热系数的影响

为了研究不同入口温度对换热系数的影响,在入口流速为10 m/s、加热功率为80 kW/m2条件下,改变入口温度,分别为600、900、1 200、1 500 K。

图8(a)、(b)为不同入口速度下,2种冷却剂工质各自的对流换热系数沿轴向的分布。由图8可知,当入口温度发生变化时,氦氙混合气体的局部对流换热系数变化较为剧烈,而超临界二氧化碳趋势与之前相似,仍有较为明显的入口段、充分发展段和出口段。

图8 不同入口温度下的氦氙混合气体及超临界二氧化碳对流换热系数沿轴向分布Fig.8 Axial distribution of convective heat transfer coefficients of He-Xe mixture and S-CO2 under different inlet temperature

不同入口温度下的平均对流换热系数如图9所示。由图9可知,随入口温度的增加,工质对流换热系数呈不同趋势。氦氙混合气体随入口温度增加,管内平均换热系数先增加后下降,在1 000～1 100 K左右有换热系数峰值,约800 K/(m2·k),而超临界二氧化碳管内平均换热系数则随入口温度的增加持续下降,对于所计算工况,在600 K时有换热系数最大值,且在入口温度600～800 K内,换热系数下降速度大于800～1 200 K。该现象出现的原因主要与工质物性随温度变化有关,对于2种冷却剂工质,当入口温度增大,导热系数增加,工质流动的雷诺数降低,而氦氙混合气体普朗特数Pr略有增加,超临界二氧化碳Pr有所减小,综合效果为局部换热系数的降低,但Pr变化的差异导致两者曲线呈现不同形态。

图9 氦氙混合气体与超临界二氧化碳对流换热系数随入口温度的变化Fig.9 Variation of convective heat transfer coefficient of He-Xe mixture and S-CO2 with inlet temperature

两者对比可知,超临界二氧化碳在较低入口温度下的换热效果良好,而氦氙混合气体的换热系数随温度先上升后下降,在1 000 K左右达到峰值;而在数值上,超临界二氧化碳对流换热系数始终高于氦氙混合气体,对流换热能力较强。

通过以上因素的分析可知,对于2种冷却剂工质而言,为了达到较优的堆内换热效果,在设计时应选择高加热功率、大流速,现有超临界二氧化碳反应堆概念设计中[5,18],堆内冷却剂流速取值在5～7 m/s,氦氙混合气体冷却剂流速约20 m/s; 入口温度的选择,氦氙混合气体适合选择在1 000～1 200 K,超临界二氧化碳则适合选择较低的入口温度,现有概念设计在550～650 K左右。

3 对流换热模型修正

尽管以牛顿冷却公式为基础计算出的对流换热系数能够对冷却剂管内流动的对流换热现象做出较好描述,但该参数无法对流动过程做出准确描述,并且基于温差的计算过程使其讨论过程中存在诸多的不确定性。以无量纲参数为基础、根据相似理论发展而来的传热关系式是研究传热问题的重要手段。其典型代表即D-B关系式[19]与Petukhov关系式[20],D-B关系式(dittus-boelter correlation)为:

Nu=0.023ReaPrb

(1)

式中加热流体时a=0.8,b=0.4。

原始型式的别祖霍夫(Petukhov)关系式如下:

(2)

ξ=(1.82lgRe-1.64)-2

K1(ξ)=1+3.4ξ,K2(Pr)=11.7+1.8Pr-1/3

D-B关系式以大量实验值为基础进行拟合,适合于流体普朗特数Pr≥0.6、流体与壁面具有中等温度的换热情况;Petukhov关系式基于常物性结合边界层理论推导而来,能够适应Pr≥0.5及变物性的情况,为简化计算,K1、K2也可分别取常数1.07、12.7,此可视为原始关系式的修正式。

以式(1)～(2)为基础,添加物性修正项的思想能够使之适合不同的实验或数值模拟条件。诸多研究表明,基础物性的差距将导致流动换热特性的差异,因此,有必要对氦氙混合气体及超临界二氧化碳的传热过程开展差异化分析。

对于冷却剂工质,物性变化、流道结构是提出不同型式流动换热关系式的主要考量因素。本研究立足小尺寸圆管内流动,对不同雷诺数下的氦氙混合气体及超临界二氧化碳在定常加热情况下的流动进行广泛计算,将计算所得数据点与现有公式进行比较,得到适合单一工质的较优流动换热关系式。计算过程采用进出口平均温度作为定性温度,在Nu的计算中,对流换热系数计算方法与第2节所述相同。

对于两型冷却剂工质,选择5个具有代表性的经验关系式,开展Nu的计算误差比较。由图10可知,当Re在1×104～5×104的范围内时,Petukhov关系式的预测结果偏差小于20%,而采用常数修正的关系式预测精度相对更高。D-B关系式的一般型式与计算值偏差较大,但其偏离趋势近似为线性,因此考虑对其系数进行调整,使之符合计算结果。图10中,实心点为氦氙混合气体;空心点为超临界二氧化碳。

图10 经验关系式与模拟结果的Nu数比较Fig.10 Comparison of Nu number between correlations and simulation results

对于氦氙混合气体,采用D-B关系式、调整系数进行修正的结果如图11(a)所示,图中实线表示的拟合曲线为:

图11 氦氙混合气体与超临界二氧化碳D-B型公式拟合Fig.11 D-B type formula fitting of He-Xe mixture and S-CO2

Nu=0.34Re0.5Pr0.4

(3)

式中:Re∈[3×103,105];Pr∈[0.15,0.5]。

当Re>104时,修正后曲线与计算数据的符合程度良好,预测精度能够达到95%以上。其原因与流动结构有关,在该计算模型尺寸下,Re在7×103～104内,Nu/Pr0.4随Re变化的斜率出现明显改变,因此,该段可被认为是氦氙混合气体管内流动的转捩区,当Re<7×103时,流动状态为层流,与之对应,Re>104对应湍流的充分发展段。

超临界二氧化碳的D-B关系式修正如图11(b)所示,公式为:

Nu=0.1Re0.65Pr0.4

(4)

式中:Re∈[3×103,1.5×105];Pr>0.5。

需要指出,对于超临界二氧化碳,研究指出浮升力在对流换热过程中具有不可忽略的影响,而重要的判断依据之一即Bo数[21]:

(5)

当Bo满足式(5)时,浮升力的影响即可忽略不计。

本文的计算中,考虑了重力与流动方向相同、相反、垂直3种情况。计算结果表明,不论对于何种工况,Bo≪5.6×10-7。对于系统压力为17 MPa下的4 mm直径管内流动,浮升力对换热过程施加的影响微乎其微。因此,图11(b)中未对3类计算数据加以区分。

对比2种冷却剂的计算数据及修正关系式可知,氦氙混合气体与超临界二氧化碳均可通过修改经验关系式中的系数使得流动符合D-B关系式的一般型式;而在湍流充分发展段,无论是原始型式还是修正型式,Petukhov关系式均能够较好地描述流动换热过程。需要指出,Petukhov关系式的传统应用范围要求Pr>0.5,在计算中可以发现,由于超临界二氧化碳Pr约为0.74,符合应用条件,因此其计算误差相对较小;氦氙混合气体作为低Pr工质,Pr≈0.25,超出了公式设定的应用范围,因此计算误差相对较大,若使用Petukhov关系式则需对Pr做出修正。

4 结论

1)在高加热功率、大入口流速时,氦氙混合气体及超临界二氧化碳具有换热系数极大值,在氦氙混合气体入口温度1 000～1 200 K、超临界二氧化碳600 K时有最大对流换热系数,该参数适合作为反应堆设计时的参考。

2)以D-B关系式的一般型式修正了氦氙混合气体与超临界二氧化碳换热模型,在高Re区域(Re>104),拟合模型与计算值符合良好。