二维方腔内HITEC熔盐熔化非稳态传热数值模拟

Unsteady Heat Transfer of HITEC Molten Salt in Two-Dimensional Square Cavity

摘要: 本文以混合熔盐HITEC为研究对象，在内置圆形加热棒的封闭方腔内，考虑液相自然对流的情况下，建立熔盐熔化传热的数学模型。FLUENT软件被用来模拟传热过程。考虑不同功率、不同加热棒尺寸对熔盐进行加热，分别获得不同情况下熔盐熔化需要的时间。并分析在某一加热功率下，方腔内温度分布。数值模拟结果表明，自然对流不能忽视，它促进方腔上部的传热，缩短了熔盐熔化需要的时间。随着液相区不断扩大，涡流的影响范围也不断扩大。熔盐熔化过程经过三个阶段：纯导热阶段，相变阶段，自然对流主导阶段，不同阶段方腔内的温度分布差别明显。另一方面，随着加热功率增大，或者增大加热棒尺寸，液相比愈来愈呈现随时间非线性变化的规律。

Abstract: A numerical study has been conducted on the heat transfer of the molten HITEC salt, in which the natural convection heat transfer of the liquid salt was considered in the closed square cavity equipped with a heating bar. The two-dimensional heat transfer model of the molten salt was es-tablished and the problem was numerically solved with the FLUENT software. The dynamic change of the molten region as well as the melting fraction in the cavity was addressed in the condition of different heating power. The numerical result has revealed that it is considered mainly caused by the buoyancy-driven fluid flow of the molten salt and the vortex would enlarge with the molten range increasing. As the result of the coupled heat transfer of the natural convection and the conduction, the process of the salt melting can be divided into three stages, i.e. the heat conduction dominated stage at the beginning, the phase changing stage in the middle period, and the convection dominated stage in the last. The temperature distribution in the square cavity is obviously different in different stage. Either increasing the heating power or enlarge the heating bar size are all helpful to the salt melting process. Nevertheless, changing the heating bar size seems to be more effective for the melting dynamics in comparison.

1. 引言

由于化石燃料等不可再生资源的短缺，人们日益重视能源的使用朝多元化、清洁化、可持续的方向发展。根据国际能源署(IEA)发布的2016年世界能源展望，未来20多年里世界能源需求预期增长30% [1] 。当前我国实现全面建成小康社会和现代化目标，能源消费将持续增长。人类面临能源革命机遇，掀起了一系列的创新的技术，其中储能技术在能源的供应和能源可持续发展中发挥重要作用，普遍应用于太阳能高温储能、电力负荷转移、工业余热回收。发展高效储能技术，有利于提高能源的使用效率 [2] [3] 。

在储能技术发展中，熔盐作为高温载热蓄热流体在冶金领域、化学产业和金属处理行业已有一定的应用。从20世纪80年代，随着太阳能热发电技术的不断发展，熔盐蓄热储能的研究不断获得进展 [4] 。优化换热单元的几何结构是一种提高蓄热系统传热性能的重要改进方式，文献 [5] 认为，矩形结构材料熔化时间少于圆形结构的时间。

同时蓄热材料的自然对流对蓄热的影响也不容忽视。文献 [6] 给出自然对流的存在加强方腔蓄热装置相变材料的内部传热，因此促进材料熔化。文献 [7] 分析在充满空气的方腔内置圆柱均匀加热对自然对流的影响。不同相变材料由于物理性质不同，相变熔化过程也不同。文献 [8] 、 [9] 分别研究多孔介质、纯金属镓在封闭方腔内相变熔化蓄热过程。方腔壁面初始条件和方腔内部热源的有无、分布、大小等都会对自然对流造成影响 [10] [11] [12] 。但是文献 [10] [11] [12] 只研究了空气的自然对流过程，没有涉及方腔内相变材料的固液相转变的熔化过程。综上文献所述，前人研究相变蓄热材料没有涉及在矩形结构中熔盐材料熔化蓄热问题，更没有考虑熔盐熔化过程中液相的自然对流对蓄热过程的影响。

