基于OpenFOAM的高温熔融锡液水力碎化数值模拟研究

Numerical Simulation of Hydraulic Fragmentation of High-Temperature Molten Tin Based on OpenFOAM

摘要: 为揭示高温熔融锡液与水相互作用过程中的碎化机理，基于开源软件包OpenFOAM®-v1806下的求解器icoReactingMultiphaseInterFoam模拟熔融锡液入水后的水力碎化过程，分析瑞利–泰勒不稳定性、开尔文–亥姆霍兹不稳定性及界面剪切力对熔融锡液柱碎化过程作用过程，探究水力碎化对液柱直径及液柱入水初速度对碎化剧烈程度的影响。研究表明液柱碎化剧烈程度、液柱单位时间内入水深度与液柱直径成反比，液柱入水后碎化剧烈程度、锡凝固速度与液柱入水时初速度成正比。

Abstract: In order to reveal the mechanism of fragmentation during the interaction between high-temperature molten tin and water, the icoReactingMultiphaseInterFoam solver within Open-FOAM®-v1806 was used to simulate the hydraulic fragmentation process of molten tin entering water. The effects of Rayleigh-Taylor instability, Kelvin-Helmholtz instability and interfacial shear force on the fragmentation process of molten tin were analyzed. And the sensitivity of the degree of fragmentation to the diameter of the liquid column and the initial velocity of the liquid column entering the water was explored. Research shows that the degree of fragmentation of the liquid column and underwater penetration in unit time are inversely proportional to the diameter of the liquid column. Meanwhile, the degree of fragmentation and the formation rate of solid tin are proportional to the initial speed of the liquid column entering water.

1. 引言

在工业生产中水作为最常见的冷却剂，与高温熔融金属接触反应过程中常伴随剧烈的传热传质过程、熔融金属碎化过程以及多种形式的能量转化过程。在高温熔融金属遇水碎化过程数值模拟研究中，周源 ‎[1] 认为水力学碎化过程发生在压力峰值的两端。袁明豪 ‎[2] 基于熔融液滴水力碎化模拟研究发现液滴变形与液滴表面的水力不稳定性是液滴碎化的主要机理。杨志 ‎[3] 利用MPS数值模拟程序将熔滴碎化区域分为凝固变形区、破碎快速增长区、破碎缓慢增长区，碎化程度随温度及Weber数的增长呈先增后减的趋势。钟明君 ‎[4] 根据METRIC程序对熔融物水力学碎化行为及其界面特性进行数值计算，研究表明在高Weber数下液柱碎化主要受侧面不稳定波和剥离效应主导。沈正祥等 ‎[5] 采用有限体积法与水平集(Level-set)方法捕捉低速射流碎化分解过程中到界面的不稳定发展、变形、流场内部速度场及压力场分布等；研究获悉当入射速度较低时，流场内轴向速度随入射深度增加而变大，射流表面碎化范围内径向速度变化比较剧烈，对应的压力场为负压，射流的碎化长度与雷诺数(Re/We^0.5)呈指数型关系。

为揭示熔融金属与水相互作用机理，前人多利用计算流体动力学相关软件，如ANSYS CFD、ANSYS CFX、MPS、MC3D等来进行数值模拟研究。其中，OpenFOAM软件运用较少，但其作为一种开源软件，可编译多种符合实际应用的求解器，用于计算分析不可压缩或可压缩流体的相变换热过程、燃烧过程、浮力驱动过程、固体应力分析等。本研究选取金属锡为研究材料，基于开源软件包OpenFOAM®-v1806下的icoReactingMultiphaseInterFoam求解器模拟熔融锡液柱入水后的水力碎化过程，分析瑞利–泰勒(Rayleigh-Taylor, R-T)不稳定性、开尔文–亥姆霍兹(Kelvin-Helmholtz, K-H)不稳定性及界面剪切力对熔融锡液柱碎化过程作用结果，探究水力碎化对液柱直径及液柱入水初速度对碎化剧烈程度的影响。

2. 数值模拟方法

2.1. VOF模型

