THE NUMERICAL SIMULATION OF MELTING PROCESS IN A LATERAL HEATED CAVITY UNDER TRANSVERSE MAGNETIC FIELDS
-
摘要: 磁场下的固−液相变过程在电磁冶金和增材制造等工程应用中广泛存在, 其中的熔化过程和流动机理尚未完全探究清楚. 方腔熔化是研究固−液相变过程的基础模型, 具有良好的普适性, 研究磁场对其流动的影响可以为其他复杂相变过程提供参考. 本文基于焓方法开展了固−液相变的数值模拟研究, 得到了垂直主环流方向的横向磁场对侧壁加热方腔中流动、传热和熔化过程的影响. 首先, 对于无磁场时的方腔熔化问题, 通过与已有的实验结果和数值结果进行对比, 证实了文献中方腔宽度对固−液界面的形状及位置影响不能忽视的结论. 随后, 对小磁场情况下的三维工况进行了直接数值模拟, 发现此时磁场效应主要表现为对混乱三维流动产生整流作用, 使流动趋于二维化. 但由于固−液界面的存在, 主流区的速度在趋于一致的同时也会反作用于界面, 其形状随磁场增强而逐渐转变为二维结构. 最后, 本文采用准二维模型分析了更强磁场时的情形, 讨论了不同参数对传热效率及界面形状的影响, 并发现了横向磁场作用下的垂直最大速度仍满足磁对流中的无量纲参数标度律关系.Abstract: The solid-liquid phase change process in the presence of the magnetic field has extensive and important applications in industrial engineering, i.e. the electromagnetic metallurgy and additive manufacturing, where the melting process and flow mechanisms have not been fully explored. The cavity melting model, as a basic problem, has good universality to study the solid-liquid phase change process, and it can also provide some basic information for the magnetohydrodynamics influence on the phase change problems. In this paper, the solid-liquid phase change process is simulated based on the enthalpy method, and the physical model of the cavity heated from the left wall is considered, with a transverse magnetic field perpendicular to the main circulation direction. The flow field, heat transfer and melting processes are investigated, focusing to the influences of different factors. At first, for the melting problem of a square cavity without a magnetic field, we compared our numerical solution with the experimental results reported by some other literatures, and it is confirmed that the influence of the cavity width on the profile and position of the solid-liquid interface cannot be ignored. Subsequently, we use the 3D model to simulate the cases in the presence of small magnetic fields, and it is found that the Lorentz force mainly acts as a rectification effect on the chaotic 3D flow, making the flow tend to be quasi-2D (Q2D). Meanwhile owing to the existence of the solid-liquid interface, the velocity field in the mainstream region tends to be more uniform under the external magnetic field, and thus the shape of the solid-liquid interface also transforms into the 2D structure correspondingly. Finally, the Q2D model is used to study the cases under greater magnetic fields, and the influences of different parameters on the heat transfer efficiency and interface shapes are discussed. Moreover, the scaling law to describe the relation between the max vertical velocity and other dimensionless parameters is also proposed, in order to quantitatively characterize the melting process under the transverse magnetic field.
-
引 言
磁场由于其非接触性以及对流动的调控能力, 被广泛应用于包含导电材料的固−液相变过程中, 如电磁冶金[1-2]、晶体生长[3-4]、增材制造[5-6]等等. 而方腔熔化模型作为研究磁场作用下的固−液相变问题的基础问题, 通常作为研究这类问题的一种重要方式. 研究磁场作用下的方腔熔化问题, 对相关工程实际应用具有重要的指导意义.
首先, 热源作用下的方腔熔化问题作为一种典型工况, 是探究相变材料熔化现象的一种有效手段, 已经有学者开展了众多实验和数值模拟方面的研究[7-11]. 其中, Zhu等[12]通过数值模拟计算了侧壁加热的方腔熔化问题, 发现了在三维和二维情形下固−液界面存在显著区别. Hu等[13]通过尺度分析和数值模拟验证, 发现只有
$ Pr\ll 1 $ 且长宽比远小于1的情况下, 三维方腔熔化问题才可以简化为二维问题. 并且文中给出的熔化体积分数标度关系对二维和三维结果给出了良好的预测.近年来, 由于电磁冶金行业的发展, 各种类型的电磁场在冶金工艺中的应用日益广泛. 磁场对液态金属的凝固/熔化的影响得到了越来越多的关注. 对于外加磁场下的方腔熔化问题, Dulikravich等[14]考虑了三维不可压缩黏性流体在磁场和不同重力条件下的凝固和熔化过程, 发现磁场的强度和方向会显著影响流动的速度和涡量. Veilleux等[15]和Zhang等[16]通过数值模拟和实验研究了低重力环境下的侧壁加热方腔熔化问题. 他们发现磁场可以用来模拟在实际的低重力环境中发现的关键熔化特征. 文献[17]通过数值模拟研究了在低温侧壁面作用下三维方腔内的自然对流和凝固过程. 发现不同方向磁场对流动和传热的影响存在显著区别, 尤其是平行温度梯度方向的磁场对流动的制稳作用最强. Feng等[18]采用LBM方法计算了沿平面内不同方向磁场作用下底面加热二维方腔内的固态镓(Ga)的熔化. 研究发现小
$ Ha $ 下($Ha={B}_{0}H\sqrt{\dfrac{{\sigma }_{{\rm{e}}}}{\mu }}$ , 用于描述洛伦兹力和黏性力的比值, 其中$ {B}_{0} $ 表示磁场强度,$ H $ 表示方腔尺寸,${\sigma }_{{\rm{e}}}$ 表示金属电导率,$ \; \mu $ 表示流体黏性), 竖直磁场和水平磁场对熔化过程的阻碍作用相似; 而在大$ Ha $ 下, 在熔化前期竖直磁场的阻碍作用更强而在熔化后期水平磁场的阻碍作用更强. Ghalambaz等[19]从实验和数值角度详细讨论了洛伦兹力对倾斜方腔内熔化流体流动的抑制作用. 文献[20]考虑了方腔底面中心方形截面条状热源加热熔化的模型, 讨论了大磁场强度下洛伦兹力对流动及传热的抑制作用. 以上的讨论主要关注均匀磁场对固−液相变的影响, 对于其他类型磁场(如移动磁场[21-22]、四极磁场[23]、非均匀磁场[24]等)也有相应的文献研究.此外, 对于不同方向磁场的影响, 在MHD(magnetohydrodynamics)自然对流中也有许多文献进行过研究分析. 文献[25]通过实验测量了不同加热速率和不同方向、大小磁场作用下熔融金属镓的自然对流传热效率. 揭示了垂直于环流方向的磁场对流动的阻滞作用明显弱于另外两个方向. 文献[26-27]分别进行了不同磁场方向下的三维自然对流数值模拟并发现了施加垂直环流磁场能够提高传热效率. 同时发现此方向磁场对速度的阻碍作用呈
$ Gr/Ha $ ($ Gr=\dfrac{g\beta \mathrm{\Delta }T{H}^{3}}{\kappa \nu } $ , 用于描述自然对流浮力和黏性力之比, 其中$ g $ 为重力加速度,$ \beta $ 为热膨胀系数,$ \mathrm{\Delta }T $ 为冷热源温差,$ \kappa $ 为热扩散系数,$ \nu $ 为运动黏性系数. )的量级而非其他情形下的$ Gr/H{a}^{2} $ 的量级. Chen等[28]在此基础上, 运用文献[29]提出的准二维模型与三维模型计算结果进行对比, 发现垂直环流方向磁场作用下, 在小$ Ha $ 时二次流动因受洛伦兹力的抑制而减小, 从而使得主流区的速度增强; 在大$ Ha $ 下三维模型与准二维模型的计算结果接近一致, 并且通过准二维模型计算验证了磁场对流动和传热的阻碍作用的标度关系.本工作研究垂直主环流方向的横向磁场影响, 并通过加热左侧壁面来驱动方腔熔化产生的自然对流. 本文仅考虑Rayleigh数
$Ra\leqslant {10}^{6}$ 的情况($Ra= \dfrac{g\beta \mathrm{\Delta }T{H}^{3}}{{\nu }^{2}}$ , 表示热传导和热对流之间的比值), 此时方腔内存在流动不稳定性但仍未出现湍流, 故壁面湍流的复杂情形[30]不在本文的讨论范围内. 另外熔化过程中可能存在杂质颗粒在流体中的流动, 类似于椭球的低速流动[31], 由于其对流动整体影响较小故本文亦不做考虑. 由于磁场的存在, 熔融金属流体的环形流动将产生感应电流和洛伦兹力, 使得主流区的速度受到抑制并沿磁场方向变得更加均匀. 然而, 当磁场足够强时, 垂直磁场方向的两侧壁面附近的剪切层具有更大的速度梯度, 一般将此剪切层称为Hartmann层. 在大$ Ha $ 下熔融液态金属的自然对流三维数值模拟中, 若要精确解析Hartmann层内的电流分布, 需要非常密集的网格和大量的计算资源, 此时采用准二维模型进行等效模拟可以大大提高计算效率. 本文对于固−液相变的模拟采用焓方法处理[32]而对于感应电势采用相容守恒方法处理[33]. 通过三维数值模拟发现小$ Ha $ 下增大磁场会使得自然对流增强和熔化速率加快, 同时流动及界面形状趋向于准二维. 此外, 本文采用准二维模型分析了大$ Ha $ 下, 磁场大小和自然对流强度对方腔内的流动形态、熔化速率以及固−液界面形状的影响. 并且通过尺度分析得到了大$ Ha $ 下, 熔化速率和垂直最大流速随不同$ Ra $ 和$ Ha $ 的变化关系.1. 数值模型
1.1 控制方程
本文所研究的横向磁场作用下侧壁加热熔化的三维方腔计算模型如图1所示. 立方体方腔边长为
$ H $ , 各壁面均为绝缘壁面. 左侧高温壁面温度为$ {T}_{{\rm{h}}} $ , 方腔内固相初始温度和右侧壁面温度均为$ {T}_{{\rm{i}}} $ 略低于熔点温度$ {T}_{{\rm{m}}} $ . 重力方向竖直向下, 均匀磁场$ {B}_{0} $ 垂直于主环流且方向沿$ z $ 轴正向.表 1 物性参数表Table 1. Material physical propertiesPhysical properties of gallium Units Values density (liquid) $ {\rho }_{\text{l}} $ $ {\text{kg/m}}^{\text{3}} $ 6095 density (solid) $ {\rho }_{\text{s}} $ $ {\text{kg/m}}^{\text{3}} $ 5907 specific heat ${c}_{{\rm{p}}}$ ${ {\rm{J} }/}({{\rm{kg}}\cdot {\rm{K}}}^{-1})$ 397.6 thermal expansivity $ \beta $ ${{\rm{K} } }^{-1}$ 1.27 × 10−4 dynamic viscosity $ \mu $ $ \text{kg/}\left(\text{m}\cdot \text{s}\right) $ 1.92 × 10−3 thermal conductivity $ \lambda $ $ \text{W/}\left(\text{m}\cdot \text{K}\right) $ 31.24 electrical conductivity (liquid) $ {\sigma }_{\text{e,l}} $ $ \text{S/m} $ 3.85 × 106 electrical conductivity (solid) $ {\sigma }_{\text{e,s}} $ $ \text{S/m} $ 6.64 × 106 thermal diffusivity (liquid) $ {\kappa }_{\text{l}} $ $ {\text{m}}^{\text{2}}\text{/s} $ 1.29 × 10−5 thermal diffusivity (solid) $ {\kappa }_{\text{s}} $ $ {\text{m}}^{\text{2}}\text{/s} $ 1.33 × 10−5 latent heat $ L $ $ \text{J/kg} $ 80160 本文采用液态金属镓(Ga)作为工质, 材料属性参数如表1所示. 该液态金属流体可视为单一组分的牛顿流体, 其流动状态为非定常不可压缩层流流动. 由于金属镓的固−液两相在熔点温度附近的密度差很小(见表1, 约为3%), 因此除浮力项采用Boussinesq假设来表征密度差引起的对流外, 方程中其他位置中的密度均视为常数
$ {\rho }_{\text{l}} $ , 数值模拟过程中不考虑相变过程的体积膨胀. 注意此简化在前人的数值模拟中已广泛采用[15,18]. 由于液态金属流动的低磁雷诺数属性(${Re}_{{{\rm{m}}}}={\mu }_{{\rm{m}}}{\sigma }_{{\rm{e}}}Lu\ll 1$ ,$\; {\mu }_{{\rm{m}}}$ 表示导磁率), 导电流体流动所产生的感应磁场亦无需考虑. 和壁面热源相比, 由焦耳耗散所产生的热量可以忽略不计. 综上, 本文中采用的无量纲形式的控制方程组如下$$ \nabla \cdot {{{\boldsymbol{u}}}} = 0 $$ (1) $$ \begin{split} & \frac{{\partial {{{\boldsymbol{u}}}}}}{{\partial t}} + \nabla \cdot \left( {{{{\boldsymbol{uu}}}}} \right) = - \nabla p + Pr {\nabla ^2}{{{\boldsymbol{u}}}} + RaPr \theta {{{{\boldsymbol{e}}}}_y}+ \\ & \qquad A{{{\boldsymbol{u}}}} + H{a^2}Pr \left( {{{{\boldsymbol{j}}}} \times {{{{\boldsymbol{e}}}}_B}} \right) \end{split} $$ (2) $$ {{{\boldsymbol{j}}}} = - \nabla \varPhi + {{{\boldsymbol{u}}}} \times {{{{\boldsymbol{e}}}}_B} $$ (3) $$ \nabla \cdot {{{\boldsymbol{j}}}} = 0 $$ (4) $$ \frac{{\partial \theta }}{{\partial t}} + {{{\boldsymbol{u}}}} \cdot \nabla \theta = {\nabla ^2}\theta - Ste\frac{{\partial f}}{{\partial t}} $$ (5) 其中,
${{\boldsymbol{u}}}$ ,${{\boldsymbol{j}}}$ ,$\varPhi$ ,$ \theta $ 分别表示无量纲形式的速度, 电流密度, 电势和温度. 方程中各变量对应的特征尺度分别为长度$ H $ , 质量$ \rho {H}^{3} $ , 速度$ \kappa /H $ , 时间$ {H}^{2}/\kappa $ , 磁场强度$ {B}_{0} $ , 电势$ \kappa {B}_{0} $ , 电流密度$ {\sigma }_{\text{e}}\kappa {B}_{0}/H $ , 相应地, 各无量纲数的定义分别为: Fourier数$Fo=\dfrac{\kappa t}{{H}^{2}}$ , 用于描述非稳态导热过程的无量纲时间, Prandtl数$Pr=\dfrac{\nu }{\kappa }$ , 用于描述动量和热扩散系数之比, Stefan数$Ste=\dfrac{{c}_{p}\mathrm{\Delta }T}{L}$ , 用于描述显热变化量与相变潜热的比值, Hartmann数Ha如前面定义.需要说明的是, 对于相变过程, 本文将采用焓方法进行模拟, 将固−液两相视为单一相流体共同计算温度场和速度场. 对于温度场, 利用工质内总比焓与显热和相变(熔化)潜热的关系式
$ h={c}_{p}\mathrm{\Delta }T + \mathrm{\Delta }H $ , 其中$ \mathrm{\Delta }H=fL $ ,$ f $ 为液相所占的体积分数,$ L $ 为相变潜热. 并将该关系式代入有量纲形式的能量方程, 可将方程改写为$$ \begin{split} & \rho \left( {\frac{{\partial h}}{{\partial t}} + {{{\boldsymbol{u}}}} \cdot \nabla h} \right) = \lambda {\nabla ^2}T \Rightarrow \\ &\qquad \frac{{\partial T}}{{\partial t}} + {{{\boldsymbol{u}}}} \cdot \nabla T = \kappa {\nabla ^2}T - \frac{L}{{{c_p}}}\frac{{\partial f}}{{\partial t}} \end{split} $$ (6) 其中, 对流项中关于相变潜热的分项只在界面附近有效, 且界面附近速度值很小固可将该项忽略. 而对于速度场, 则采用Carman-Kozeny关系式
$A=-\dfrac{C{\left(1-f\right)}^{2}}{{f}^{3} + \varepsilon }$ 对固相以及混合网格施加一个等效黏性, 使得速度趋近于0[34]. 其中$ C $ 为Carman-Kozeny常数, 其值取为$ {10}^{6} $ ,$ \varepsilon $ 为防止分母除0的计算常数, 取为$ {10}^{-3} $ .1.2 准二维模型
根据文献[29]得到的准二维模型, 可以得到考虑固液相变的控制方程组的无量纲形式如下
$$ {\nabla _ \bot } \cdot {{{{\boldsymbol{u}}}}_ \bot } = 0 $$ (7) $$ \begin{split} & \frac{{\partial {{{{\boldsymbol{u}}}}_ \bot }}}{{\partial t}} + {\nabla _ \bot } \cdot \left( {{{{{\boldsymbol{u}}}}_ \bot }{{{{\boldsymbol{u}}}}_ \bot }} \right) = - \nabla p + Pr {\nabla _ \bot }^2{{{{\boldsymbol{u}}}}_ \bot }+ \\ &\qquad Ra Pr \theta {{{{\boldsymbol{e}}}}_y} + A{{{{\boldsymbol{u}}}}_ \bot } - 2Ha Pr {{{{\boldsymbol{u}}}}_ \bot } \end{split} $$ (8) $$ \frac{{\partial \theta }}{{\partial t}} + {{\boldsymbol{u}}_ \bot } \cdot {\nabla _ \bot }\theta = {\nabla _ \bot }^2\theta - Ste\frac{{\partial f}}{{\partial t}} $$ (9) 上述方程中, 下标
$ \perp $ 表示该方程的求解域为主流区中任一垂直于磁场方向的截面. 通常来说, 当$ Ha\gg 1 $ 时, Hartmann层的厚度$ H/Ha $ 显著小于黏性边界层厚度, 层内的黏性耗散也由于速度梯度的增大而增强并可将其等效为对主流区流动的线性阻力(如式(8)的最后一项所示). 同时, 主流区的流场沿磁场方向趋于均匀,$ {H}{a}\gg 1 $ 流动由三维状态转变为准二维状态.1.3 边界条件
对于三维和准二维模型, 各壁面速度均采用无滑移无穿透边界条件. 左右两侧壁面温度采用Dirichlet边界条件, 其余壁面温度采用绝热边界条件. 各壁面的电势采用绝缘边界条件. 横向均匀磁场产生的洛伦兹力项作为体积力作用在全场. 对于三维模型, 洛伦兹力的计算需要通过先行求解关于
$ \varPhi $ 的电势泊松方程来得到. 而对于准二维模型, 此时洛伦兹力项作为线性项可直接代入动量方程计算.2. 数值算法及验证
本文中的数值模拟研究基于开源自适应计算平台Basilisk[35], N-S方程的计算采用两步近似投影法[36], 电势泊松方程的计算采用相容守恒格式. Hartmann层厚度为
$ H/Ha $ 且在本文计算中至少包含了6个网格, 以精确解析电流边界层. 通过网格无关性验证, 可以发现为使液相体积分数趋于收敛, 主流区网格尺寸应满足$\Delta r \leqslant H / 128$ , 固−液界面附近网格尺寸应满足$\Delta r \leqslant H / 256$ .2.1 二维/三维无磁场侧壁加热方腔熔化验证
本文首先验证了固−液相变过程模拟的精确性. 图2比较了三维矩形腔中, 左侧壁面加热熔化作用下, 不同时刻的固−液界面形状以及位置. 该工况中矩形腔的长高比为
$ D/H=1.4 $ , 宽高比为$ W/H=0.6 $ , 各无量纲参数为$ Pr=0.024\;4 $ ,$ Ra=6.6\times {10}^{5} $ , 以及$Ste= 0.04$ . 本文模拟的三维结果与文献[7]的实验测量结果和文献[12]的三维模拟结果符合良好. 同时还针对二维模拟结果与文献[12]的结果进行了比较, 发现由三维和二维模型计算得到的界面形状存在显著差异, 这种差异主要是由三维情形下沿第三个方向产生的二次流动造成的.![]()
2.2 三维MHD方腔自然对流验证
随后, 本文验证了磁场作用下三维方腔内的自然流动问题. 图3展示了
$ Ra=1\times {10}^{5} $ 时, 不同$ Ha $ 下的平均竖直速度($\overline{{u}_{y}}=\dfrac{1}{{t}_{\mathrm{t}\mathrm{o}\mathrm{t}\mathrm{a}\mathrm{l}}}{\displaystyle\int }_{0}^{{t}_{\mathrm{t}\mathrm{o}\mathrm{t}\mathrm{a}\mathrm{l}}}{u}_{y}\mathrm{d}t$ )在靠近加热壁面中线($ x=0.05 $ ,$ y=0.5 $ )的分布情况, 各$ Ha $ 下的结果均与文献[28]的模拟结果符合.![]()
图 3$ Ra=1\times {10}^{5} $ , 不同$ Ha $ 下位于$ x=0.05 $ ,$ y=0.5 $ 线段处的平均竖直速度($ \overline{{u}_{y}} $ )分布图. 点状线为文献[28]的计算结果,实线为本文的计算结果Figure 3. The distribution of mean vertical velocity ($ \overline{{u}_{y}} $ ) along$ x=0.05 $ ,$ y=0.5 $ for$ Ra=1\times {10}^{5} $ with different$ Ha $ . Dotted lines are the results simulated by Ref. [28], solid lines are the present numerical results此外, 本文还比较了
$ Ra=1\times {10}^{5} $ 时, 三维和准二维情形下不同$ Ha $ 下的时空平均$ Nu $ , 如表2所示 ($\overline{Nu}=\dfrac{1}{{t}_{\mathrm{t}\mathrm{o}\mathrm{t}\mathrm{a}\mathrm{l}}}\dfrac{1}{A}{\displaystyle\int }_{0}^{{t}_{\mathrm{t}\mathrm{o}\mathrm{t}\mathrm{a}\mathrm{l}}}{\displaystyle\int }_{0}^{l}{\displaystyle\int }_{0}^{l}\dfrac{\partial T}{\partial x}\mathrm{d}y\mathrm{d}z\mathrm{d}t$ ), 发现各个工况均与文献[28]的计算结果吻合良好. 由此充分证明了本文计算模型和算法的准确性.Table 2. Comparison of$ \overline{Nu} $ against the numerical results provided by Ref. [28] at different Ha3. 结果分析与讨论
3.1 小
$ \mathit{H}\mathit{a} $ 下的三维侧壁加热方腔熔化文献[26-27]分别进行了不同磁场方向下的三维自然对流数值模拟并发现了施加垂直环流磁场能够提高传热效率. 这是由于小
$ Ha $ 下洛伦兹力的整流作用造成的. 在此$ Ha $ 区间, 平行磁场方向的二次流动受到了抑制使得流动被迫向主环流发展, 同时Hartmann层中的黏性耗散没有明显增强, 流动仍然由浮升力主导. 本节中三维情形下的计算均以$ Ra={10}^{5} $ ,$ Ste= 0.05 $ 时的工况进行分析讨论.图4展示了
$ Fo=1.6, $ 在靠近加热壁面中线($ x= 0.05 $ ,$ y=0.5 $ )的竖直速度和竖直洛伦兹力分布. 当$ Ha=0 $ 时, 流动结构受到三维效应的影响变得混乱, 主环流的速度存在较大波动. 当$ Ha=50 $ 时, 靠近壁面位置的速度由于洛伦兹力的作用有所增大, 主流区的洛伦兹力接近于0, 此时流动受到的洛伦兹力在Hartmann层内的驱动作用大于在主流区中的焦耳耗散, 使得速度整体增大并且趋于均匀.![]()
图 4$ Fo=1.6, x=0.05 $ ,$ y=0.5 $ 处的瞬时竖直速度$ {u}_{y} $ (黑线)和洛伦兹力分量$ {Fl}_{y} $ (蓝线)的分布Figure 4. The distribution of the instantaneous vertical velocity (black lines) and y-component Lorentz force (blue line) at$Fo=1.6$ ,$ x=0.05 \;{\rm{and}}\; y=0.5 $ 图5对比展示了
$ Fo=1.6\,\mathrm{时}, Ha=0 $ ($z\leqslant 0.5$ )和$ Ha=50 $ ($z\geqslant 0.5$ )两工况中, 靠近壁面热源的平面($ x=0.05 $ )内的瞬时竖直速度云图和平面速度矢量. 从平面速度矢量的变化可以看出$ Ha=50 $ 时横向的二次流动明显受到了抑制, 从瞬时竖直速度云图可以发现$ Ha=50 $ 时主环流中的最大速度明显增大且分布更均匀. 这是由于洛伦兹力对Hartmann层流动的驱动使得液相的流动更多地被导向主环流方向从而使得横向二次流动减少, 因此整体上表现出对流场的整流效果. 为了确定整个流场中速度的变化情况, 本文进一步统计了全场总动能随时间的变化情况, 如图6所示. 可以看到各工况下总动能随时间均存在一定波动, 而随着液相区域的增大,$ Ha=50 $ 工况下的总动能明显大于$ Ha=0 $ 工况下的总动能, 这说明小$ Ha $ 下洛伦兹力的作用以Hartmann层内的驱动力为主导, 对流体做正功, 使得流动速度整体增大.![]()
图 5$ Fo=1.6, x=0.05 $ 处的瞬时竖直速度云图, 分别对应$ Ha=0 $ ($z\leqslant 0.5$ )和$ Ha=50 $ ($z\geqslant 0.5$ ). 图中矢量表示$ z-y $ 平面内的速度Figure 5. The contour maps of the instantaneous vertical velocity ($ {u}_{y} $ ) for$ Ha=0 $ ($z\leqslant 0.5$ ) and$ Ha=50 $ ($z\geqslant 0.5$ ) at$ x=0.05 $ . The vectors present the velocity in the$ z-y $ plane![]()
图 6$ Ha=0 $ 和$ Ha=50 $ 时, 全场总动能随时间的变化情况Figure 6. Evolution of total kinetic energy in the whole domain for$ Ha=0 $ and$ Ha=50 $ 为了确认小
$ Ha $ 下磁场对传热的增强作用, 比较了不同$ Ha $ 下液相体积分数和$ Nu $ 随时间的变化. 图7展示了$ Ha=0 $ ,$ Ha=25 $ 和$ Ha=50 $ 时液相体积分数随时间的变化, 进一步验证了小$ Ha $ 下洛伦兹力对流场结构的整流作用. 主流区速度受的洛伦兹力的抑制作用较小, 同时由于剪切边界层内洛伦兹力做正功加快了附近流体的速度从而使得主流区增大, 而因此出现传热和熔化增强的现象.![]()
图 7 不同$ Ha $ 下, (a)液相体积分数随时间的变化和(b)$ Nu $ 随时间的变化Figure 7. (a) Time evolution of the liquid volume fraction and (b) time evolution of$ Nu $ at different$ Ha $ 图8通过流线图对比
$ Ha=0 $ 和$ Ha=50 $ 时流线和固−液界面形状, 发现随着磁场增大, 流动逐渐趋于稳定和二维化, 同时, 从界面上云图的$ x $ 轴坐标值可以看到界面形状也趋于二维化. 因此采用准二维模型模拟大横向磁场下自然对流熔化问题是有效且可靠的.![]()
图 8$ Fo=1.6 $ , 界面形状及流线图比较, 界面上云图的颜色表示$ x $ 轴坐标值Figure 8. Comparison of the solid-liquid interface and the streamlines at$Fo=1.6 $ The contours on the interface are the$ x $ coordinate3.2 大
$ \mathit{H}\mathit{a} $ 下的准二维侧壁加热方腔熔化由上文可知, 随着
$ Ha $ 的增大, 方腔内熔融流体的流动形态逐渐从混乱(不稳定)的三维流动转变为稳定的准二维环流. 同时, 固−液界面也从中间凸起的三维界面转变为准二维界面. 运用提出的准二维模型, 本节分别对$ Ste=0.05 $ 时,${10}^{4}\leqslant Ra\leqslant {10}^{6}$ 和$100\leqslant Ha\leqslant 6400$ 范围内的方腔熔化自然对流过程中流动、传热和界面形状进行了分析.图9展示了
$ Ha=400 $ 时, 不同时刻和$ Ra $ 下界面的形状、流场和温度场. 可以看到$ Fo=0.2 $ 时, 熔融流体区域所占宽度只有方腔尺寸的15%, 此时流场对应的瞬时$ Ra $ 很小, 传热和熔化受热传导主导. 而其中$ Ra=1\times {10}^{6} $ 时已经出现4个涡流结构使得对流换热的作用开始增强. 而对于$ Fo=1 $ 时,$ Ra $ 较小的两个工况中均出现单一主环流, 而$ Ra= 1\times {10}^{6} $ 中的流动由于受到浮力的驱动作用更强仍然表现出主环流和多个小环流共存的情况. 当$ Fo=3 $ 时, 对流作用对界面形状的影响更为明显.图10展示了
$ Ra=1\times {10}^{6} $ 时, 固−液界面随$ Ha $ 增大的变化趋势, 所取时刻均为$ Fo=2 $ . 可以看到当$ Ha $ 逐渐增大时, 由于自然对流受到的抑制增强从而使得固−液界面形态逐渐由大曲率界面向平直界面转化.![]()
图 9$ Ha= $ 400时, 不同$ Ra $ 下的温度云图及流线图比较(从左至右依次为$ Ra=1\times {10}^{4} $ ,$ Ra=1\times {10}^{5}, \;Ra=1\times {10}^{6} $ )Figure 9. Comparison of the temperature contour and the streamlines at different instants ($ Ra=1\times {10}^{4} $ ,$ Ra=1\times {10}^{5},\;Ra=1\times {10}^{6} $ from left to right) at$ Ha= $ 400![]()
9 Ha = 400时, 不同Ra下的温度云图及流线图比较(从左至右依次为$Ra=1\times {10}^{4} $ ,$Ra=1\times {10}^{5},\; Ra=1\times {10}^{6} $ )(续)9. Comparison of the temperature contour and the streamlines at different instants ($Ra=1\times {10}^{4} $ ,$Ra=1\times {10}^{5},\;Ra=1\times {10}^{6} $ from left to right) at$Ha= $ 400 (continued)![]()
图 10$ Ra=1\times {10}^{6} $ ,$ Fo=2 $ 固−液界面随$ Ha $ 增大的变化趋势Figure 10. Comparison of the solid-liquid interfaces at$ Ra=1\times {10}^{6} $ and$ Fo=2 $ when Ha is varied图11展示了在不同Fo, Ha和Ra下, 熔融金属液相体积分数f和传热Nu的演化情况. 如图11(a)所示, 磁场强度对熔化速度的影响与
$ Ra $ 有关. 当$ Ra=1\times {10}^{4} $ 时, 液相在熔化过程中受到的浮升力过小, 熔化界面随时间的变化与一维热传导的解析解(图中空心点)的结果基本一致, 说明此时熔化完全由热传导主导. Hu等[13]通过尺度分析得到热传导占主导时的标度关系为$ {f}_{\mathrm{l}}\sim {\left(SteFo\right)}^{\frac{1}{2}} $ , 对流换热占主导时的标度关系$ {f}_{\mathrm{l}}\sim SteFo $ . 因此当$ Ste $ 为常数时,$ {f}_{\mathrm{l}} $ 与$ Fo $ 之间仍然应当满足上述标度关系. 当$Ra=1\times {10}^{5}$ 和$ Ra=1\times {10}^{6} $ 时, 由于液相中对流作用的增强使得在$ Ha $ 较小时仍能在对流阶段保持液相体积分数与$ Fo $ 的一次方关系, 而随着$ Ha $ 逐渐增大, 凝固速率由于洛伦兹力对主环流的阻滞作用逐渐下降, 并分别在$ Ha=3200 $ 和$ Ha=12\;800 $ 时回复到由热传导主导的液相体积分数与$ Fo $ 的$ 1/2 $ 次方关系.图11(b)展示了不同
$ Ra $ 下, 磁场强度对$ Nu $ 的影响. 当$Ra=1\times {10}^{4}$ 时, 不同$ Ha $ 下$ Nu $ 基本没有变化, 当$Ra=1\times {10}^{5}$ 和$ Ra=1\times {10}^{6} $ 时, 由于自然对流随$ Ha $ 而减弱, 使得$ Nu $ 逐渐减小. 值得注意的是当$Ra= 1\times {10}^{6}$ 时, 尽管流动已趋于二维, 在小$ Ha $ 下平面内的强烈且不稳定的自然对流仍然会引起壁面处$ Nu $ 的震荡.![]()
图 11 不同$ Ha $ 下, (a)液相体积分数和(b)左侧壁面$ Nu $ 随时间的变化, 从左至右依次为$ Ra=1\times {10}^{4} $ ,$ Ra=1\times {10}^{5}, \;Ra=1\times {10}^{6} $ Figure 11. Time evolution of (a) the liquid volume fraction and (b) the Nusselt number at different Hartmann number, while the sub-panels correspond to$ Ra=1\times {10}^{4} $ ,$ Ra=1\times {10}^{5} $ and$ Ra=1\times {10}^{6} $ from left to right如Tagawa等[27]的文中所分析的, 对于磁场方向垂直主环流方向的MHD问题, 侧层中的电流是用来表征洛伦兹力的阻滞作用的一个重要方式. 并且由电流与速度的关系可以得知方腔侧壁剪切层内的最大速度应在
$ Ha/Gr $ , 即$ HaPr/Ra $ 的量级. 图12展示了熔化过程中熔融液体的最大竖向速度$ {u}_{y,\mathrm{m}\mathrm{a}\mathrm{x}} $ 与无量纲参数组合$ HaPr/Ra $ 的变化关系. 可以看到在$ Ha $ 足够大的情况下, 本文计算中所考虑的不同$ Ha $ 及$ Ra $ 下各工况中的$ {u}_{y,\mathrm{m}\mathrm{a}\mathrm{x}} $ 均满足同一标度关系.![]()
图 12 熔化过程中的最大速度($ {u}_{y,\mathrm{m}\mathrm{a}\mathrm{x}} $ )与无量纲参数组合$ HaPr/Ra $ 的关系Figure 12. The relation between the maximum velocity (${u}_{y,{\rm{max}}}$ ) during the melting process and the combination of dimensionless numbers of$ HaPr/Ra $ 4. 结 论
本文通过数值模拟的方法研究了垂直环流方向横向磁场对侧壁加热方腔内熔化和自然对流的影响, 并分别对小
$ Ha $ 和大$ Ha $ 两个区间的工况采用三维模型和准二维模型进行计算和分析, 得到了以下主要结论.本文采用三维数值模拟分析了小
$ Ha $ 下的情形. 在该区间下的磁场强度较小, 洛伦兹力对流动的影响有限, 其对流动的影响主要表现为对流动结构的整流作用, 而熔化和流动过程仍然受浮升力主导. 随着磁场增大, 二次流动受到抑制, 且主环流速度有所增大, 从而出现熔化和传热速率增大的情况, 同时熔化过程中固−液界面的存在也使得流动形态进一步趋向于二维化, 而二维的流动形态也反过来决定了固−液界面也将趋于二维化.对于大
$ Ha $ 下的情形, 随着磁场增大, 方腔内的流动形态和界面形状均趋于二维化. 经验证在此区间准二维模型的计算结果与三维模型非常接近, 采用准二维模型计算能在保证准确的前提下大大缩短计算时间. 可以发现不同$ Ra $ 下的对流强度都受到洛伦兹力的抑制而显著减小, 固−液界面形态也趋于平直. 液相的熔化速率和壁面处的传热效率都呈单调下降趋势. 而熔化过程中的最大速度仍然满足Tagawa等[27]提出的无量纲参数组合$ HaPr/Ra $ 之间的标度关系. -
![]()
![]()
图 3
$ Ra=1\times {10}^{5} $ , 不同$ Ha $ 下位于$ x=0.05 $ ,$ y=0.5 $ 线段处的平均竖直速度($ \overline{{u}_{y}} $ )分布图. 点状线为文献[28]的计算结果,实线为本文的计算结果Figure 3. The distribution of mean vertical velocity (
$ \overline{{u}_{y}} $ ) along$ x=0.05 $ ,$ y=0.5 $ for$ Ra=1\times {10}^{5} $ with different$ Ha $ . Dotted lines are the results simulated by Ref. [28], solid lines are the present numerical results![]()
图 4
$ Fo=1.6, x=0.05 $ ,$ y=0.5 $ 处的瞬时竖直速度$ {u}_{y} $ (黑线)和洛伦兹力分量$ {Fl}_{y} $ (蓝线)的分布Figure 4. The distribution of the instantaneous vertical velocity (black lines) and y-component Lorentz force (blue line) at
$Fo=1.6$ ,$ x=0.05 \;{\rm{and}}\; y=0.5 $ ![]()
图 5
$ Fo=1.6, x=0.05 $ 处的瞬时竖直速度云图, 分别对应$ Ha=0 $ ($z\leqslant 0.5$ )和$ Ha=50 $ ($z\geqslant 0.5$ ). 图中矢量表示$ z-y $ 平面内的速度Figure 5. The contour maps of the instantaneous vertical velocity (
$ {u}_{y} $ ) for$ Ha=0 $ ($z\leqslant 0.5$ ) and$ Ha=50 $ ($z\geqslant 0.5$ ) at$ x=0.05 $ . The vectors present the velocity in the$ z-y $ plane![]()
图 6
$ Ha=0 $ 和$ Ha=50 $ 时, 全场总动能随时间的变化情况Figure 6. Evolution of total kinetic energy in the whole domain for
$ Ha=0 $ and$ Ha=50 $ ![]()
图 7 不同
$ Ha $ 下, (a)液相体积分数随时间的变化和(b)$ Nu $ 随时间的变化Figure 7. (a) Time evolution of the liquid volume fraction and (b) time evolution of
$ Nu $ at different$ Ha $ ![]()
图 8
$ Fo=1.6 $ , 界面形状及流线图比较, 界面上云图的颜色表示$ x $ 轴坐标值Figure 8. Comparison of the solid-liquid interface and the streamlines at
$Fo=1.6 $ The contours on the interface are the$ x $ coordinate![]()
图 9
$ Ha= $ 400时, 不同$ Ra $ 下的温度云图及流线图比较(从左至右依次为$ Ra=1\times {10}^{4} $ ,$ Ra=1\times {10}^{5}, \;Ra=1\times {10}^{6} $ )Figure 9. Comparison of the temperature contour and the streamlines at different instants (
$ Ra=1\times {10}^{4} $ ,$ Ra=1\times {10}^{5},\;Ra=1\times {10}^{6} $ from left to right) at$ Ha= $ 400![]()
9 Ha = 400时, 不同Ra下的温度云图及流线图比较(从左至右依次为
$Ra=1\times {10}^{4} $ ,$Ra=1\times {10}^{5},\; Ra=1\times {10}^{6} $ )(续)9. Comparison of the temperature contour and the streamlines at different instants (
$Ra=1\times {10}^{4} $ ,$Ra=1\times {10}^{5},\;Ra=1\times {10}^{6} $ from left to right) at$Ha= $ 400 (continued)![]()
图 10
$ Ra=1\times {10}^{6} $ ,$ Fo=2 $ 固−液界面随$ Ha $ 增大的变化趋势Figure 10. Comparison of the solid-liquid interfaces at
$ Ra=1\times {10}^{6} $ and$ Fo=2 $ when Ha is varied![]()
图 11 不同
$ Ha $ 下, (a)液相体积分数和(b)左侧壁面$ Nu $ 随时间的变化, 从左至右依次为$ Ra=1\times {10}^{4} $ ,$ Ra=1\times {10}^{5}, \;Ra=1\times {10}^{6} $ Figure 11. Time evolution of (a) the liquid volume fraction and (b) the Nusselt number at different Hartmann number, while the sub-panels correspond to
$ Ra=1\times {10}^{4} $ ,$ Ra=1\times {10}^{5} $ and$ Ra=1\times {10}^{6} $ from left to right![]()
图 12 熔化过程中的最大速度(
$ {u}_{y,\mathrm{m}\mathrm{a}\mathrm{x}} $ )与无量纲参数组合$ HaPr/Ra $ 的关系Figure 12. The relation between the maximum velocity (
${u}_{y,{\rm{max}}}$ ) during the melting process and the combination of dimensionless numbers of$ HaPr/Ra $ 表 1 物性参数表
Table 1 Material physical properties
Physical properties of gallium Units Values density (liquid) $ {\rho }_{\text{l}} $ $ {\text{kg/m}}^{\text{3}} $ 6095 density (solid) $ {\rho }_{\text{s}} $ $ {\text{kg/m}}^{\text{3}} $ 5907 specific heat ${c}_{{\rm{p}}}$ ${ {\rm{J} }/}({{\rm{kg}}\cdot {\rm{K}}}^{-1})$ 397.6 thermal expansivity $ \beta $ ${{\rm{K} } }^{-1}$ 1.27 × 10−4 dynamic viscosity $ \mu $ $ \text{kg/}\left(\text{m}\cdot \text{s}\right) $ 1.92 × 10−3 thermal conductivity $ \lambda $ $ \text{W/}\left(\text{m}\cdot \text{K}\right) $ 31.24 electrical conductivity (liquid) $ {\sigma }_{\text{e,l}} $ $ \text{S/m} $ 3.85 × 106 electrical conductivity (solid) $ {\sigma }_{\text{e,s}} $ $ \text{S/m} $ 6.64 × 106 thermal diffusivity (liquid) $ {\kappa }_{\text{l}} $ $ {\text{m}}^{\text{2}}\text{/s} $ 1.29 × 10−5 thermal diffusivity (solid) $ {\kappa }_{\text{s}} $ $ {\text{m}}^{\text{2}}\text{/s} $ 1.33 × 10−5 latent heat $ L $ $ \text{J/kg} $ 80160 
下载: 导出CSV
-
[1] Battira M, Bessaih R. Three-dimensional natural convection in the horizontal Bridgman configuration under various wall electrical conductivity and magnetic field. Numerical Heat Transfer, Part A:Applications, 2008, 55(1): 58-76 doi: 10.1080/10407780802603113
[2] Hadid HB, Henry D. Numerical study of convection in the horizontal Bridgman configuration under the action of a constant magnetic field. Part 2. Three-dimensional flow. Journal of Fluid Mechanics, 1997, 333: 57-83
[3] Grants I, Pal J, Gerbeth G. Physical modelling of Czochralski crystal growth in horizontal magnetic field. Journal of Crystal Growth, 2017, 470: 58-65 doi: 10.1016/j.jcrysgro.2017.04.021
[4] Dennis BH, Dulikravich GS. Magnetic field suppression of melt flow in crystal growth//CHT'01-Advances in Computational Heat Transfer II. Proceedings of a Second Symposium. Begel House Inc. , 2001
[5] Wang L, Yao J, Hu Y, et al. Suppression effect of a steady magnetic field on molten pool during laser remelting. Applied Surface Science, 2015, 351: 794-802 doi: 10.1016/j.apsusc.2015.05.179
[6] Wang L, Yao J, Hu Y, et al. Influence of electric-magnetic compound field on the WC particles distribution in laser melt injection. Surface and Coatings Technology, 2017, 315: 32-43 doi: 10.1016/j.surfcoat.2017.01.116
[7] Gau C, Viskanta R. Melting and solidification of a pure metal on a vertical wall. Journal of Heat Transfer, 1986, 108(1): 174-181 doi: 10.1115/1.3246884
[8] Ben-David O, Levy A, Mikhailovich B, et al. 3 D numerical and experimental study of gallium melting in a rectangular container. International Journal of Heat and Mass Transfer, 2013, 67: 260-271 doi: 10.1016/j.ijheatmasstransfer.2013.07.058
[9] Dai R, Bian Q, Wang Q, et al. Evolution of natural convection melting inside cavity heated from different sides using enthalpy based lattice Boltzmann method. International Journal of Heat and Mass Transfer, 2018, 121: 715-725 doi: 10.1016/j.ijheatmasstransfer.2018.01.032
[10] Favier B, Purseed J, Duchemin L. Rayleigh–Bénard convection with a melting boundary. Journal of Fluid Mechanics, 2019, 858: 437-473 doi: 10.1017/jfm.2018.773
[11] Hu Y, Li D, Shu S, et al. Lattice Boltzmann simulation for three-dimensional natural convection with solid-liquid phase change. International Journal of Heat and Mass Transfer, 2017, 113: 1168-1178 doi: 10.1016/j.ijheatmasstransfer.2017.05.116
[12] Zhu W, Wang M, Chen H. 2 D and 3 D lattice Boltzmann simulation for natural convection melting. International Journal of Thermal Sciences, 2017, 117: 239-250 doi: 10.1016/j.ijthermalsci.2017.03.025
[13] Hu N, Fan LW, Zhu ZQ. Can the numerical simulations of melting in a differentially-heated rectangular cavity be rationally reduced to 2 D? A comparative study between 2 D and 3 D simulation results. International Journal of Heat and Mass Transfer, 2021, 166: 120751 doi: 10.1016/j.ijheatmasstransfer.2020.120751
[14] Dulikravich GS, Ahuja V, Seungsoo L. Modeling three-dimensional solidification with magnetic fields and reduced gravity. International Journal of Heat and Mass Transfer, 1994, 37(5): 837-853 doi: 10.1016/0017-9310(94)90120-1
[15] Veilleux DL, Faghri M, Asako Y, et al. Convection enhancement in melting by electromagnetic fields in a low-gravity environment: Side wall heating. Numerical Heat Transfer, Part A:Applications, 2007, 51(2): 129-158 doi: 10.1080/10407780600710425
[16] Zhang H, Charmchi M, Veilleux D, et al. Numerical and experimental investigation of melting in the presence of a magnetic field: Simulation of low-gravity environment. Journal of Heat Transfer, 2007, 129(4): 568-576 doi: 10.1115/1.2709961
[17] Bouabdallah S, Bessaїh R. Effect of magnetic field on 3 D flow and heat transfer during solidification from a melt. International Journal of Heat and Fluid Flow, 2012, 37: 154-166 doi: 10.1016/j.ijheatfluidflow.2012.07.002
[18] Feng Y, Li H, Li L, et al. Investigation of the effect of magnetic field on melting of solid gallium in a bottom-heated rectangular cavity using the lattice Boltzmann method. Numerical Heat Transfer, Part A:Applications, 2016, 69(11): 1263-1279 doi: 10.1080/10407782.2015.1127732
[19] Ghalambaz M, Doostanidezfuli A, Zargartalebi H, et al. MHD phase change heat transfer in an inclined enclosure: effect of a magnetic field and cavity inclination. Numerical Heat Transfer, Part A:Applications, 2017, 71(1): 91-109 doi: 10.1080/10407782.2016.1244397
[20] Bondareva NS, Sheremet MA. Natural convection heat transfer combined with melting process in a cubical cavity under the effects of uniform inclined magnetic field and local heat source. International Journal of Heat and Mass Transfer, 2017, 108: 1057-1067 doi: 10.1016/j.ijheatmasstransfer.2016.12.108
[21] Dadzis K, Lukin G, Meier D, et al. Directional melting and solidification of gallium in a traveling magnetic field as a model experiment for silicon processes. Journal of Crystal Growth, 2016, 445: 90-100 doi: 10.1016/j.jcrysgro.2016.03.037
[22] Ben-David O, Levy A, Mikhailovich B, et al. Impact of traveling permanent magnets on low temperature metal melting in a cuboid. International Journal of Heat and Mass Transfer, 2016, 99: 882-894 doi: 10.1016/j.ijheatmasstransfer.2016.04.017
[23] Dai R, Bian Q, Wang Y, et al. Lattice Boltzmann simulation for melting control through an extra magnetic quadrupole field. Numerical Heat Transfer, Part A:Applications, 2019, 75(4): 254-270 doi: 10.1080/10407782.2019.1583952
[24] Mehryan SAM, Tahmasebi A, Izadi M, et al. Melting behavior of phase change materials in the presence of a non-uniform magnetic-field due to two variable magnetic sources. International Journal of Heat and Mass Transfer, 2020, 149: 119-184
[25] Okada K, Ozoe H. Experimental heat transfer rates of natural convection of molten gallium suppressed under an external magnetic field in either the X, Y, or Z direction. Journal of Heat Transfer, 1992, 114(1): 107-114 doi: 10.1115/1.2911234
[26] Ozoe H, Okada K. The effect of the direction of the external magnetic field on the three-dimensional natural convection in a cubical enclosure. International Journal of Heat and Mass Transfer, 1989, 32(10): 1939-1954 doi: 10.1016/0017-9310(89)90163-4
[27] Tagawa T, Authié G, Moreau R. Buoyant flow in long vertical enclosures in the presence of a strong horizontal magnetic field. Part 1. Fully-established flow. European Journal of Mechanics-B/Fluids, 2002, 21(4): 383-398
[28] Chen L, Liu BQ, Ni MJ. Study of natural convection in a heated cavity with magnetic fields normal to the main circulation. International Journal of Heat and Mass Transfer, 2018, 127: 267-277 doi: 10.1016/j.ijheatmasstransfer.2018.06.019
[29] Sommeria JL, Moreau R. Why, how, and when, MHD turbulence becomes two-dimensional. Journal of Fluid Mechanics, 1982, 115: 507-518
[30] 郑晓静, 王国华. 高雷诺数壁湍流的研究进展及挑战. 力学进展, 2020, 50: 202001 Zheng X J, Wang G H. Progresses and challenges of high Reynolds number wall-bounded turbulence. Advances in Mechanics, 2020, 50: 202001（in Chinese）
[31] 丛成华, 邓小刚, 毛枚良. 绕椭球的低速流动研究. 力学进展, 2021, 51(3): 467-619 Cong C H, Deng X G, Mao M L. Advances in complex low speed flow around a prolate spheroid. Advances in Mechanics, 2021, 51(3): 467-619（in Chinese）
[32] Voller VR. Fast implicit finite-difference method for the analysis of phase change problems. Numerical Heat Transfer, 1990, 17(2): 155-169
[33] Zhang J, Ni MJ. Direct simulation of multi-phase MHD flows on an unstructured Cartesian adaptive system. Journal of Computational Physics, 2014, 270: 345-365 doi: 10.1016/j.jcp.2014.03.030
[34] Voller VR, Cross M, Markatos NC. An enthalpy method for convection/diffusion phase change. International Journal for Numerical Methods in Engineering, 1987, 24(1): 271-284 doi: 10.1002/nme.1620240119
[35] Popinet S. An accurate adaptive solver for surface-tension-driven interfacial flows. Journal of Computational Physics, 2009, 228(16): 5838-5866 doi: 10.1016/j.jcp.2009.04.042
[36] Popinet S. Gerris: a tree-based adaptive solver for the incompressible Euler equations in complex geometries. Journal of Computational Physics, 2003, 190(2): 572-600 doi: 10.1016/S0021-9991(03)00298-5