本文的研究对象是一种以混合熔盐HITEC作为相变蓄热材料的方腔式蓄热结构，在方腔内布置加热源，且四周绝热。利用FLUENT软件，研究相变材料从固相到液相的熔化相变传热过程，液相的自然对流对相变蓄热过程温度变化的影响，以及考虑内置热源功率的变化和加热棒尺寸的变化等因素对熔盐相变过程的影响。采用FLUENT软件求解封闭方腔内自然对流问题在文献 [13] 中已有报道。

2. 物理模型

图1所示，在正方形方腔内有三元混合熔盐HITEC，金属棒作电加热元件对熔盐进行加热，方腔四周绝热的条件下，用数值模拟的方法研究熔盐熔化过程的变化情况。方腔的边长a = 30 cm，金属棒直径d = 10 cm。熔盐的初始温度为25℃(298 K)，金属棒的电加热功率为q = 500~3000 W/m²。

从传热学方面看，熔盐熔化过程固液相问题的实质是相变的传热问题，又被称作斯蒂芬(Stefan)问题 [14] ，对于混合熔盐相变材料，在固、液共存的两相间有一个分界面，它是具有一定宽度的两相糊状区，相界面随着熔盐熔化或凝固的进行，伴有潜热吸收或者释放。熔盐相变时液相会产生自然对流，这是因为温度变化造成的。

为了分析方便，对物理模型作以下假设：

1) 熔盐不可压缩，各向同性，与外界环境不发生任何物质交换；

2) 在材料熔化区间，熔盐的物性参数采用温度的函数，其具体形式在下文叙述；

3) 液相的流动为层流；

4) 熔盐内部初始温度分布均匀。

3. 数值模型

FLUENT软件中凝固熔化模型(Solidification/Melting Model)是由Voller和Prakash等提出，它被称为“计算材料凝固熔化的焓–多孔法” [15] [16] 。在该方法中，熔融界面并不被跟踪，而是在数值分析相变时，计算区域采用材料的多孔性进行表征，其特征在于引入液相系数β，其定义为

$\beta =\{\begin{array}{ll}0\hfill & T<T_s\hfill \\ \frac{T-T_s}{T_l-T_s}\hfill & T_s<T<T_l\hfill \\ 1\hfill & T_l<T\hfill \end{array}$ (1)

Figure 1. Geometric model of square cavity in the simulation

图1. 问题的方腔模型示意图

其中，T_s为固相熔化温度，T_l为液相凝固温度。

3.1. 数学描述

1) 连续性方程

$\frac{\partial \rho }{\partial t}+\nabla \cdot \left(\rho V\right)=0$ (2)

其中，ρ为密度， $V$ 为流体速度。

2) 动量方程

$\frac{\partial \left(\rho u\right)}{\partial t}+\nabla \cdot \left(\rho uV\right)=\nabla \cdot \left(\mu \nabla u\right)+S_u$ (3)

$\frac{\partial \left(\rho v\right)}{\partial t}+\nabla \cdot \left(\rho vV\right)=\nabla \cdot \left(\mu \nabla v\right)+S_v$ (4)

其中，u、v分别为速度 $V$ 在x、y方向上的速度分量；S_u、S_v分别为 $V$ 取u、v时动量方程源项，它们与液相系数β有关。

由于在相变过程中只有液态熔盐才可以流动，根据焓–多孔理论，动量的损失是因为两相区孔隙率减少形成的，则x方向动量方程源项为：

$S_u=\frac{{\left(1-\beta \right)}^2}{{\beta }^3+\epsilon }{A}_{mush}u-\frac{\partial p}{\partial x}$ (5)

其中，β液相系数，ε为小于0.001的一个常数，避免被0除。A_mush 为糊状区的一个连续数，通常取到10⁵。

在计算浮升力引起的流动时，采用Boussinesq假设，则y方向动量方程源项为：