VOF模型对计算域内不互溶的多种流体求解同一动量方程，网格单元内各流体的体积分数由和为1的条件来约束。a_i为0或1时，分别表示网格内无该流体或均为该流体，介于之间则表示存在其他流体。

${\displaystyle \underset{i=1}{\overset{j}{\sum }}{\alpha }_i=1}$ (1)

因多相流系统内存在相变，以液相向汽相运输为例，相体积运输方程如下：

$\frac{\partial {\alpha }_{liq}}{\partial t}+\nabla \cdot \left(u{\alpha }_{liq}\right)=\stackrel{\dot{}}{m}\left[\frac{1}{{\rho }_{liq}}-{\alpha }_{liq}\left(\frac{1}{{\rho }_{liq}}-\frac{1}{{\rho }_{vap}}\right)\right]$ (2)

密度、粘度、比热容等相关物理量与网格内介质的体积分数相关，物性表达公式如下：

$\rho ={\displaystyle \underset{i=1}{\overset{j}{\sum }}{\alpha }_i}{\rho }_i$ (3)

$\mu ={\displaystyle \underset{i=1}{\overset{j}{\sum }}{\alpha }_i}{\mu }_i$ (4)

$\lambda ={\displaystyle \underset{i=1}{\overset{j}{\sum }}{\alpha }_i}{\lambda }_i$ (5)

$C_p={\displaystyle \underset{i=1}{\overset{j}{\sum }}{\alpha }_i}{C}_p{}_i$ (6)

式中，a为体积分数；t为时间，s；u为流体矢量流速，m/s；r为密度，kg/m³；m为动力粘度，kg/(m·s)；l为导热系数，W/(m·K)；C_p为定压比热容，J/(kg·K)；j为计算域内的相数；i为计算域内每一相对应的标号；下标liq、vap分别指液相、汽相。

在setFieldsDict中定义初始相域，所有单元均需被定义且同一域中组分不可叠加，要求满足下式：

${\displaystyle \underset{i=1}{\overset{j}{\sum }}{Y}_i}=1$ (7)

式中，Y为组分分数。

连续性方程是质量守恒定律在流体力学中的具体应用，使用前提是将流体看作连续介质，应用连续介质模型，认为速度和密度均为空间坐标和时间的连续可微函数 ‎[6] 。其中，不可压流体的连续性方程为：

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

Navier-Stokes方程是牛顿第二定律的运用，核心是动量守恒。因此，Navier-Stokes方程最基本的假设是把流体看作连续介质，从宏观方位分析时，可以忽略流体的随机分子运动，以描述流体的多个宏观物理量 ‎[7] 。结合连续性方程，考虑重力、表面张力两种体积力的守恒型不可压缩动量方程如下：

$\frac{\partial \left(\rho u\right)}{\partial t}+\nabla \cdot \left(\rho uu\right)=-\nabla p+\nabla \cdot \left[\mu \left(\nabla u+\nabla u^{\text{T}}\right)\right]+\rho g+f$ (9)

式中，p为正应力，Pa；g为重力加速度矢量，m/s²；f为通过连续表面应力模型(Continuum Surface Force, CSF)计算所得的表面张力，N/m。

由于相变过程中伴随着能量的传递，需要建立并求解能量方程。能量方程基于热力学第一定律，在OpenFOAM的多相流求解器中，不可压缩流体的能量守恒方程表达形式如下：

$\frac{\partial T}{\partial t}+\nabla \cdot \left(uT\right)-\nabla \cdot \left(\frac{\lambda }{\rho C_p}\nabla T\right)=-\frac{\stackrel{\dot{}}{m}L}{\rho C_p}$ (10)

式中，T为温度，K；L为单位质量蒸汽冷凝所释放的相变潜热，J/kg； $\stackrel{\dot{}}{m}$ 为单位体积蒸发率，kg/(m³·s)。

2.2. 相变模型

求解器中编译有两种相变模型——Lee模型和KineticGasEvaporation model模型，可针对每一类型的相变过程引入特定相变模型。本次研究涉及金属凝固过程，选取Lee相变驱动模型。

1) 由Lee模型驱动的相间传质方程

当C > 0时为熔化/蒸发过程：

