基于格子玻尔兹曼模型的微米级铝颗粒氧化层传热熔化特性

Heat transfer and melting characteristics of micronmeter-sized aluminum particle oxide layers based on lattice Boltzmnn model

摘要

建立了基于焓的二维格子玻尔兹曼模型（LBM），研究了不同边界条件以及不同颗粒粒径下微米级铝颗粒氧化层在熔化过程中的温度分布、平均液相分数和熔化结束时间。结果表明：自然对流浮升力对熔化特性有着显著影响；在单个加热壁加热时，自然对流的浮升力使下壁面加热时的熔化进程最快，上壁面时的最慢；当两个加热壁共同作用时，上下加热壁条件下的熔化最快，与左右加热壁的条件相比，熔化时间缩短了18.05%；随着加热壁的数量增多，熔化效率的提升效果减弱。通过计算发现，与单加热壁相比，双加热壁、三加热壁和四加热壁的熔化效率分别提升了21%、73%和75%。冷壁面的存在使熔化进程变慢，但影响仅存在于熔化进程的后半段。此外，颗粒粒径的增大可以加速熔化，但这一特征也会随着粒径的增大逐渐减弱。

Abstract

In this paper， a two-dimensional lattice Boltzmann model based on enthalpy was established to study the temperature distribution， average liquid fraction and melting end time of the oxide layer of micronmeter-sized aluminum particles under different boundary conditions and particle sizes. The results show that the natural convection buoyancy has a significant effect on the melting characteristics. When a single heating wall is heated， due to the influence of natural convection buoyancy， the melting process of the lower wall is the fastest， and the melting speed of the upper wall is the slowest. When the two heating walls act together， the melting speed of the upper and lower heating walls is the fastest， and the melting time is shortened by 18.05% compared to the conditions of the left and right heating walls. As the number of heating walls increases， the improvement effect of melting efficiency weakens. Through calculation， it was found that compared with a single heating wall， the melting efficiency of double heating walls， three heating walls， and four heating walls increased by 21%， 73%， and 75%， respectively. The presence of a cold wall slows down the melting process， but the impact only exists in the latter half of the melting process. In addition， an increase in particle size can accelerate melting， but this characteristic will gradually weaken as the particle size increases.

Graphical abstract

关键词

工程热物理 / 格子玻尔兹曼 / 液相分数 / 熔化特性

Key words

engineering thermophysics / lattice Boltzmann / liquid fraction / melting characteristics

0 引　言

金属燃烧剂是固体复合推进剂的基本组分之一^［1］。金属颗粒经过燃烧后能够释放出大量的热，将其加入固体推进剂后可以显著增加推进剂的密度和暴热。在金属颗粒中，铝颗粒由于其具有较低的成本，获取方式简单，在燃烧时会显著提高燃烧热，被普遍应用于提高推进剂的燃烧热和密度^［2-4］。

铝颗粒燃烧过程十分复杂，包括复杂的相变和传热传质过程，其中伴随着常温氧化、点火、熔化、燃烧以及燃烧产物的凝聚^［5，6］。由于铝的活性比较高，微米级铝粉颗粒在室温下会即刻发生氧化反应。氧化后的铝颗粒外围包裹着一层无定形氧化铝，这些氧化铝会阻碍铝粉的进一步氧化和燃烧点火时间，从而导致点火延迟。因此，有关氧化层的实验和模拟研究对后续铝氧燃烧及点火温度至关重要^［7，8］。在实验方面，Merzhanov等^［9］研究了气体氧化剂对细铝线电加热点火的实验结果，假设点火温度和A1₂O₃的熔点相一致，为23 00 K，并预计其活化能为17 kcal/mol。Breiter等^［10］基于Glassman的工作得出金属可以在扩散状态下点火，在低放热率条件下，着火过程主要发生在动力学控制区域，与金属及其氧化物的物理性质关系不大；而在高放热率条件下，着火则主要由金属氧化物在金属中的溶解度或相变的开始所决定。Ermakov等^［11］得出加热过程中会产生热机械应力，这使在加热过程中，温度还未达到氧化层的熔化温度时氧化层的完整性已经破裂。在数值模拟方面，Desjardin 等^［12］建立了单个铝颗粒在空气中燃烧的机理模型，提出了一种新的守恒标量公式，该公式考虑了金属氧化物在熔融铝表面的沉积，并且该模型与实验结果有良好的一致性。Yang等^［13］提出了一个悬浮在空气中的铝颗粒燃烧的简化模型，考虑对流和辐射加热，重点研究了氧化铝对铝颗粒燃烧的影响，结果表明：氧化膜内部的铝由于非均相反应以及来自环境空气的热传递而熔化，直到达到氧化层的熔点开始燃烧。Volkov^［14］分析了固体火箭发动机燃烧室内湍流气体颗粒流动，建立了含有氧化层的单个铝液滴两相流动的数学模型，并与实验有良好的一致性。从目前的研究来看，国内外学者针对铝粉燃烧开展了大量的理论和实验研究，但主要集中在铝粉颗粒上，对颗粒表面的氧化物研究较少。