$S_v=\frac{{\left(1-\beta \right)}^2}{{\beta }^3+\epsilon }{A}_{mush}v-\frac{\partial p}{\partial y}-{\rho }_{ref}g\left[1-\alpha \left(T-T_{ref}\right)\right]$ (6)

其中，T_ref表示参考温度，ρ_ref为参考温度时的密度，α为体膨胀系数，g为重力加速度。

3) 能量方程

凝固熔化问题能量方程可表示为：

$\frac{\partial }{\partial t}\left(\rho H\right)+\nabla \cdot \left(\rho VH\right)=\nabla \cdot \left(k\nabla T\right)$ (7)

其中，H为物质的焓，由材料的显热h和潜热ΔH组成：

$H=h+\Delta H$ (8)

$h=h_{ref}+{\displaystyle {\int }_{T_{ref}}^T{C}_p\text{d}T}$ (9)

其中，h_ref为参考焓，T_ref表示参考温度，C_p表示等压比热容。潜热部分可从材料的熔化潜热值L求得：

$\Delta H=\beta L$ (10)

3.2. 初始条件与边界条件

本研究的计算区域为金属棒外壁与方腔之间的夹层部分，并不包括金属棒本身。为了便于计算，本文简化了模型的初始条件与边界条件。整个区域由圆形金属棒加热。

初始条件： $t=0$ ， $u=v=0$ ， $T=298\text{K}$ ；

边界条件： $x=\pm 0.15\text{m}$ ， $\frac{\partial T}{\partial x}=0$ ； $y=\pm 0.15\text{m}$ ， $\frac{\partial T}{\partial y}=0$ ；

$r=0.0\text{5m}$ 的圆表面上，

$-k\frac{\partial T}{\partial r}=q$ (11)

其中，q为电加热棒功率。

3.3. 熔盐HITEC的物性参数

混合硝酸盐是一种低熔点的水溶无机盐混合物，它是由硝酸钠、亚硝酸钠和硝酸钾按照质量比7:40:53配制而成。它的熔点415 K，沸点953 K。

熔盐处于液态时物性参数拟合公式 [17]

密度 ρ /(kg/m³)：

$\rho =2287.7993-0.748T$ (12)

定压比热容C_p/(J/( kg·K))：

$C_p=1507-0.1T$ (13)

导热系数l /(W/( m·K))：

$\lambda =0.8099-7.827\times {10}^{-4}T+1.043\times {10}^{-7}{T}^2$ (14)

动力粘度 η/(Pa·s)：

$\eta =0.06248-2.253\times {10}^{-4}T+2.7903\times {10}^{-7}{T}^2-1.1738\times {10}^{-10}{T}^3$ (15)

根据熔盐的凝固温度和熔化温度，确定 $T_s=\text{415}\text{K}$ ， $T_l=\text{416}\text{K}$ ，熔化潜热 $L=\text{8}0\text{kJ}/\text{kg}$ 。

3.4. 网格划分

利用FLUENT前处理软件GAMBIT划分网格。正式计算之前，考虑网格尺寸和时间步长的无关性验证。为保证收敛效果，采用时间步长为1 s。然后，针对不同数量的网格9061、13,034及20,433，取某一点在某时刻的温度为参考变量进行比较。选取点T (0.05 m，0.05 m)，在加热功率q = 3000 W/m²的情况下加热到10⁴ s时，记录该点的温度，三种网格下分别为399.444 K、399.635 K、399.826 K，三种网格数下温度相对误差为0.048%，相差非常小。所以最终网格数量确定取13,034。典型的网格划分见图2。

4. 计算结果分析

4.1. 相界面变化

由于熔盐的固液相密度不同，所以导致相变过程中液相出现对流现象。在熔化开始阶段，固态熔盐最内层与加热棒外壁面直接接触，受到加热之后产生一层较薄的均匀液态熔融盐。随着加热之进行，液相熔化面积扩大，开始出现熔盐的自然对流。自然对流逐渐影响相界面，进一步加快熔盐熔化进程。