$\stackrel{\dot{}}{m}=C\rho \alpha \frac{T-T_{\text{activate}}}{T_{\text{activate}}}\text{for}T>T_{\text{activate}}$ (11)

当C < 0时为凝固/冷凝过程：

$\stackrel{\dot{}}{m}=-C\rho \alpha \frac{T-T_{\text{activate}}}{T_{\text{activate}}}\text{for}T<T_{\text{activate}}$ (12)

式中，C为控制相变强度的模型系数；T_activate为相变温度，K。

2) 由KineticGasEvaporation模型驱动的相间传质方程

KineticGasEvaporation模型基于动力学气体理论，结合Hertz-Knudsen公式 ‎[8] ，给出了平坦界面的蒸发–冷凝质量通量：

$F=C\left(p-p_{\text{sat}}\right)\cdot \sqrt{\frac{M}{2\pi R{T}_{\text{activate}}}}$ (13)

Clapeyron-Clausius方程 ‎[9] 将压力与饱和条件下的温度联系在一起：

$\frac{\text{d}p}{\text{d}T}=-\frac{L}{T\left({\nu }_{\text{vap}}-{\nu }_{\text{liq}}\right)}$ (14)

假设液体和蒸汽处于平衡状态，即界面的调节系数是等效的。结合上两式得到如下相变方程：

$F=2\frac{C}{2-C}\sqrt{\frac{M}{2\pi R{T}_{\text{activate}}}}L\frac{{\rho }_{\text{l}}{}_{\text{iq}}{\rho }_{\text{v}}{}_{\text{ap}}}{\left({\rho }_{\text{l}}{}_{\text{iq}}-{\rho }_{\text{vap}}\right)}\frac{\left(T-T_{\text{activate}}\right)}{T_{\text{activate}}}$ (15)

式中，F为单位体积相变率，kg/(m³·s)；M为分子量；R为气体常数，J/(mol·K)；p为蒸汽，Pa；p_sat为饱和压力，Pa；n为运动粘度，m²/s；当C > 0时为蒸发过程，当C < 0时为冷凝过程。

2.3. 湍流模型

在icoReactingMultiphaseInterFoam的turbulenceProperties中，可以设置用于求解的流动类型。这里有三种不同的模型以供选择，包括laminar、LES、RAS。LES和RAS是由多个子类模型组成的。其中，LES分为“等容的LES湍流模型”和“各向异性的LES湍流模型” ‎[10] 。RAS分为两类，一类是不可压缩的，一类是可压缩的。本次研究指定RAS中的k-e模型，流动方程见下式。

1) 湍流动能方程

$\frac{D}{Dt}\left(\rho \kappa \right)=\nabla \cdot \left(\rho D_{\kappa }\nabla \kappa \right)+P_{\kappa }-\rho \epsilon$ (16)

式中，k为湍流动能，m²/s²；D_k为湍流动能的有效扩散系数；P_k为湍流动能产生率，m²/s³；e为湍流动能耗散率，m²/s³。

2) 湍流动能的耗散率方程

$\frac{D}{Dt}\left(\rho \epsilon \right)=\nabla \cdot \left(\rho D_{\epsilon }\nabla \epsilon \right)+\frac{C_1\epsilon }{\kappa }\left(P_{\kappa }+C_3\frac{2}{3}\kappa \nabla \cdot u\right)-C_2\rho \frac{{\epsilon }^2}{\kappa }$ (17)

式中，D_e为湍流动能耗散率的有效扩散系数；C₁、C₂和C₃是模型系数。

3) 湍流粘度方程

$v=C_{\nu }\frac{{\kappa }^2}{\epsilon }$ (18)

式中，n为湍流粘度，m²/s；C_n为湍流粘度的模型系数。

4) OpenFOAM实现

OpenFOAM软件对上述方程做了进一步变化，以实现湍流模型在各类求解器中的调用；修正后的湍流动能方程和湍流动能的耗散率方程如下：

$\frac{\partial }{\partial t}\left(\alpha \rho \kappa \right)+\nabla \cdot \left(\alpha \rho u\kappa \right)-{\nabla }^2\left(\alpha \rho D_{\kappa }\kappa \right)=\alpha \rho G-\left(\frac{2}{3}\alpha \rho \nabla \cdot u\kappa \right)-\left(\alpha \rho \frac{\epsilon }{\kappa }\kappa \right)+S_{\kappa }+S_{\text{fvOptions}}$ (19)