近年来，格子玻尔兹曼方法已经成功地被应用于介观尺度下流体的流动、传热及化学反应，成为传统数值方法的有力替代^［15］。Jourabian等^［16］采用单相流模型，使用基于焓的LBM模型，研究了相变材料在倾斜椭圆环中的约束熔化，忽略了相变材料的过冷问题，确定了增强传热的方法，发现倾斜不会导致任何熔化过程中液体分数的变化。Ghasemi等^［17］开发了基于焓的格子玻尔兹曼程序，对以嵌入相变材料颗粒为代表的圆形多孔区域内的对流熔化进行研究，得到了孔隙率、渗透率、有效导热系数和熔化潜热的影响规律。Chen等^［18］将轴对称LBM模型与焓更新模型相结合，采用双分布函数和多重弛豫时间建立了多孔介质中瞬态相变固液相变LBM模型，用来模拟多孔介质中的流场和温度场，并验证了该模型的准确性。He等^［19］建立了基于焓方法的二维格子玻尔兹曼模型，研究了气体扩散层中冰的熔化过程，得到了气体扩散层熔化位置、气孔率、瑞利数和斯特凡数对熔化速率的影响。可以看出，格子玻尔兹曼方法已经成功用于相变材料的熔化、传热、传质等，与宏观方法相比，LBM有良好的并行计算能力和便于进行复杂边界计算等优点，可以很好地用于微米级铝粉氧化层熔化机理的研究，但研究者发现，由于重力和自然对流引起浮升力的作用，导致不同壁面对熔化的作用机理不同，且目前对加热壁面数量及不同位置对熔化的影响研究还不是很全面，导致这方面的物理过程，尤其是冷壁的作用机理还没有得到具体的解释。

在本文研究中，使用Matlab软件建立了二维流动LBM模型和基于焓的热LBM模型，通过更新总焓来跟踪相界面，研究了不同颗粒直径和不同加热壁数目与位置对微米级铝粉颗粒氧化层熔化行为的影响。

1 数学模型

1.1　物理模型

本文主要研究微米级铝颗粉粒氧化层的熔化过程。在计算中，图1（a）为铝颗粒的SEM图片，图1（b）为采用计算机随机算法构建氧化铝熔化的二维物理模型，铝颗粒外围包裹一层固体氧化铝。固体氧化物作为氧化反应的产物，厚度是不均匀的，为了简化模型，方便进行计算，在建模过程中将其视为均匀厚度。本文设定铝颗粒半径为32 μm，在格子模型构建中采用正态分布，颗粒在平面内随机分布。固体氧化铝包裹在铝颗粒外围，厚度为4 μm。

在研究不同加热壁面位置及数量时，本文建立了图2所示的9种模型。其中模型1~3为不同位置的单壁加热，主要比较重力和热浮升力带来的差异，模型4~6为双壁面加热，主要分为相邻双加热壁、相对双加热壁，模型7和8分别给定了3个加热壁面和4个加热壁面，用来比较不同数量的加热壁面对熔化速度的促进效果，模型9设定右壁面为冷壁，比较模型1和9可以研究冷壁在熔化过程中的作用机理，该模型充分考虑了不同加热壁数目和位置对熔化的影响。在计算中，本文建立了双分布函数，一个用于速度场，一个用于温度场，采用D₂Q₉格子的BGK近似方法模拟了固液相界面的温度场和速度场，并且基于焓法跟踪固液相界面。

1.2　格子玻尔兹曼模型

在计算中建立了双分布函数，一个用于速度场，一个用于温度场，接下来将介绍本研究中的LBM模型。

考虑浮力影响的密度分布函数f_i （x，t）由式（1）给出^［15］：