图3所示为加热棒热流密度q = 3000 W/m²的条件下，不同时刻混合熔盐熔化过程的相界面变化情况。图3(a)因导热作用固相熔盐升温，2.8小时之后固相与加热壁面之间开始形成很薄的液相区域，相界面为规则圆形。图3(b)液相区域渐渐扩大，相界面开始变得不规则，向上略微凸起，开始出现自然对流。图3(c)、图3(d)、图3(e)随着液相区域继续扩大，相界面不断移动，熔化后的熔盐内部温度由于分布不均匀而造成密度的变化。在浮升力的作用下密度较小的液态熔盐向上移动，换热过程以自然对流为主，热传导为辅。这样一来促进方腔上部分熔盐先于下部分熔化。图3(f)相界面消失，表示固相熔盐已全部熔化成液相。

Figure 2. Schematic of the computational grid

图2. 方腔网格划分示意图

(a) 2.8 h (b) 8.4 h (c) 11.2 h (d) 14.0 h (e) 19.6 h (f) 22.4 h

Figure 3. The interface change over time when q = 3000 W/m², Red presents liquid state, blue presents solid state

图3. 固液相变界面随时间变化图(q = 3000 W/m²时) (红色代表熔融状态，蓝色代表固态)

4.2. 不同时刻温度分布

图4给出不同时刻方腔内熔盐熔化过程的温度分布情况。图4(a)固态熔盐受热后温度开始上升，等温线之间呈一定梯度均匀递减，表明此时方腔的传热以导热为主。图4(b)靠近加热壁面的液态区域等温线向上略有波浪状，而固态部分的等温线仍是均匀分布，表明开始存在自然对流，但是仍以导热为主。图4(c)反映在浮升力作用下，液态熔盐受热，逐渐产生由热源中部向上运动，然后分别沿着左右壁面向下的流动。液相区也随着流动方向先至上而后分别沿左右向下不断扩大，自然对流使得等温区域增大而且等温线变得弯曲，在固相区等温线之间还存在较大的温度梯度。图4(d)反映，液相区占据了方腔的大部分区域，等温区域也随之扩大。根据方腔内熔盐温度分布图可知，由于固液相的熔盐存在密度差，受重力场下产生的浮升力的影响，熔盐出现热分层现象，方腔上部分区域熔盐温度比下部分熔盐温度高。

(a) 2.8 h (b) 8.4 h (c) 11.2 h (d) 19.6 h

Figure 4. Temperature distribution of the molten salt at different times

图4. 不同时刻的熔盐温度分布图

为了更清楚分析熔盐在熔化蓄热过程中温度变化规律，计算时分别设置了监测点A(0 m, 0.1 m)、B(0 m, −0.1 m)、C(0.1 m, 0 m)、D(0.1 m, 0.1 m)，其位置见图5中小插图。图5分别给出了A、B、C、D四点在整个熔化过程中温度的变化情况，根据曲线的整体趋势，大致可以将熔盐熔化蓄热过程分为三个阶段：固相熔盐显热蓄热阶段，随着加热时间增加，熔盐从初始温度298 K升温至熔化温度415 K，熔盐温度上升曲线斜率较大。接着进入熔盐相变蓄热阶段，随着熔化过程的进行，熔盐开始发生相变，吸收大量的潜热，此阶段熔盐会在415 K左右小幅度浮动，温度增加减缓，温升曲线的斜率变小。相变吸热持续一段时间后，进入液相熔盐显热蓄热阶段，熔盐液相区扩大，液相的熔盐继续吸热升温直至加热结束。

分析图5可知，因为熔化的液相熔盐密度小于固相熔盐，受浮升力影响液相熔盐开始向上流动，引起自然对流。点A在方腔加热棒上部，很快受自然对流影响，所以升温快，完全熔化后持续受热保持较高温度。点D离加热中心比A、B、C远，所以在开始阶段因导热传热升温较慢，温度比较低，但是因点D处于方腔上部，在11小时以后液相扩大到此处，受自然对流影响，温升速率加快，先于方腔中下部分的熔盐完成熔化。

4.3. 不同时刻的流场