$\frac{\partial }{\partial t}\left(\alpha \rho \epsilon \right)+\nabla \cdot \left(\alpha \rho u\epsilon \right)-{\nabla }^2\left(\alpha \rho D_{\epsilon }\epsilon \right)=C_1\alpha \rho G\frac{\epsilon }{\kappa }-\left(\left(\frac{2}{3}{C}_1-C_3\right)\alpha \rho \nabla \cdot u\epsilon \right)-\left(C_2\alpha \rho \frac{\epsilon }{\kappa }\epsilon \right)+S_{\epsilon }+S_{\text{fvOptions}}$ (20)

式中，G是由于雷诺应力张量的各向异性部分而产生的湍流动能产生率，m²/s³；在turbulenceProperties字典中，模型系数值分别为C_n = 0.09、C₁ = 1.44、C₂ = 1.92、C₃ = 0.0、s_k = 1.0、s_e = 1.3。

3. 建模

3.1. 数值计算模型

Figure 1. Numerical model for hydraulic fragmentation of molten metal

图1. 熔融金属水力碎化数值计算模型

为研究熔融锡液柱水力碎化过程的界面行为及相关物理现象，做出以下假设：熔融锡液柱周围冷却剂不发生相变，引入熔融锡液凝固相变模型，即忽略多相流中水的蒸发冷凝效应，只考虑水力碎化作用。高温熔融锡液柱水力碎化模拟的物理模型及网格划分如图1所示。建立二维矩形计算域，计算域的长、宽分别为700 mm ´ 140 mm，进口处为20 mm。计算模型的两侧及底部黑色边界类型定义为无滑移、零速度梯度的壁面(wall)，用于约束流体及固体域；顶部两端的灰色边界定义为压力出口(pressure-outlet)，当迭代期间出现回流时会收敛速度更好；顶部的红色边界定义为熔融锡液柱的速度入口(velocity-inlet)，锡液作为不可压缩流体，速度入口可进一步控制流量。整个模拟不考虑空气的夹带效应，流动模型为层流。熔融锡液柱经顶部进口以一定速度流入水域内，如图1所示红色域为高温熔融锡液，蓝色域为水；在网格划分时，对锡液柱及液柱附近的网格进行局部加密。在常温常压环境下，锡液初始温度为1050℃、初始速度为3 m/s；计算域内水的初始温度为27℃。边界条件根据模拟情况设置如下：进口处压力为calculated边界条件，压力由其他边界条件及参数来定义；温度与速度为固定值，因此均为fixedValue边界条件；相体积分数为inletOutlet边界条件，该边界条件由mixed边界条件衍生而来，当流体流出时指定为zeroGradient条件，当流体流入时则是fixedValue条件。出口处速度为pressureInletOutletVelocity边界，流体流出时指定zeroGradient条件，流入时则根据计算块(patch)法线方向上的通量分配速度，边界条件详细信息见表1。

Table 1. Set up a series of boundary condition

表1. 边界条件设置

3.2. 网格相关性验证

为准确捕获界面，同时节约计算时间成本，设置三种网格以进行网格相关性研究，这三种四面体网格质量均为1。监测点设置在水平中心0.07 m、纵向高度0.6 m处；如图2所示，比较监测点处温度随时间的变化曲线，确定后续计算分析所用网格单元格总数为276,202个。

Figure 2. Variation curves of temperature at measuring points according to different grid numbers

图2. 不同网格数下监测点处温度变化曲线

4. 仿真结果与分析

4.1. 水力碎化过程分析

Figure 3. Simulation diagrams of molten tin fragmentation at diameter is 20 mm and velocity is 3 m/s

图3. 液柱直径为20 mm、初速度为3 m/s熔融锡液柱碎化相图

