固态熔盐堆又被称为氟盐冷却高温堆[1](fluoride salt cooled high temperature reactor，FHR)，是一种结合了熔盐堆(molten salt reactor，MSR)、高温气冷堆(high temperature gas-cooled reactor，HTGR)和钠冷快堆(sodium fast reactor，SFR)等技术优势的概念堆型。FHR采用和高温气冷堆相同的以石墨为基体、含有三结构同向性型(tri-structural isotropic, TRISO)包覆燃料颗粒的燃料元件，但使用熔盐代替氦气作为一回路冷却剂。由于熔盐具有较高的沸点、较大的热容量和良好的热传导性能，并结合TRISO燃料颗粒的耐高温特性，使得FHR能够在常压、高温、高功率密度条件下运行，有效提高了经济性和安全性。核热泉(nuclear hot spring，NHS)[2]熔盐球床堆是以满功率自然循环为特征的FHR。
核热泉堆的系统如图 1所示。环形的球床堆芯完全浸没在熔盐池内，熔盐冷却剂通过堆芯外侧孔板后横向流过球床反射层，然后进入堆芯，被燃料球加热后，经内孔板进入中心通道，在浮升力的作用下沿着上升通道向上流动，经过熔盐池上部的低阻力中间换热器将热量传给二次侧，最后沿着外侧环形下降通道回到堆芯。
核热泉堆要在满功率下依靠自然循环运行，必须保证整个一回路的总阻力低于上升和下降通道间的浮升力。除了堆芯阻力较大之外，中间换热器的一次侧是主要的阻力来源。降低中间换热器的阻力就要降低流过换热器一次侧的冷却剂流速，但是低流速流体要将反应堆中如此巨大的热量传递给中间回路，除了采用导热率高的冷却剂外，就必须利用巨大的换热面积，于是单位体积换热面积较大的板翅式换热器成为了核热泉堆的首选。图 2所示为这种板翅式超低阻力换热器芯体结构，主要由隔板和翅片堆叠而成，冷热流体流过隔板间翅片形成的通道来交换热量。整个换热器完全浸没在堆顶冷却剂内，一回路冷却剂直接流过换热器，在一回路冷却剂流道两端没有设置任何的进出口管道、联箱或分配器，因而流体流过换热器时只有很小的沿程阻力和摩擦阻力。
|
| 图 2 板翅式换热器芯体结构 |
文[3]依据经验通过简单的流动与传热计算，在满足设定的换热和阻力要求下，初步给定了中间换热器的翅片类型和尺寸，但是这样得到的换热器并不是最佳方案，因此需要对换热器结构进行优化。在板翅式换热器优化过程中经常会遇到不止一个且相互矛盾的目标函数，多目标优化方法可获得满足设计需求的优化解集，为设计提供更多参考，因此在板翅式换热器优化设计中得到了广泛应用。文[4]以换热器总传热率和年度总成本为目标函数，文[5]以最低总成本和最大有效度为目标函数，文[6]以最小熵增和成本最低为目标函数，文[7]以最大有效度和最低压降为目标函数，采用多目标遗传算法(multi-objective genetic algorithm，MOGA)[8]对板翅式换热器进行优化。上述这些优化方法都是采用以空气为工质的经验关系式来计算翅片的传热和阻力特性，而对于熔盐这种高Prandtl数的流体，存在着很大的误差。
本文将采用计算流体力学(computational fluid dynamics，CFD)方法，对不同尺寸结构的翅片通道中的阻力特性和对流换热特性进行数值模拟，并在Ansys Workbench平台上基于响应面方法[9]以翅片阻力因子、Colburn因子和综合评价因子为目标函数，对换热器翅片的高度、翅片间距和翅片厚度进行多目标优化。
1 计算模型及数值方法 1.1 物理模型文[3]中初步设计的中间换热器选取了三角形翅片，本文以此设计为基础，进行优化。图 3为三角形翅片的平面示意图，其中翅片参数包括翅片高度h、翅片间距s和翅片厚度t。针对文[10]中三角形翅片主要参数的适用范围，本文选取的翅片优化范围是：h为4~8.7 mm，s为3~5.5 mm，t为0.05~0.25 mm。
|
| 图 3 三角型翅片示意图 |
在换热器内每一层的翅片都是对称布置的，并且流道通过隔板阻隔相互独立，因此选取如图 4所示的计算单元，这不仅能够大大减少计算量，并且能够有效地反映出换热器的传热特性。
|
| 图 4 模型的计算单元 |
1.2 基本方程及数值方法
针对本文换热器低流动阻力的特性，对计算单元内的流体流动作如下假设：1)工作流体为不可压缩的常物性流体；2)稳态、层流流动；3)忽略流体的体积力和黏性耗散。基于以上假设，可得到控制方程：
连续性方程为
| $ \frac{\partial }{{\partial {x_i}}}\left( {\rho {u_i}} \right) = 0. $ | (1) |
动量方程为
| $ \frac{\partial }{{\partial {x_i}}}\left( {\rho {u_i}{u_j}} \right) =-\frac{{\partial p}}{{\partial {x_i}}} + \frac{\partial }{{\partial {x_i}}}\left( {\mu \frac{{\partial {u_j}}}{{\partial {x_i}}}} \right). $ | (2) |
能量方程为
| $ \frac{\partial }{{\partial {x_i}}}\left( {\rho {u_i}T} \right) = \frac{\partial }{{\partial {x_i}}}\left[{\left( {\frac{\lambda }{{{c_p}}}} \right)\frac{{\partial {u_j}}}{{\partial {x_i}}}} \right]. $ | (3) |
式中：ρ为流体密度，u为流道内速度，下标i和j表示不同方向分量，T为温度，λ为热导率，cp为定压比热容，μ为流体动力黏度。
在计算模型中，边界条件如图 4所示：
1) 入口：ux=uin, uy=uz=0, T=Tin。
2) 垂直于y方向的两个外表面：∂ux/∂y=0, uy=0, ∂uz/∂y=0, ∂T/∂y=0。
3) 垂直于z方向的上下两个表面：T=Tw, ux=uy=uz=0。
4) 翅片和流体的接触表面：ux=uy=uz=0, Tw|s=Tw|f, qw|s=qw|f, -λ(∂T/∂n)w|s=h(Tw-Tf)|f。其中：n为单位法向量，q为热流密度。下标s代表固体，f代表流体，w代表接触表面。下标in、out分别表示进、出口。
本文使用商用的CFD软件Ansys Fluent 16.0进行计算，采用SIMPLE算法对压力和速度进行耦合求解，压力方程采用PRESTO!方法进行离散，动量方程和能量方程都使用二阶迎风格式进行离散。连续性方程的求解残差为10-4，动量方程的求解残差为10-6，能量方程的求解残差为10-10。
1.3 确定换热器特征参数对于不同的翅片，其传热特性和流动阻力特性分别通过Colburn因子j和摩擦因子f来表征：
| $ j = St \cdot \mathit{P}{\mathit{r}^{{\rm{2/3}}}}, $ | (4) |
| $ f = \frac{{\Delta p{D_{\rm{h}}}}}{{2\rho {u^2}L}}. $ | (5) |
其中：
| $ St = \frac{{Nu}}{{\mathit{Re} \cdot \mathit{Pr}}}, $ | (6) |
| $ \mathit{Re} = \frac{{\rho u{D_{\rm{h}}}}}{\mu }, $ | (7) |
| $ \mathit{Pr} = \frac{{\mu {c_p}}}{\lambda }, $ | (8) |
| $ Nu = \frac{{h{D_{\rm{h}}}}}{\lambda }, $ | (9) |
| $ h = \frac{Q}{{A\Delta {T_m}}}, $ | (10) |
| $ Q = m{c_p}\left( {{T_{{\rm{out}}}}-{T_{{\rm{in}}}}} \right), $ | (11) |
| $ \Delta {T_{\rm{m}}} = \left( {{T_{{\rm{out}}}}-{T_{{\rm{in}}}}} \right)/\ln \left( {\frac{{{T_{\rm{W}}}-{T_{{\rm{in}}}}}}{{{T_{\rm{W}}}-{T_{{\rm{out}}}}}}} \right). $ | (12) |
式中：Δp为流道内压降，L为翅片长度，Dh为水力直径，h为对流换热系数，Q换热器交换的热量，ΔTm为对数平均温度，m为质量流量，TW为壁面温度。
1.4 模型的验证选取型号为12.00T翅片[11]的流动和传热的实验结果作为参考值, 表 1列出了该翅片的尺寸。
使用Ansys mesh模块对计算模型进行了网格划分，网格如图 5所示，在靠近翅片的周围采用采用局部网格加密处理。在Reynolds数为500的条件下进行了网格无关性验证，结果见表 2。可见，不同网格数下阻力系数f和Colburn系数j相差非常小，并且当网格数达到801 196以后，f和j变化小于1%。考虑到计算精度和计算效率，选择网格3进行后面的优化计算。
|
| 图 5 计算模型的网格划分 |
| 网格数 | j | f | |
| 实验值 | 0.003 33 | 0.009 8 | |
| 网格1 | 209 880 | 0.003 61 | 0.009 24 |
| 网格2 | 411 492 | 0.003 59 | 0.009 36 |
| 网格3 | 801 996 | 0.003 57 | 0.009 46 |
| 网格4 | 1 625 600 | 0.003 55 | 0.009 50 |
| 网格5 | 3 081 655 | 0.003 54 | 0.009 54 |
为了验证数值计算的可靠性，将数值计算结果与文[11]中三角形翅片的实验值进行对比，如图 6所示。可见，j因子和f因子与实验结果非常吻合，最大的误差不超过10%，同时由于模拟计算过程中忽略了翅片表面的粗糙度而导致数值计算得到的f因子比实验结果偏小，可以认为数值计算的结果是准确的。
|
| 图 6 数值计算与实验结果对比 |
2 三角形平翅片的优化 2.1 优化方法
本文对换热器优化的目的是使中间换热器的流动阻力尽可能低以满足反应堆满功率自然循环的要求，同时体积要尽可能小使整个反应堆系统更加紧凑。在设计时为了权衡阻力和传热，引入了综合评价因子j/f1/3 [12]。在Ansys Workbench平台上，通过拉丁超立方体抽样设计[13]，利用CFD方法求解不同设计点的翅片单元，建立j/f1/3、j因子和f因子的响应面，选用MOGA优化方法，以j因子和j/f1/3最大、f因子最小为目标函数，得到优化的翅片参数。
2.2 翅片尺寸的敏感性通过构建响应面得到了翅片的尺寸参数对其性能的影响。图 7为翅片几何尺寸对换热器性能的敏感性分析。可以发现，在流量固定的情况下，仅翅片厚度t与j因子、f因子和j/f1/3负相关，且影响较小，j因子、f因子和j/f1/3随着翅片高度h和翅片间距s增加而逐渐增大，显然在流量固定时，h和s的增加增大了流通截面，导致通道内的流速降低，Reynolds数减小，因此j因子、f因子和j/f1/3增大，且h的影响比s的稍微大一些。
|
| 图 7 翅片几何尺寸对换热器性能的敏感性 |
2.3 优化结果
以翅片的高度h、间距s、厚度t作为设计变量，以j因子和j/f1/3最大、f因子最小为目标函数，由于自然循环过程对换热器的流动阻力要求较高，因此赋予f因子更高的权重。计算得到3组推荐的设计点见表 3。根据阻力因子最小选择A点作为优化设计点，并对其尺寸进行圆整，得到的优化后的翅片尺寸为：高5.51 mm、间隔2.57 mm、厚0.1 mm。
| A | B | C | |
| h/mm | 5.510 1 | 5.284 4 | 6.256 0 |
| s/mm | 2.569 9 | 2.608 7 | 2.343 7 |
| t/mm | 0.098 8 | 0.067 5 | 0.119 2 |
| j | 0.009 42 | 0.009 41 | 0.009 44 |
| f | 0.032 08 | 0.032 20 | 0.032 19 |
| j/f1/3 | 0.029 65 | 0.029 62 | 0.029 49 |
利用优化后的翅片尺寸重新进行CFD计算，并将结果与初始值进行比较，见表 4。与初始值对比，优化后翅片的f因子减小、j因子和j/f1/3增大，翅片的性能得到了提高。
| 初始值 | 优化值 | 对比/% | |
| j | 0.009 188 | 0.009 416 | 2.48 |
| f | 0.032 40 | 0.032 07 | -1.27 |
| j/f1/3 | 0.028 87 | 0.029 65 | 2.70 |
3 换热器设计
换热器运行参数见表 5。传递2 500 MW的热量总共需要8台换热器，假设换热器的冷热通道满足热平衡，中间回路和一回路都选用熔盐7LiFBeF2(简称氟锂铍，Flibe)作为冷却剂[14]，因其具有较好的慢化能力以及负的空泡反应性系数。各种物性参数见表 6。为了保证冷却剂始终处于熔融状态，二次侧入口温度选取550 ℃，换热器材料为耐高温的碳化硅[15]，换热器的宽度固定为3 m。为了满足自然循环要求，热回路的压降要低于2 kPa，并且尽可能低。
通过商用软件Aspen的换热器设计组件，对中间换热器进行设计。优化前后的换热器尺寸见表 7。可以看到，优化后的换热器比初始设计体积减小了30%，使一回路系统更加紧凑。
| 流体 | 温度 | 最高压降 | 流量 | ||||
| ℃ | kPa | kg·s-1 | |||||
| 热通道 | 入口 出口 |
Flibe | 1 000 800 |
2 | 646.8 | ||
| 冷通道 | 入口 出口 |
Flibe | 550 950 |
2 | 323.4 |
4 结论
中间换热器是自然循环熔盐球床堆的关键部件之一，在设计时为了权衡阻力和传热引入了综合评价因子j/f1/3，并选定j因子、f因子和j/f1/3作为评价因子，通过基于CFD方法的结构优化平台构建响应面来确定翅片结构参数对这3个评价因子的影响。在流量固定的情况下，翅片高度h和间距s对j因子、f因子和j/f1/3为正影响，且h的影响比s的稍微大一些，翅片厚度t为负影响且影响较小。采用MOGA优化方法，对翅片进行了优化，得到最佳的换热器单元尺寸：翅片高5.51 mm、间隔2.57 mm、厚0.1 mm。优化后的翅片与初始选取的翅片相比，减小了阻力但提高了换热性能。根据此翅片尺寸进行了换热器优化设计，换热器优化后比未优化前的设计体积减小了30%，优化效果显著。此方法可以在类似设计中适当推广。
| [1] | FORSBERG C W, PICKARD P S, PETERSON P F. Molten salt cooled advanced high temperature reactor for production of hydrogen and electricity[J]. Nuclear Technology, 2003, 144: 289–302. DOI:10.13182/NT03-1 |
| [2] |
吕应中. 在任何功率下长期自动运行生产高温核能的方法: 201010145086. 9. 2010-08-25. LÜ Y Z. A progress for producing high temperature nuclear energy with full power natural circulation operation: 201010145086. 9. 2010-08-25. (in Chinese) |
| [3] |
吕应中.
大规模替代化石及其他有限能源的固有安全高温核动力[J]. 核科学与工程, 2011, 31(1): 1–8.
LÜ Y Z. An inherently-safe high-temperature nuclear energy producing process for the replacement of the fossil and other depletive energy on a large scale[J]. Chinese Journal of Nuclear Science and Engineering, 2011, 31(1): 1–8. (in Chinese) |
| [4] | NAJAFI H, NAJAFI B, HOSEINPOORI P. Energy and cost optimization of a plate and fin heat exchanger using genetic algorithm[J]. Applied Thermal Engineering, 2011, 31: 1839–1847. DOI:10.1016/j.applthermaleng.2011.02.031 |
| [5] | SANAYE S, HAJABDOLLAHI H. Thermal-economic multi-objective optimization of plate fin heat exchanger using genetic algorithm[J]. Applied Energy, 2009, 87: 1893–1902. |
| [6] | AHMADI P, HAJABDOLLAHI H, DINCER I. Cost and entropy generation minimization of a cross-flow plate fin heat exchanger using multi-objective genetic algorithm[J]. Journal of Heat Transfer, 2011, 133(2): 21801–21810. DOI:10.1115/1.4002599 |
| [7] | HAJABDOLLAHI H, TAHANI M, SHOJAEEFARD M H. CFD modeling and multi-objective optimization of compact heat exchanger using CAN method[J]. Applied Thermal Engineering, 2011, 31: 2597–2604. DOI:10.1016/j.applthermaleng.2011.04.027 |
| [8] | ABDULLAH K, DAVID W C, ALICE E S. Multi-objective optimization using genetic algorithms:A tutorial[J]. Reliability Engineering & System Safety, 2006, 91(9): 992–1007. |
| [9] | MYERS R H, MONTGOMERY D C. Response surface methodology:Process and product in optimization using designed experiments[M]. New York, USA: John Wiley & Sons, 1995. |
| [10] | SHAH R K, DUSAN P, SEKULIC D P. Fundamentals of heat exchanger design[M]. New York, USA: John Wiley & Sons, 2003. |
| [11] | KAYS W M, LONDON A L. Compact heat exchanger[M]. 3rd ed, New York, USA: McGraw-Hill, 1984. |
| [12] | BHOWMIK H, LEE K S. Analysis of heat transfer and pressure drop characteristics in an offset strip fin heat exchanger[J]. International Communications in Heat and Mass Transfer, 2009, 36(3): 259–263. DOI:10.1016/j.icheatmasstransfer.2008.11.001 |
| [13] | IMAN R L, CONOVER W J. Small sample sensitivity analysis techniques for computer models, with an application to risk assessment communications in statistics:Part A[J]. Theory and Methods, 1980, 17: 1749–1842. |
| [14] | SOHAL M S, EBNER M A, SABHARWALL P, et al. Engineering database of liquid salt thermos physical and thermochemical properties[R]. Idaho Falls, USA: Idaho National Laboratory, 2010. http://digital.library.unt.edu/ark:/67531/metadc840052/ |
| [15] | SOMMERS A, WANG Q, HAN X, et al. Ceramics and ceramic matrix composites for heat exchangers in advanced thermal systems:A review[J]. Applied Thermal Engineering, 2010, 30: 1277–1291. DOI:10.1016/j.applthermaleng.2010.02.018 |