图6为不同时刻混合熔盐熔化过程流场矢量图。由图可以看出，图6(a)靠近加热壁面附近的熔盐刚刚开始形成液相，因此速度不明显。图6(b)熔盐液相区开始扩大，流体经热源中部向上，然后分别沿着左右壁面向下。图6(c)在液相区域开始出现两个明显且对称的涡流，最大流速出现在涡流的交界处。图6(d)随着液相区继续扩大，涡流动范围也扩大，液相中涡流速也不断增大，方腔上部分的液态流速大于下部分。

4.4. 熔盐液相比变化

液相比反映熔盐材料的液相体积份额在两相状态时总体积份额中的占比。图7、图8分别是变加热

Figure 5. Temperature variation of the monitored points during the melting process at q = 3000 W/m²

图5. 方腔内部分监测点的温度变化(q = 3000 W/m²时)

功率和不同加热棒直径时熔盐液相比变化情况。

图7所示为当加热棒直径d =10 cm时，不同的加热功率下熔盐的液相比–时间曲线图。从曲线整体趋势看，加热功率越大，熔盐完全熔化需要的时间越短。当加热功率为500 W/m²时，在监测前28小时内熔盐依然处于固态，不出现熔化现象。当功率为1000 W/m²时，从13.5 h时刻左右熔盐开始出现熔化现象。当功率逐渐增大，熔盐开始熔化的时间明显提前。当功率为2000 W/m²时，熔盐很快就有了明显的熔化现象，在28 h时刻几乎完全熔化。当功率为3000 W/m²时，混合熔盐完全熔化的时间约为22.2小时，熔化的速率是2000 W/m²时速率的1.25倍，熔化时间明显缩短。

图8所示为当加热功率q = 3000 W/m²时，不同加热棒直径下的熔盐熔化情况。从曲线整体趋势看，加热棒直径越大，熔盐完全熔化需要的时间越短。当加热棒直径为5 cm时，从曲线明显地看出熔盐很长

(a) 2.8 h (b) 8.4 h (c) 11.2 h (d) 19.6 h

Figure 6. The flow patterns of the molten salt at different heating time

图6. 熔盐熔化过程流场矢量图

Figure 7. The profile of the melting fraction for different heating power at d = 10 cm

图7. 不同加热功率，液相比随时间变化曲线(d = 10 cm时)

Figure 8. The change of melting fraction for different heating bar size at q = 3000 W/m²

图8. 不同加热棒直径，液相比随时间变化曲线(q = 3000 W/m²时)

一段时间基本上处于固态，没有明显熔化的现象。随着加热进行，从14小时以后，熔化开始加速，液相比明显增大。与其他两组比较，d = 5 cm情况下熔盐从固相转化到液相的时间较长。当加热棒直径为10 cm时，正是前述变功率的情形，熔盐全部熔化的时间约为22.2小时。当加热棒直径为15 cm时，熔盐完全熔化的时间约为14.0小时，熔化的速率是直径10 cm时的1.6倍，熔盐完全熔化所需时间明显缩短。但是，在热流密度q不变时增大加热棒直径，实质是总加热功率增大了。

5. 结论

通过对内置圆形加热棒，四周绝热的方腔内混合熔盐熔化动态传热过程的非稳态数值模拟，得到如下结论。1) 熔盐熔化过程中，液相比受加热功率和加热棒尺寸影响。在文章研究的参数范围内，当计算工况的加热棒直径为15 cm，加热功率为3000 W/m²，熔盐熔化至少需要14小时。2) 熔盐从固相加热熔化至液相熔融状态，经过热传导阶段，相变阶段，和自然对流主导阶段。3) 自然对流使得液相等温区域逐渐扩大，等温线变得弯曲，且液相温度梯度小于固相。4) 受液相熔盐涡流影响，方腔上部分的熔盐先于下部分熔盐完成熔化。涡流的最大速度出现在对称涡流交界面处。5) 方腔内熔化后的液相熔融盐出现热分层现象。

基金项目

国家重点基础研究发展规划(973)基金资助项目(编号：NO.2015CB251303)。