如图3所示为液柱直径为20 mm、初速度为3 m/s熔融金属锡液碎化相图。在0.015 s时，由于顶部相对高压区域及与金属锡和水密度差相关的R-T不稳定性的综合作用，金属锡液柱前端发生形变、横向扩张形成伞状；伞状顶部两侧金属细丝向外翻转，受边界层剪切效应发生断裂。由于熔融锡液与水刚开始直接接触，锡液总量小，此时碎化效应并不明显。从0.015 s开始，金属锡液柱持续以3 m/s的速度流入水中，液柱前端横向扩张的伞状现象因熔融金属锡液柱入水及液柱碎化造成的扰动不再显著，受与速度差相关的K-H不稳定波动影响而造成的碎化现象更加突出。在0.015 s~0.039 s内，熔融锡液柱侧面出现明显的水力不稳定波动，随着该不稳定波持续生长，从0.03 s开始可以明显地观察到液柱两侧出现多个不对称的倒钩状金属细丝。金属细丝在0.033 s时发生明显的断裂现象，剥离液柱表面。

Figure 4. Disturbance near the wire on the metal side

图4. 液柱侧面金属细丝附近的扰流

图4为0.039 s时熔融锡液附近的扰流速度矢量图，因K-H不稳定性波动形成的倒钩状金属细丝附近往往伴随着流体涡旋流动。该扰流对金属碎化的影响表现在扰流能促使液柱表面金属细丝拉长而加剧细丝断裂、碎化。断裂首先发生在倒钩状细丝尖端，该部分剥离液柱主体的同时，由于熔融锡液内部粘性力的存在，剩余部分将继续发展，再次经历拉长、断裂、碎化的过程，最终使液柱发生整体断裂。

熔融锡液遇水后的碎化过程为高密度比(液柱密度与水的密度之比)模拟，对比如图5所示陈竞覃等 ‎[11] 高密度比下水银射流碎化过程拍摄图、原理示意图以及如图6所示Iwasawa等 ‎[12] 低熔点合金液柱射流碎化实验拍摄图。由于模拟不考虑空气夹带和蒸发相变，观察图5中液柱入水后侧面的不稳定性波动情况及图6中液柱前端的横向扩张和液柱侧面的剥离现象，与数值模拟结果基本一致。

Figure 5. Photo and schematic diagram of jet breakup behavior

图5. 射流碎化过程拍摄图及原理示意图

Figure 6. Photo of jet breakup behavior

图6.射流碎化实验拍摄图

Figure 7. Simulation diagrams of molten tin’s solidification at diameter is 20 mm and velocity is 3 m/s

图7. 液柱直径为20 mm、初速度为3 m/s时熔融锡液柱凝固相图

如图7所示，熔融锡液柱与水直接接触，液柱顶部及两侧率先发生凝固。0.039 s时对比凝固相图与熔融锡液柱碎化相图，明显发现脱离液柱的锡液滴在失去液柱内部锡液的热量供给后，凝固速度加快，凝固相体积分数较大。图8为固态锡体积分数随时间的变化曲线，固态锡生成速率随反应的深入而提升，这是由于锡液碎化剧烈，形成众多小尺寸金属颗粒，这些颗粒储热少，极容易被周围冷却水冷却凝固。

Figure 8. Curve for theaveragealphas of solid tin per unit time

图8. 固态锡体积分数均值随时间的变化曲线

4.2. 液柱直径对碎化过程的影响

Figure 9. Fragmentation processes of molten tin at diameter is 10 mm

图9. 液柱直径为10 mm时熔融锡液柱碎化过程

Figure 10. Fragmentation processes of molten tin at diameter is 30 mm

图10. 液柱直径为30 mm时熔融锡液柱碎化过程