f_i^{\mathrm{e}\mathrm{q}}\left(\mathit{x}+{\mathit{e}}_{\mathit{i}}\mathrm{\Delta }t,t+\mathrm{\Delta }t\right)=f_i\left(\mathit{x},t\right)-

\left.\frac{\mathrm{1}}{{\tau }_{\mathrm{f}}}\left[f_i\left(\mathit{x},t\right)-\right.f_i^{\mathrm{e}\mathrm{q}}\left(\mathit{x},t\right)\right]+\mathrm{\Delta }t{F}_i（1）

式中： x 为离散点的位置； e_i 为离散速度， x + e_i ∆t为它在某一时刻沿着某一方向迁移的距离；τ_f为运动黏度相关的弛豫时间；Δt为时间增量；f_i^eq（x，t）为密度平衡分布函数。

在本研究中使用D₂Q₉模型，离散速度 e_i 可写为式（2）^［15］：

{\mathit{e}}_i=\left\{\begin{array}{l}\left(\mathrm{0,0}\right),i=\mathrm{0}\\ \left(\mathrm{c}\mathrm{o}\mathrm{s}\left[\left(i-\mathrm{1}\right)\frac{\mathrm{\pi }}{\mathrm{2}}\right],\mathrm{s}\mathrm{i}\mathrm{n}\left[\left(i-\mathrm{1}\right)\frac{\mathrm{\pi }}{\mathrm{2}}\right]\right),\\ i=\mathrm{1,2},\mathrm{3,4}\\ \sqrt[]{\mathrm{2}}{\left(\mathrm{c}\mathrm{o}\mathrm{s}\left[\left(\mathrm{2}i-\mathrm{9}\right)\frac{\mathrm{\pi }}{\mathrm{4}}\right],\mathrm{s}\mathrm{i}\mathrm{n}\left[\left(\mathrm{2}i-\mathrm{9}\right)\frac{\mathrm{\pi }}{\mathrm{4}}\right]\right)}_{},\\ i=\mathrm{5,6},\mathrm{7,8}\end{array}\right.（2）

密度平衡分布函数f_i^eq（x，t）由式（3）给出^［15］：

f_i^{\mathrm{e}\mathrm{q}}={\omega }_i\rho \left[\mathrm{1}+\frac{{\mathit{e}}_{\mathit{i}}\cdot \mathit{u}}{c_{\mathrm{s}}^{\mathrm{2}}}+\frac{{\left({\mathit{e}}_i\cdot \mathit{u}\right)}^{\mathrm{2}}}{\mathrm{2}{c}_{\mathrm{s}}^{\mathrm{4}}}-\frac{{\left|\mathit{u}\right|}^{\mathrm{2}}}{\mathrm{2}{c}_{\mathrm{s}}^{\mathrm{2}}}\right]（3）

式中：c_s为晶格单位表示的声速；c_s=c/\sqrt[]{\mathrm{3}}，c为晶格单位表示的速度。c=∆x/∆t，其中，对于D₂Q₉模型，式（4）给出了权重系数ω_i^［19］：

{\omega }_i=\left\{\begin{array}{l}\mathrm{4}/\mathrm{9},i=\mathrm{0}\\ \mathrm{1}/\mathrm{9},i=\mathrm{1,2},\mathrm{3,4}\\ \mathrm{1}/\mathrm{36},i=\mathrm{5,6},\mathrm{7,8}\end{array}\right.（4）

体力项F_i 可以定义为式（5）^［16］：

F_i={\omega }_i\left(\mathrm{1}-\frac{\mathrm{1}}{\mathrm{2}{\tau }_{\mathrm{f}}}\right)\left(\frac{{\mathit{e}}_{\mathit{i}}-\mathit{u}}{c_{\mathrm{s}}^{\mathrm{2}}}+\frac{{\mathit{e}}_{\mathit{i}}\cdot \mathit{u}}{c_{\mathrm{s}}^{\mathrm{2}}}\right)·\left(\rho f\right)（5）

式中：f为流体流动的热浮升力，它可以由式（6）计算^［16］。

f=\rho g\beta \left(T-T_{\mathrm{r}\mathrm{e}\mathrm{f}}\right)（6）

速度场的松弛时间τ_f可以由式（7）计算^［15］：

{\tau }_{\mathrm{f}}=\frac{\mathrm{3}\upsilon }{c^{\mathrm{2}}\mathrm{\Delta }t}+\mathrm{0.5}（7）

则宏观量可以由式（8）（9）计算出来^［15］：

\rho =\sum _i{f}_i（8）

\rho u=\sum _i{f}_i{\mathit{e}}_{\mathit{i}}+\frac{\mathrm{\Delta }t}{\mathrm{2}}f（9）

对于流体和固体区域的换热，温度分布函数g_i （x，t）可由式（10）计算^［19］：

\begin{array}{c}{g}_i^{\mathrm{e}\mathrm{q}}\left(\mathit{x}+{\mathbf{e}}_{\mathit{i}}\mathrm{\Delta }t,t+\mathrm{\Delta }t\right)=g_i\left(\mathit{x},t\right)-\\ \frac{\mathrm{1}}{{\tau }_g}\left[g_i\left(\mathit{x},t\right)-g_i^{\mathrm{e}\mathrm{q}}\left(\mathit{x},t\right)\right]\end{array}（10）

式（11）给出了密度平衡分布函数g_i^eq（x，t）^［19］：

g_i^{\mathrm{e}\mathrm{q}}={\omega }_{i}T\left[\mathrm{1}+\frac{{\mathit{e}}_{\mathit{i}}\cdot \mathit{u}}{c_{\mathrm{s}}^{\mathrm{2}}}+\frac{{\left({\mathit{e}}_{\mathit{i}}\cdot \mathit{u}\right)}^{\mathrm{2}}}{\mathrm{2}{c}_{\mathrm{s}}^{\mathrm{4}}}-\frac{{\left|\mathit{u}\right|}^{\mathrm{2}}}{\mathrm{2}{c}_{\mathrm{s}}^{\mathrm{2}}}\right]（11）

温度可以被式（12）计算出^［17］：

T=\sum _i{g}_i（12）

为了处理固液相界面，采用基于焓的方法，局部热焓计算公式为式（13）^［16］：

E{n}^k=C_{\mathrm{p}}{T}^k+L_{\mathrm{f}}{f}_{\mathrm{l}}^{k-\mathrm{1}}（13）

式中：En^k 和T^k 分别为当前时间步长的热焓和温度；f_l^k-1为前一个步长的液体分数。

在固液相界面计算过程中，温度分布函数可以被修正为式（14）^［19］：

\begin{array}{c}{g}_i\left(\mathit{x}+{\mathit{e}}_{\mathit{i}}\mathrm{\Delta }t,t+\mathrm{\Delta }t\right)-g_i\left(\mathit{x},t\right)=\\ -\frac{\mathrm{1}}{{\tau }_g}\left[g_i\left(\mathit{x},t\right)-g_i^{\mathrm{e}\mathrm{q}}\left(\mathit{x},t\right)\right]-\end{array}

{\omega }_i\frac{L_{\mathrm{f}}}{C_{\mathrm{p}}}\left[f_{\mathrm{l}}^k\left(\mathit{x}\right)-f_{\mathrm{l}}^{k-\mathrm{1}}\left(\mathit{x}\right)\right]（14）

其中，局部热焓可以由式（15）确定^［18］：

（15）

固液相界面更新是通过局部焓变化来追踪的。在处理对流熔化时，f_l 取值为0~1，计算仅对固液相进行捕捉，并且相间完全分离，因此，不考虑固液相变的糊状区。使用0.5作为判断标准可近似实现相变和相分离的过程。当液体分数f_l≥0.5时，在计算中认定为液体，发生流体流动模型中的碰撞步骤。另外，当液体分数f_l＜0.5时，在计算中认定为固体，对流体流动采用反弹边界^［19］。

2 模型验证

为了验证基于焓的LBM方法处理问题的适用性和准确性，首先需要进行模型验证，将方腔内纯物质的导热熔化及自然对流熔化解析解和模拟解进行比较^{［19，20］}。在计算过程中，空腔内由纯固体填充，纯固体的初始温度为T₀，然后在其左壁施加恒温T_b，随着左壁的不断加热，纯物质开始熔化。

验证模型分为导热熔化过程和对流熔化过程。在纯物质导热熔化模型中，左侧壁面设置高温，其余3个壁面均采用绝热反弹边界，如图3（a）所示。

对于纯物质的对流熔化，根据实验条件，左壁面设置高温，右壁面为熔化温度，上下壁面为绝热壁面。壁面设置为无滑移边界，采用反弹格式，相同的设置被应用于固体颗粒。

对于多孔结构的对流熔化计算，边界条件及模型如图3（b）所示，处理方法与验证模型中纯物质的对流熔化相同。对于壁面和铝颗粒均采用回弹处理。对于对流和传导融化过程使用了计算域为500×500的网格尺寸，网格分辨率为1 μm。这两个网格已经被证明是网格无关的。

首先，对纯物质的导热熔化进行模拟计算，这是一个一维问题，存在一个解析解。在导热熔化过程中，主要观察固液相界面，将LBM计算值与解析解固液相界面随时间变化情况进行比较。

图4（a）给出了导热熔化过程中LBM计算值与解析解的无量纲温度比较。其中，T*为温度与加热温度的比值；X*为监测点的位置与总长度的比值。从计算可知，解析解与模拟值具有良好的一致性。

导热的热传递过程可用式（16）描述^［19］：

\frac{\partial T}{\partial t}=\alpha \left(\frac{{\partial }^{\mathrm{2}}T}{\partial x^{\mathrm{2}}}+\frac{{\partial }^{\mathrm{2}}T}{\partial y^{\mathrm{2}}}\right)（16）

式中：T、t和α分别为温度、时间和热扩散率。对于二维计算域，方程（17）可被用于固液相界面^［19］：

\alpha {\left(\frac{\partial T}{\partial x}\right)}_{x=x_{\mathrm{m}}}=\frac{L_{\mathrm{f}}}{C_{\mathrm{p}}}\frac{\mathrm{d}{x}_{\mathrm{m}}}{\mathrm{d}t}（17）

式中：L_f、C_p、α和x_m分别为熔融潜热、比热、热扩散系数和固液相界面位置。

当体系的初始温度为熔融温度时，方程（16）在界面满足式（19）时的解析解可由式（18）计算^［17］：

（18）

x_{\mathrm{m}}\left(t\right)=\mathrm{2}\lambda \sqrt[]{\alpha t}（19）

式中：λ由式（20）确定^［19］。

\frac{c_{\mathrm{p}\mathrm{l}}\left(T_{\mathrm{h}}-T_{\mathrm{m}}\right)}{L_{\mathrm{f}}\mathrm{e}\mathrm{x}\mathrm{p}({\lambda }^{\mathrm{2}})\mathrm{e}\mathrm{r}\mathrm{f}\left(\lambda \right)}-\frac{c_{\mathrm{p}\mathrm{s}}(T_{\mathrm{m}}-T_i)}{L_{\mathrm{f}}\mathrm{e}\mathrm{x}\mathrm{p}({\lambda }^{\mathrm{2}})\mathrm{e}\mathrm{r}\mathrm{f}c\left(\lambda \right)}=\lambda \sqrt[]{\mathrm{\pi }}（20）

在LBM模拟中，氧化铝的物理性质为：C_p=1 250 J/（kg∙K），ρ=3.9g/cm³，T_m=2 054 ℃，L_f=16.435×10⁶。在计算中，将真实物理单位转化为LBM单位，转化结果如表1所示。图4（a）显示了LBM计算结果与分析结果之间的比较。可以看出，基于焓的LBM模型可以用于后续的计算过程。

对于对流问题，平均Nusselt number由式（21）计算^［21］：

\overline{Nu}=\frac{\overline{h}H}{k}（21）

式中：\overline{Nu}为平均换热系数；k为流体的导热系数。

对流熔化不同于传导熔化，并且没有线性特征，因此，使用平均努塞尔数进行验证，并计算相关方程（22）的结果进行比较^{［22，23］}：

Nu=\frac{{\left(\mathrm{\pi }Fo\right)}^{-\mathrm{0.5}}}{\mathrm{e}\mathrm{r}\mathrm{f}\left(\lambda \right)}+\left[\mathrm{0.29}R{a}^{\mathrm{0.25}}-\frac{{\left(\mathrm{\pi }Fo\right)}^{-\mathrm{0.5}}}{\mathrm{e}\mathrm{r}\mathrm{f}\left(\lambda \right)}\right]\cdot

{\left\{\mathrm{1}+{\left[\mathrm{0.0175}R{a}^{\mathrm{0.75}}{\left(\frac{\alpha t}{l^{\mathrm{2}}Ste}\right)}^{\mathrm{1.5}}\right]}^{-\mathrm{2}}\right\}}^{-\mathrm{0.5}}（22）

式中：Fo表示傅里叶数；瑞利数Ra是与浮力驱动流相关的无量纲数，用来区分自然对流与强制对流；斯特凡数Ste代表熔化过程中显热与潜热之比。

从图4（b）可以看出，基于焓的LBM模型的计算值与用于对流熔化的解析解非常吻合。结果表明：LBM模型具有良好的计算精度，可用于后续计算。

3 结果与讨论

不同加热壁和不同颗粒尺寸对金属熔化有着重要影响，在计算中，本文在不同条件下分析了温度分布、熔化进程以及熔化结束时间。当熔化刚开始时，计算域内固相分数为1，液相分数为0，随着时间增加，液相分数开始增大。同时，传热性能的差异直接影响熔化结束时间，因此，本文通过分析温度分布和熔化结束时间判断不同条件下的传热性能。

3.1　不同加热壁面对熔化的影响

不同加热条件对氧化层熔化有着重要影响，在这一部分，将讨论不同加热墙的位置及数目对熔化的影响。

3.1.1　单层加热壁对熔化的影响

本节讨论了处在不同位置下的单加热壁（见图5）对熔化进程的影响。在计算过程中，本文引入了瑞利数（Ra）、斯特凡数（Ste）和普朗特数（Pr）来描述熔化行为。普朗特数表示流体中能量和动量迁移过程相互影响的无因次组合数，表明温度边界层和流动边界层的关系。无量纲数设置为：Ra=10⁵，Ste=1，Pr=0.73，孔隙率为20%。同时，傅里叶数（Fo）是一个用来描述非稳态热传导及分子扩散的无因次量，本文引入它来表示无因次化的时间。熔融过程中部分固体点发生相变转化为液相，完全可以使用某个时刻下发生相变的格点数量来衡量熔融的效率。在后续的计算中，本文引入了熔化速度来衡量熔化进程，即发生相变的固体点数量与总固体点数量之比，熔化速度是无量纲化的。图5给出了熔化过程中的温度分布图和熔化进程图。

从图5可以看出，在熔化初期，整个热传递过程是导热主导的，可以看到清晰的相界面，随着时间的延长，对流在熔化过程中发挥主要作用。模型（c）的熔化速度最快，这是由于浮升力的存在促进了熔化进程，同时使得模型（b）的温度扩散减缓。

图6（a）给出了不同单层加热壁的熔化结束时间，根据计算，下加热壁的熔化结束时间比上加热壁快90%，比左加热壁快64.18%，这说明在熔化过程中，对流引起的浮升力对熔化的影响十分巨大。

3.1.2　不同双层加热壁对熔化的影响

这一部分研究了不同位置双层加热壁（模型（d）（e）（f））对温度分布以及平均液体分数的影响，相应的边界条件如图2所示。图7显示了不同双层加热壁下的温度分布和熔化进程。从图中可以看出，在Fo=0.01时，模型（d）（e）（f）的温度变化和熔化进程变化都不大，这是由于在熔化初期，热传导是主要的熔化机制，各种外力的作用并不明显，当Fo＞0.03时，相界面均开始弯曲，且腔内下部物质熔化快于上部，这一现象在模型（e）中尤为明显，这主要是由于模型（e）下壁面受到自然对流热浮升力的影响，使温度变化和熔化进程加快。

随着Fo数的增加，熔化进程受对流的影响越来越明显，方腔内上下部的熔化差异越来越大。图6（b）为不同双壁下的熔化结束时间，从图中可以看到，模型（d）熔化速度最慢，熔化结束时间最长，这是由于模型（d）的上边界为加热边界，浮升力的存在使上边界传热减缓。通过对比模型（e）和（f）发现，熔化前期，两模型熔化速度差别不大，在Fo＞0.04后，熔化速度开始有明显差别，这说明浮升力作用开始明显起来。通过计算，模型（e）的熔化速度比模型（d）提升了61.43%，模型（f）的熔化速度比模型（d）提升了52.94%，模型（e）的熔化速度比模型（f）提升了18.05%。

3.1.3　不同加热壁数量对熔化进程的影响

本节分别讨论了单层加热壁（模型（a））、双层加热壁（模型（d））、三层加热壁（模型（g））和四层加热壁（模型（h））对熔化进程的影响，对于上述模型，无量纲数分别设置为Ra=10⁵，Ste=1。图8给出了不同加热壁面的熔化过程图和温度变化图。从图中可以明显看出，加热壁面数量的增加对熔化过程有着明显的促进作用，随着加热壁的数量增多，熔化效率的提升效果减弱。尤其当从双加热壁增加至三加热壁时，熔化结束时间缩短了65.7%。在熔化初期，导热占据整个熔化过程的主导地位，熔化进程图和温度云图存在清晰的相界面。随着熔化时间的延长，相界面开始变得弯曲，自然对流开始主导整个熔化过程。

图6（c）给出了不同加热壁数量的熔化结束时间。从图中可以看到，加热壁面的数量对熔化过程有着显著影响，通过计算，模型（d）比模型（a）的熔化时间缩短了21%，模型（g）比模型（a）的熔化时间缩短了73%，模型（h）比模型（a）的熔化时间缩短了75%。

3.1.4　冷壁对熔化进程的影响

冷壁在熔化过程中的影响是很重要的。在之前的研究中，很少有人探究冷壁存在对熔化的影响，本节将讨论冷壁对熔化的影响。模型（a）和（i）的差异如图1所示，计算无量纲参数设置为Ra=10⁵，Ste=1。图9给出了模型（a）和（i）在Fo为0.01~0.2的温度变化和熔化进程图。从图中可以看出，在Fo＜0.15时，温度变化和平均液相分数几乎没有差异，随着时间的延长，冷壁的存在使自然对流的效果更加明显。图6（d）给出了单层绝热壁和冷壁的熔化结束时间，从图中可以更加明显地看出，冷壁的存在对熔化的影响仅存在于熔化进程的后半段，冷壁使熔化进程变慢。

3.2　不同颗粒尺寸对熔化进程的影响

不同颗粒尺寸对熔化进程有着重要的影响。在本节计算中，Ra=10⁵，Ste=1，孔隙率为20%，计算了颗粒粒径分别为4、6、8、10 μm时的温度分布变化和熔化进程。图10给出了不同颗粒粒径下的温度分布和熔化进程，从图10中可以看到，当颗粒粒径变化时，温度分布图变化不明显。当Fo数较小时，导热主导着整个熔化进程，温度分布图和熔化进程图均存在明显的相界面。当Fo>0.03时，固液相界面开始弯曲，此时，对流在熔化中起主导作用。

图11给出了不同颗粒粒径下的熔化结束时间，从图中可以看出，随着粒径的增大，熔化速度随之加快，且在小粒径中尤为明显。这是由于，随着颗粒粒径的增大，颗粒与颗粒之间的间隙增大，这使自然对流的作用加强，同时，随着颗粒粒径的增大，颗粒与空气的接触面积增大，这也加快了熔化进程。

4 结　论

（1）自然对流浮升力的影响在熔化过程中起着重要作用。在单个加热壁，由于自然对流浮升力的影响，下壁面熔化速度最快，上壁面熔化速度最慢。当两个加热壁存在时，上下两加热壁熔化速度最快。

（2）对于不同数量的加热壁，双加热壁比单加热壁熔化进程加快了21%，三加热壁比双加热壁熔化进程加快了65.7%，四加热壁比三加热壁熔化进程加快了7.7%，加热壁数量的增加使熔化进程加快，且随着加热壁的数量增多，熔化效率的提升效果减弱。

（3）当将右侧绝热壁换成冷壁时，冷壁的存在使熔化进程变慢，但影响仅存在于熔化进程的后半段。

（4）当颗粒粒径增大时，熔化速率加快。这一规律在小粒径时更为明显。

（5）根据研究所得规律，综合考虑熔化效率和加热成本，建议使用3个加热壁，其中由于浮升力的强化换热作用需要包含底部加热面。此外，选择适当的颗粒粒径（8 μm）也可促进固体复合推进剂的熔化。

参考文献

[1] 田入园. 复合固体推进剂中Al颗粒燃烧行为数值模拟[D]. 南京: 南京理工大学能源与动力工程学院, 2017.

[2] Tian Ru-yuan. Numerical simulation of combustion behavior of aluminum particle in composite solid propellants[D]. Nanjing: School of Energy and Power Engineering, Nanjing University of Science and Technology, 2017.

[3] 胡松启, 李葆萱. 固体火箭发动机燃烧基础[M]. 西安: 西北工业大学出版社, 2015.

[4] 王克秀. 固体火箭推进剂及燃烧[M]. 北京:国防工业出版社, 1983.

[5] 萨顿·比布拉兹. 火箭发动机基础[M]. 洪鑫, 张宝炯, 等译. 北京: 科学出版社, 2003.

[6] Fang C, Fen L S. Experimental studies on effects of AP content and particle size in NEPE propellant[J]. Journal of Solid Rocket Technology, 2001, 24(3): 47-53.

[7] 唐伟强, 杨荣杰, 李建民, 高铝固体推进剂中氟化物促进铝燃烧研究进展[J]. 固体火箭技术, 2020,43(6): 679-686.

[8] Tang Wei-qiang, Yang Rong-jie, Li Jian-min et al. Research progress of fluorides in high aluminum solid propellant to promote aluminum combustion[J]. Solid Rocket Technology, 2020, 43(6): 679-686.

[9] Jeurgens L P H, Sloof W G, Tichelaar F D, et al. Thermodynamic stability of amorphous oxide films on metals: application to aluminum oxide films on aluminum substrates[J]. Physical Review B, 2000, 62(7): 4707-4719.

[10] Jeurgens L P H, Sloof W G, Tichelaar F D, et al. Structure and morphology of aluminiumoxide films formed by thermal oxidation of aluminium[J]. Thin solid Films, 2002, 418(2): 89-101.

[11] Merzhanov A G, Grigorjev Y M, Galchenko Y A. Aluminum ignition[J]. Combustion and Flame, 1977, 29: 1-14.

[12] Breiter A L, Mal´Tsev V M, Popov E I. Models of metal ignition[J]. Combustion, Explosion and Shock Waves, 1977, 13(4): 475-484.

[13] Ermakov V A, Razdobreev A A, Skorik A I, et a1. Temperature of aluminum particles at the time of ignition and combustion[J]. Combustion, Explosion and Shock Waves, 1982, 18(2): 256-257.

[14] Desjardin P E, Felske J D, Carrara M D. Mechanistic model for aluminum particle ignition and combustion in air[J]. Journal of Propulsion and Power, 2005, 21(3): 478-485.

[15] Yang H, Lee J, Kim K, et a1. Simplified model for single aluminum particle combustion[C]∥The 47th AIAA Aerospace Science Meeting Including, The New Horizons Forum and Aerospace Exposition, Orlando,USA, 2009: 1-9.

[16] Volkov K N. Combustion of Single Aluminum Droplet in Two-phase Flow[M]. Hauppauge: Nova Science, 2011.

[17] Mohammadian S K, Zhang Y. Convection heat transfer with internal heat generation in porous media: implementation of thermal lattice Boltzmann method[J]. Numerical Heat Transfer Part A: Applications, 2019, 76(3): 101-114.

[18] Jourabian M, Darzi A, Akbari O A, et al. The enthalpy-based lattice Boltzmann method(LBM) for simulation of NePCM melting in inclined elliptical annulus[J]. Physica A: Statistical Mechanics and its Applications, 2019, 548: No.123887.

[19] Ghasemi K, Tasnim S, Mahmud S. Shape-stabilized phase change material convective melting by considering porous configuration effects[J]. Journal of Molecular Liquids, 2022, 355: No.118956.

[20] Chen D, Riaz A, Aute V C, et al. A solid–liquid model based on lattice Boltzmann method For phase change material melting with porous media in cylindrical heat exchangers[J]. Applied Thermal Engineering, 2022, 207: No.118080.

[21] He P, Chen L, Mu Y, et al. Lattice Boltzmann method simulation of ice melting process in the gas diffusion layer of fuel cell[J]. International Journal of Heat and Mass Transfer, 2020, 149(2): No.119121.

[22] Eshraghi M, Felicelli S D. An implicit lattice Boltzmann model for heat conduction with phase change[J]. International Journal of Heat and Mass Transfer, 2012, 55(9, 10): 2420-2428.

[23] Zhao C Y, Dai L N, Tang G H, et al. Numerical study of natural convection in porous media (metals) using lattice Boltzmann method(LBM)[J]. International Journal of Heat and Fluid Flow, 2010, 31(5): 925-934.

[24] Jiang Z H, Yang G G, Li S A, et al. Investigation of the ice melting process in a simplified gas diffusion layer of fuel cell by the lattice Boltzmann method[J]. Energy and Fuel, 2022, 36(10): 5403-5414.

[25] Jany P, Bejan A. Scaling theory of melting with natural convection in an enclosure[J]. International Journal of Heat and Mass Transfer, 1988, 31(6): 1221-1235.

基金资助

国家自然科学基金项目(11872318)