为研究金属液柱直径对水力碎化的影响，除了进口尺寸20 mm外，另设置了10 mm和30 mm两种进口尺寸，分别进行同一数值工况的水力碎化模拟研究，如图9、图10所示。R-T不稳定性波动与相对速度无关故而保持不变，因此相图重点关注由K-H不稳定性波动造成的液柱侧面的水力碎化现象。当熔融锡液直径为10 mm时，0.0135 s时液柱已经下落至水深0.19 m处；而熔融液柱直径为20 mm时，下落至水深0.19 m处需要耗时0.039 s。熔融锡液柱在水内下落过程中，由于液柱直径较小，内部粘性力弱，液柱侧面受到K-H不稳定性波动的影响，多处出现断裂现象并发生剧烈的碎化。当熔融锡液直径为30 mm时，由于进口处水压的存在，直径为30 mm的液柱下落高度进口处所受阻力较大，熔融锡液柱下落高度耗时0.056 s仅下落了0.13 m，低于相同时间内直径为10 mm和20 mm时熔融锡液柱下落高度。反应持续至0.054 s时，熔融锡液在水下形成直径较为稳定的液柱，液柱两侧受K-H不稳定性波动的影响，两侧出现明显的弯钩状金属细丝并断裂，液柱碎化现象较之前显著提升。

因此，液柱入水初速度相同时，随着液柱直径的增大，液柱入水用于克服入水阻力的动能消耗越大，在相同时间内液柱入水深度越小，如图11所示。同时，相同入水时间内，液柱化程度与液柱的直径成反比，液柱直径越小，熔融锡液柱碎化得越剧烈。

Figure 11. Curve for height of the top of liquid column entering water with time

图11. 液柱顶部入水高度随时间的变化曲线

4.3. 液柱入水初速度对碎化过程的影响

Figure 12. Fragmentation processes of molten tin at velocity is 5 m/s

图12. 初速度为5 m/s时熔融锡液柱碎化过程

Figure 13. Fragmentation processes of molten tin at velocity is 10 m/s

图13. 初速度为10 m/s时熔融锡液柱碎化过程

为研究金属液柱下落速度对水力碎化的影响，除了初速度3 m/s外，另设置了5 m/s和10 m/s两种进口速度，分别进行同一数值工况的水力碎化模拟研究，如图12、图13所示。对比分析初速度为3 m/s、5 m/s、10 m/s时熔融金属液柱碎化过程，当熔融锡液柱与水直接接触反应至0.02 s时，熔融锡液柱入水初速度增大意味着熔融液柱入水时的动能越大，液柱克服顶部高压及水的阻力作用向下流动过程中，初速度大的液柱顶部碎化更剧烈。同时，由于单位时间内流入水内的锡液总量增加，在液柱直径相同的情况下，热量传递速率加快，水力扰动变大，加速了熔融锡液柱的碎化。图14为固体锡初速度10 m/s时，同一时间段内固态锡生成量高于其他两种初速度下的固态锡生成量，即初速度为10 m/s的锡凝固速度较快。

因此，熔融锡液柱入水后的碎化程度与液柱入水时的初速度成正比，液柱入水时初速度越大，水力波动越剧烈，熔融锡液柱顶端及侧面的碎化程度得越高，锡凝固速度越快。

Figure 14. Curves for the averagealphas of solid tin per unit time at different velocities

图14. 不同入水速度下固态锡体积分数均值随时间的变化曲线

5. 结论

熔融锡液水力碎化模拟研究得到以下结论：

1) 水力碎化现象方面：熔融锡液柱碎化是三种水力因素共同作用的结果。液柱顶部受与密度相关的瑞利–泰勒(Rayleigh-Taylor, R-T)不稳定性影响呈伞状向两端横向发展，液柱侧面受与相对速度相关的开尔文–亥姆霍兹(Kelvin-Helmholtz, K-H)不稳定性影响生成多个不对称的倒钩状金属细丝，受界面剪切力的影响液柱伞状顶部横向发展的两个尖端及液柱侧面倒钩状金属细丝被拉长并断裂。

2) 水力碎化剧烈程度与液柱直径、液柱入水初速度的敏感性分析：液柱入水初速度相同时，随着熔融锡液柱直径的增大，相同时间内液柱入水深度越小；液柱碎化剧烈程度与液柱直径成反比。液柱入水直径相同时，熔融锡液柱入水后的碎化程度、锡凝固速度与液柱入水时初速度成正比。

基金项目

国家重点研发计划项目(2017YFC0805101)；江苏省绿色过程装备重点实验室开放课题基金资助项目(GPE201901)；2022年江苏省研究生科研与实践创新计划项目(SJCX22_1434)。

NOTES

^*第一作者。