Numerical Simulations of Molten Breakup Behaviors of a de Laval-Type Nozzle, and the E ﬀ ects of Atomization Parameters on Particle Size Distribution

: In this work, an atomizer with a de Laval-type nozzle is designed and studied by commercial computational ﬂuid dynamics (CFD) software, and the secondary breakup process during atomization is simulated by two-way coupling and the discrete particle model (DPM) using the Euler-Lagrange method. The simulation result demonstrates that the gas ﬂow patterns greatly change with the introduction of liquid droplets, which clearly indicates that the mass loading e ﬀ ect is quite signiﬁcant as a result of the gas-droplet interactions. An hourglass shape of the cloud of disintegrating molten metal particles is observed by using a stochastic tracking model. Finally, this simulation approach is used for the quantitative evaluation of the e ﬀ ects of altering the atomizing process conditions (gas-to-melt ratio, operating pressure P , and operating gas temperature T ) and nozzle geometry (protrusion length h , half-taper angle α , and gas slit nozzle diameter D ) on the particle size distribution of the powders produced.


Introduction
High quality metal powder is the basis for applications such as powder metallurgy, metal injection molding, spray coating, laser cladding, metal welding, and additive manufacturing [1][2][3]. Rising surface and molding technologies demand finer metal powders and a narrow particle size distribution, which has greatly promoted the development of powder production technology. Due to its versatility, quality, and purity, gas atomization is currently considered the best way to produce high quality, spherical metal powders. The basic principle of gas atomization is utilizing the kinetic energy of a high-speed gas jet to breakup liquid molten metal, and the broken droplets are then spheroidized under the surface tension and undergo rapid solidification to obtain spherical metal powders. The gas atomization breakup process involves the complex multiphase coupling between the high-speed gas jet and the high-temperature molten metal stream, which occurs under the conditions of an extremely high cooling rate, momentum exchange, and heat exchange [4]. Although gas atomization has been where Ma is the gas jet Mach number, the specific heat ratio k for the nitrogen is k = 1.4, and A/A 0 varies with the Mach number, as shown in Figure 1. It can be observed from Figure 1 that when the area ratio A/A 0 is determined, the Mach number follows. Additionally, the same area ratio A/A 0 corresponds to two Mach numbers, which respectively correspond to the subsonic flow (Ma < 1) and the supersonic flow (Ma > 1) states in the de Laval tube.
Processes 2020, 8 where Ma is the gas jet Mach number, the specific heat ratio k for the nitrogen is k = 1.4, and A/A0 varies with the Mach number, as shown in Figure 1. It can be observed from Figure 1 that when the area ratio A/A0 is determined, the Mach number follows. Additionally, the same area ratio A/A0 corresponds to two Mach numbers, which respectively correspond to the subsonic flow (Ma < 1) and the supersonic flow (Ma > 1) states in the de Laval tube. The convergent wall curve in the de Laval nozzle is generally determined by the Witoszynski Curve Method, the Batchelor-Shaw Curve Method, or the High-order Curve Method [21][22][23][24]. The divergent wall curve is mainly determined by the Characteristic Method or the Foelsch Method [25].
In this paper, as shown in Figure 2, the convergent section curve AB is determined by the Highorder Curve Method, and two methods are used for the diverging section (BC and CD), which includes the Characteristic Method. The specific curve equations are as follows.  The convergent wall curve in the de Laval nozzle is generally determined by the Witoszynski Curve Method, the Batchelor-Shaw Curve Method, or the High-order Curve Method [21][22][23][24]. The divergent wall curve is mainly determined by the Characteristic Method or the Foelsch Method [25].
In this paper, as shown in Figure 2, the convergent section curve AB is determined by the High-order Curve Method, and two methods are used for the diverging section (BC and CD), which includes the Characteristic Method. The specific curve equations are as follows. where Ma is the gas jet Mach number, the specific heat ratio k for the nitrogen is k = 1.4, and A/A0 varies with the Mach number, as shown in Figure 1. It can be observed from Figure 1 that when the area ratio A/A0 is determined, the Mach number follows. Additionally, the same area ratio A/A0 corresponds to two Mach numbers, which respectively correspond to the subsonic flow (Ma < 1) and the supersonic flow (Ma > 1) states in the de Laval tube. The convergent wall curve in the de Laval nozzle is generally determined by the Witoszynski Curve Method, the Batchelor-Shaw Curve Method, or the High-order Curve Method [21][22][23][24]. The divergent wall curve is mainly determined by the Characteristic Method or the Foelsch Method [25].
In this paper, as shown in Figure 2, the convergent section curve AB is determined by the Highorder Curve Method, and two methods are used for the diverging section (BC and CD), which includes the Characteristic Method. The specific curve equations are as follows.  Convergent section curve AB: where R 0 is the de Laval throat radius, R 1 is the gas inlet radius, L 1 is the total length of the convergent section, x is the coordinate of the x-axis in Figure 2, and R is the transition arc radius corresponding to the x position. Divergent section curve BC and CD: Using the Characteristic Method to obtain the line equation and connect it smoothly.
where L 2 is the total length of the divergent section, θ is the divergent angle, r 1 and r 2 are the radius of the analytical circle, and R 2 is the gas outlet radius. In this paper, in order to obtain a higher gas velocity at the lower atomizing pressure, the designed Ma = 4.04, and the corresponding A 2 /A 0 = 11.56, R 0 = 0.3 mm, R 1 = 3.2 mm, R 2 = 1.02 mm, r 1 = 15.0 mm, r 2 = 14. 3 mm, L 1 = 8.0 mm, L 2 = 6.67 mm, and θ = 25 • .
From the above content, we confirmed the parameters of the de Laval-type gas slit. Typically, the structure of the atomizing nozzle depends on the half-taper angle α, the gas injection ring diameter D, and the protrusion length of the melt delivery tube h. These critical parameters will greatly affect the gas flow field structure and the metal melt breakup process during atomization [26].

Modelling the Fluids
The gas flow field in atomization process can be described by k-ε, k- Convergent section curve AB: where R0 is the de Laval throat radius, R1 is the gas inlet radius, L1 is the total length of the convergent section, x is the coordinate of the x-axis in Figure 2, and R is the transition arc radius corresponding to the x position. Divergent section curve BC and CD: Using the Characteristic Method to obtain the line equation and connect it smoothly.

:
, 0 : , where L2 is the total length of the divergent section, θ is the divergent angle, r1 and r2 are the radius of the analytical circle, and R2 is the gas outlet radius. In this paper, in order to obtain a higher gas velocity at the lower atomizing pressure, the designed Ma = 4.04, and the corresponding A2/A0 = 11.56, R0 = 0.3 mm, R1 = 3.2 mm, R2 = 1.02 mm, r1 = 15.0 mm, r2 = 14. 3 mm, L1 = 8.0 mm, L2 = 6.67 mm, and θ = 25°. From the above content, we confirmed the parameters of the de Laval-type gas slit. Typically, the structure of the atomizing nozzle depends on the half-taper angle α, the gas injection ring diameter D, and the protrusion length of the melt delivery tube h. These critical parameters will greatly affect the gas flow field structure and the metal melt breakup process during atomization [26].

Modelling the Fluids
The gas flow field in atomization process can be described by k-ε, ɷ and the Reynolds Averaging Naiver-Stokes (RANS). The turbulent gas flow behavior can be properly described by the k-ε model with less simulation, so we chose it as the turbulence model of the gas flow field [27].
To simulate the secondary breakup process of the droplets, the discrete particle model (DPM) is implemented in ANSYS Fluent, and the Euler-Lagrange approach is utilized to describe the melt breakup behavior. We assume that the metal melt stream becomes an ensemble of large droplets after the primary breakup, and further transforms into fine droplets after the secondary breakup process. The motion and trajectories of the particle droplets can be predicted by force balance, which equates the particle's inertia with the forces that act upon it. The equation is written in the Lagrange reference frame [28]: Here, , / is the drag force per unit particle mass, and the Fi term includes the incorporation of additional forces.
is the droplet relaxation time, defined as The droplet relative Reynolds number Re is defined as , In Equations (4)-(6), the index i represents the three coordinates of the particle position and velocity components, ui and up,i refers to the particle velocity, respectively; μ is the molecular viscosity of the gas, ρg and ρp represents the density of the gas and particle, respectively; Fi refers to the , and the Reynolds Averaging Naiver-Stokes (RANS). The turbulent gas flow behavior can be properly described by the k-ε model with less simulation, so we chose it as the turbulence model of the gas flow field [27].
To simulate the secondary breakup process of the droplets, the discrete particle model (DPM) is implemented in ANSYS Fluent, and the Euler-Lagrange approach is utilized to describe the melt breakup behavior. We assume that the metal melt stream becomes an ensemble of large droplets after the primary breakup, and further transforms into fine droplets after the secondary breakup process. The motion and trajectories of the particle droplets can be predicted by force balance, which equates the particle's inertia with the forces that act upon it. The equation is written in the Lagrange reference frame [28]: Here, u i − u p,i /t r is the drag force per unit particle mass, and the F i term includes the incorporation of additional forces. t r is the droplet relaxation time, defined as The droplet relative Reynolds number Re is defined as Processes 2020, 8, 1027

of 18
In Equations (4)-(6), the index i represents the three coordinates of the particle position and velocity components, u i and u p,i refers to the particle velocity, respectively; µ is the molecular viscosity of the gas, ρ g and ρ p represents the density of the gas and particle, respectively; F i refers to the incorporation of additional acceleration, d p is the particle diameter, and the C d is the drag coefficient [17,28].
The process of secondary breakup can be described by four breakup models: the TAB model, the WAVE model, the Stochastic Secondary Droplet (SSD) model, and the KHRT model. The TAB model is based on Taylor's analogy [29] between an oscillating and distorting droplet and a spring mass system. The WAVE model considers the breakup of the droplets to be induced by the relative velocity between the gas and liquid phases. It assumes that the time of breakup and the resulting droplet size are related to the fastest-growing Kelvin-Helmholtz (KH) instability [30]. The SSD model [31] treats breakup as a discrete random event that results in a distribution of diameter scales over a range. The KHRT model combines the effects of the Kelvin-Helmholtz model with the Rayleigh-Taylor (RT) instabilities. It assumes that there is a liquid core near the nozzle area, and the liquid core droplets are suddenly accelerated when they are ejected into the freestream, and RT instability becomes the main influence. From the Levich theory [32], we can obtain the length of the liquid core where C l is the Levich constant and d 0 is a reference nozzle diameter specified for atomizer injections in Fluent.
The KH model is based on linear stability analysis of a liquid jet, with its fastest growing mode contributing to the eventual breakup and generation of new droplets. The wave length of the fastest growing mode and growth rate was found numerically to be [33] λ KH = 9.02R 1 + 0.45Z 0.5 1 + 0.47T 0.7 1 + 0.87We 1.67 g 0.6 (8) where λ KH is the wave length of the KH wave, ω KH is its growth rate, R is the jet radius, σ is the surface tension of the melt, Z is the Ohnesorge number, defined as and T is the Taylor number, defined as The relationship between the radius of the new droplet r generated in the primary breakup process and the wave length λ KH is described below where B KH is the size constant of the KH breakup model. The KH breakup time scale t KH is computed by where C KH is the time constant of the KH breakup, which is a value between 1.732 and 20 [28]. Like the Kelvin-Helmholtz model, the RT model is founded on wave instabilities on the surface of the droplet. The frequency of the fastest growing wave is calculated as where g t is the droplet acceleration in the direction of motion. The corresponding wave number is given by Breakup occurs after RT waves have been growing for a longer time than the breakup time t RT , defined as where C t is the Rayleigh-Taylor breakup time constant, which has a default value of 0.5 [28]. When the tracked wave length 2πC RT /K RT is smaller than the local droplet diameter. The radius of the smaller child droplets is calculated as follows where C RT is the breakup radius constant, which has a default value of 0.1 [28]. The KHRT effects are calculated and considered for breakup outside the liquid core. Within the liquid core, only aerodynamic breakup computed by the WAVE breakup model is considered. Normally, the RT instability grows faster when droplet acceleration is high, and this effect is dominant for a high Weber number stream [28].

Model Implementation
Due to the symmetry of the atomizer model, a 2D axis-symmetric computational domain is performed on one half of the axial section of the atomizer chamber with a 200 mm length and 80 mm width, and the boundary conditions are illustrated in Figure 3. where gt is the droplet acceleration in the direction of motion. The corresponding wave number is given by Breakup occurs after RT waves have been growing for a longer time than the breakup time , defined as = (16) where is the Rayleigh-Taylor breakup time constant, which has a default value of 0.5 [28]. When the tracked wave length 2π / is smaller than the local droplet diameter. The radius of the smaller child droplets is calculated as follows = (17) where is the breakup radius constant, which has a default value of 0.1 [28]. The KHRT effects are calculated and considered for breakup outside the liquid core. Within the liquid core, only aerodynamic breakup computed by the WAVE breakup model is considered. Normally, the RT instability grows faster when droplet acceleration is high, and this effect is dominant for a high Weber number stream [28].

Model Implementation
Due to the symmetry of the atomizer model, a 2D axis-symmetric computational domain is performed on one half of the axial section of the atomizer chamber with a 200 mm length and 80 mm width, and the boundary conditions are illustrated in Figure 3. To evaluate the droplet particle size distribution produced by the secondary breakup, a transient two-way turbulence coupling Euler-Lagrange method was utilized. The dispersion of the particles can be expressed by the stochastic tracking model with the discrete random walk model. To ensure solution convergence and capture the molten melt breakup process, very small time steps of 1 × 10 −7 s were employed in this simulation. Nitrogen was selected as a compressible atomizing gas. The density was set to "ideal gas", and the viscosity was calculated by the Sutherland law. Other properties were taken from the Fluent database. A Fe-based self-fluxing alloy (FeNiCrBSiNbC) was selected as a molten metal material, where the obtained powder was used for laser cladding. To To evaluate the droplet particle size distribution produced by the secondary breakup, a transient two-way turbulence coupling Euler-Lagrange method was utilized. The dispersion of the particles can be expressed by the stochastic tracking model with the discrete random walk model. To ensure solution convergence and capture the molten melt breakup process, very small time steps of 1 × 10 −7 s were employed in this simulation. Nitrogen was selected as a compressible atomizing gas. The density was set to "ideal gas", and the viscosity was calculated by the Sutherland law. Other properties were taken from the Fluent database. A Fe-based self-fluxing alloy (FeNiCrBSiNbC) was selected as a molten metal material, where the obtained powder was used for laser cladding. To approximate the manufacturing process of the metal powders, this simulation utilized dynamic thermo-physical properties, as shown in Figure 4. The initial diameters of the inner parent droplets were set to 0.4~4.0 mm with a Rosin-Rammler distribution, and the droplets were injected at the top corner of the melt nozzle at a temperature of 1890 K (200 K superheat) with a mass flow rate of 0.1 kg/s, corresponding to a gas-to-melt ratio (GMR) of 3.09. The de Laval-type nozzle operated at a pressure of 2.0 MPa and temperature of 300 K during the entire simulation process unless otherwise noted. approximate the manufacturing process of the metal powders, this simulation utilized dynamic thermo-physical properties, as shown in Figure 4. The initial diameters of the inner parent droplets were set to 0.4~4.0 mm with a Rosin-Rammler distribution, and the droplets were injected at the top corner of the melt nozzle at a temperature of 1890 K (200 K superheat) with a mass flow rate of 0.1 kg/s, corresponding to a gas-to-melt ratio (GMR) of 3.09. The de Laval-type nozzle operated at a pressure of 2.0 MPa and temperature of 300 K during the entire simulation process unless otherwise noted. Considering the potential influence of mesh size, we assumed a maximum generic mesh spacing of δ. Mesh spacings in the range of 1/2δ, 3/5δ, 2/3δ, and 4/5δ were produced to verify the independence of the mesh space in the simulation. The d50 value after a sufficient number of iteration steps was used to evaluate the influence of the mesh space, and it was observed that the d50 difference between the coarse mesh spacing of 4/5δ and the finest mesh spacing of 1/2δ was well below the 2% error, as can be seen in Figure 5. Therefore, the structured grid has a cell number of 155,048, with the mesh space of 3/5δ selected due to the considerations of the simulation amount and accuracy. Figure  6 plots the influence of simulation time on the average diameter of the droplets in different mass flow rates. It was determined that a simulation time of approximately 5.0 ms was required to reach a stable d50 value; beyond that, there is no substantial change in the particle size distribution, indicating that the system reached a steady-state condition.  Considering the potential influence of mesh size, we assumed a maximum generic mesh spacing of δ. Mesh spacings in the range of 1/2δ, 3/5δ, 2/3δ, and 4/5δ were produced to verify the independence of the mesh space in the simulation. The d 50 value after a sufficient number of iteration steps was used to evaluate the influence of the mesh space, and it was observed that the d 50 difference between the coarse mesh spacing of 4/5δ and the finest mesh spacing of 1/2δ was well below the 2% error, as can be seen in Figure 5. Therefore, the structured grid has a cell number of 155,048, with the mesh space of 3/5δ selected due to the considerations of the simulation amount and accuracy. Figure 6 plots the influence of simulation time on the average diameter of the droplets in different mass flow rates. It was determined that a simulation time of approximately 5.0 ms was required to reach a stable d 50 value; beyond that, there is no substantial change in the particle size distribution, indicating that the system reached a steady-state condition. approximate the manufacturing process of the metal powders, this simulation utilized dynamic thermo-physical properties, as shown in Figure 4. The initial diameters of the inner parent droplets were set to 0.4~4.0 mm with a Rosin-Rammler distribution, and the droplets were injected at the top corner of the melt nozzle at a temperature of 1890 K (200 K superheat) with a mass flow rate of 0.1 kg/s, corresponding to a gas-to-melt ratio (GMR) of 3.09. The de Laval-type nozzle operated at a pressure of 2.0 MPa and temperature of 300 K during the entire simulation process unless otherwise noted. Considering the potential influence of mesh size, we assumed a maximum generic mesh spacing of δ. Mesh spacings in the range of 1/2δ, 3/5δ, 2/3δ, and 4/5δ were produced to verify the independence of the mesh space in the simulation. The d50 value after a sufficient number of iteration steps was used to evaluate the influence of the mesh space, and it was observed that the d50 difference between the coarse mesh spacing of 4/5δ and the finest mesh spacing of 1/2δ was well below the 2% error, as can be seen in Figure 5. Therefore, the structured grid has a cell number of 155,048, with the mesh space of 3/5δ selected due to the considerations of the simulation amount and accuracy. Figure  6 plots the influence of simulation time on the average diameter of the droplets in different mass flow rates. It was determined that a simulation time of approximately 5.0 ms was required to reach a stable d50 value; beyond that, there is no substantial change in the particle size distribution, indicating that the system reached a steady-state condition.

Flow Field Analysis
The de Laval nozzle will show different working modes under different pressures. Therefore, it will affect the particle size distribution of the droplets after breakup. The structure of the flow field was first studied before the melt was injected. Figure 7 presents the velocity contours developed in the region adjacent to the nozzle exit of the atomizing jet. The typical flow field of a nozzle slit is an under-expanded jet flow, a clear Mach disk, and a closed-wake structure [14]. Unlike the typical nozzle slit, the atomizing gas can be fully expanded in the de Laval nozzle due to its special convergence-divergence structure. The gas is constantly accelerating and reaches the speed of sound in the throat area, and then reaches supersonic speed in the divergence area. Therefore, parallel supersonic gas jets can be obtained in the nozzle exit. The maximum gas velocity that can be reached is 623 m/s, and the supersonic gas jets then meet and reflect below the delivery tube to form a series of periodic diamond-shaped reflected waves.  Figure 8 shows the influence of the particle droplets on the gas flow field. From the contours, we can see that the gas velocity on the central axis is significantly decreased due to the breakup effect and the mass loading effect. Additionally, the kinetic energy of the gas jets has dropped dramatically, and the maximum velocity on the central axis is reduced by almost 40%. When comparing the gas flow field with the two-way coupling, significant changes can be observed in the behavior of the atomizing jets, especially under the recirculation zone; the droplets converged in this area and destroyed the stability of the velocity field, which resulted in the velocity become disrupted and the disappearance of the periodic diamond-shock reflected waves. These findings show the importance

Flow Field Analysis
The de Laval nozzle will show different working modes under different pressures. Therefore, it will affect the particle size distribution of the droplets after breakup. The structure of the flow field was first studied before the melt was injected. Figure 7 presents the velocity contours developed in the region adjacent to the nozzle exit of the atomizing jet. The typical flow field of a nozzle slit is an under-expanded jet flow, a clear Mach disk, and a closed-wake structure [14]. Unlike the typical nozzle slit, the atomizing gas can be fully expanded in the de Laval nozzle due to its special convergence-divergence structure. The gas is constantly accelerating and reaches the speed of sound in the throat area, and then reaches supersonic speed in the divergence area. Therefore, parallel supersonic gas jets can be obtained in the nozzle exit. The maximum gas velocity that can be reached is 623 m/s, and the supersonic gas jets then meet and reflect below the delivery tube to form a series of periodic diamond-shaped reflected waves.

Flow Field Analysis
The de Laval nozzle will show different working modes under different pressures. Therefore, it will affect the particle size distribution of the droplets after breakup. The structure of the flow field was first studied before the melt was injected. Figure 7 presents the velocity contours developed in the region adjacent to the nozzle exit of the atomizing jet. The typical flow field of a nozzle slit is an under-expanded jet flow, a clear Mach disk, and a closed-wake structure [14]. Unlike the typical nozzle slit, the atomizing gas can be fully expanded in the de Laval nozzle due to its special convergence-divergence structure. The gas is constantly accelerating and reaches the speed of sound in the throat area, and then reaches supersonic speed in the divergence area. Therefore, parallel supersonic gas jets can be obtained in the nozzle exit. The maximum gas velocity that can be reached is 623 m/s, and the supersonic gas jets then meet and reflect below the delivery tube to form a series of periodic diamond-shaped reflected waves.  Figure 8 shows the influence of the particle droplets on the gas flow field. From the contours, we can see that the gas velocity on the central axis is significantly decreased due to the breakup effect and the mass loading effect. Additionally, the kinetic energy of the gas jets has dropped dramatically, and the maximum velocity on the central axis is reduced by almost 40%. When comparing the gas flow field with the two-way coupling, significant changes can be observed in the behavior of the atomizing jets, especially under the recirculation zone; the droplets converged in this area and destroyed the stability of the velocity field, which resulted in the velocity become disrupted and the disappearance of the periodic diamond-shock reflected waves. These findings show the importance  Figure 8 shows the influence of the particle droplets on the gas flow field. From the contours, we can see that the gas velocity on the central axis is significantly decreased due to the breakup effect and the mass loading effect. Additionally, the kinetic energy of the gas jets has dropped dramatically, and the maximum velocity on the central axis is reduced by almost 40%. When comparing the gas flow field with the two-way coupling, significant changes can be observed in the behavior of the atomizing jets, especially under the recirculation zone; the droplets converged in this area and destroyed the stability of the velocity field, which resulted in the velocity become disrupted and the disappearance of the periodic diamond-shock reflected waves. These findings show the importance of utilizing two-way coupling when simulating the gas atomizer with a realistic heavy particle loading [17]. of utilizing two-way coupling when simulating the gas atomizer with a realistic heavy particle loading [17].

Breakup Behavior
According to the research of Mates and Settles [5], the particle size of the atomized powder is mainly determined by the secondary breakup process and the primary breakup in large ligaments (dligaments/djet = 0.1~1 [34]). Therefore, the initial diameters of the inner parent droplets were set to 0.4~4.0 mm and the droplets were injected at the top corner of the melt nozzle. In this model, when the melt temperature is lower than the melting point of 1520 K, it is considered solidified and no more breakup occurs. Figure 9 shows the trajectories of the particle droplets produced by the KHRT models. The parent droplets breakup and are accelerated by supersonic gas jets, which produce a comet-like cloud of disintegrating molten metal droplets. As shown in Figure 9a, the velocities reach almost 200 m/s after acceleration, and some even reach 300 m/s. The size of the particles represents the diameter of the droplets and the color represents the speed in Figure 10. The droplets have a significant breakup effect between the recirculation zone and the high-speed jet interface, resulting in a large number of small droplets. These small droplets are accelerated by the high speed gas flow, and, converging below the recirculation zone, the broken small droplets accelerate to more than 143.3 m/s (corresponding to the yellow droplets), while the large droplets with insufficient breakup have a speed of 46.0~98.1 m/s (corresponding to the green droplets). The colors in Figure 9b represent different sizes of droplets. The red parent droplets undergo a fast breakup process after being released near the nozzle, which turns the red droplets (~1 mm) to yellow (~70 μm) and blue (~50 μm) within a short distance of 50 mm from the tube. From Figures 9c and 10, we can see that in this region, where the droplet size changes drastically, the temperature of the droplets decreases relatively slowly near the nozzle exit area. However, after the violent breakup, the temperature of the child droplets decreases more sharply because the small particle size results in a large specific surface area and a faster heat transfer rate, which will increase the cooling rate [35,36]. However, in the downstream region, the color map of the diameters of the droplets shows that the diameters remain stable, indicating that the droplets are no longer substantially broken up in this region.

Breakup Behavior
According to the research of Mates and Settles [5], the particle size of the atomized powder is mainly determined by the secondary breakup process and the primary breakup in large ligaments (d ligaments /d jet = 0.1~1 [34]). Therefore, the initial diameters of the inner parent droplets were set to 0.4~4.0 mm and the droplets were injected at the top corner of the melt nozzle. In this model, when the melt temperature is lower than the melting point of 1520 K, it is considered solidified and no more breakup occurs. Figure 9 shows the trajectories of the particle droplets produced by the KHRT models. The parent droplets breakup and are accelerated by supersonic gas jets, which produce a comet-like cloud of disintegrating molten metal droplets. As shown in Figure 9a, the velocities reach almost 200 m/s after acceleration, and some even reach 300 m/s. The size of the particles represents the diameter of the droplets and the color represents the speed in Figure 10. The droplets have a significant breakup effect between the recirculation zone and the high-speed jet interface, resulting in a large number of small droplets. These small droplets are accelerated by the high speed gas flow, and, converging below the recirculation zone, the broken small droplets accelerate to more than 143.3 m/s (corresponding to the yellow droplets), while the large droplets with insufficient breakup have a speed of 46.0~98.1 m/s (corresponding to the green droplets). The colors in Figure 9b represent different sizes of droplets. The red parent droplets undergo a fast breakup process after being released near the nozzle, which turns the red droplets (~1 mm) to yellow (~70 µm) and blue (~50 µm) within a short distance of 50 mm from the tube. From Figures 9c and 10, we can see that in this region, where the droplet size changes drastically, the temperature of the droplets decreases relatively slowly near the nozzle exit area. However, after the violent breakup, the temperature of the child droplets decreases more sharply because the small particle size results in a large specific surface area and a faster heat transfer rate, which will increase the cooling rate [35,36]. However, in the downstream region, the color map of the diameters of the droplets shows that the diameters remain stable, indicating that the droplets are no longer substantially broken up in this region.
The particle diameter and temperature data for the particle trajectories extracted from Fluent is presented in Figure 11. Due to the Rosin-Rammler distribution of the initial parent droplets, in a short period of 0.5 ms, the diameters of the droplets were spread over a wide range of 20~2000 µm, and their temperatures remained high. After that, the particle size sharply reduced to 20~100 µm because of the breakup effect, and this was accompanied by a dramatic decrease in temperature. Finally, most of the child droplets flew away from the simulation area in 2 ms. To demonstrate the detailed change of the diameters of the droplets, we sampled droplets at different axial positions as a function of the radial distance from the nozzle axis ( Figure 12). x represents the distance between the droplets and the delivery tube, and from this we found that the diameters of the droplets at x = 2 mm vary significantly from 150~2000 µm, which means that only a few droplets had broken. At x = 4 mm, parent droplets appear together with small-sized child droplets, and increasingly more droplets breakup in this region. At x = 6~8 mm, the particle size decreased sharply to 20~200 µm, which is attributable to the very high gradient of the velocity between the gas and droplets; the RT instability increases faster and dominates the main breakup process. When x > 8 mm, the vertical distance of the droplets is gradually decreased to a negative value, and the diameter remains in the range of 20~100 µm, indicating that the droplets travel through the axis, and there is essentially no more breakup in the downstream area. In order to quantitatively evaluate the diameter reduction of the droplets during the atomizing process, we sampled droplets at the exit boundary and extracted their sizes. Detailed analysis is provided in the next section.
Processes 2020, 8, x FOR PEER REVIEW 10 of 19 Figure 9. The particle trajectories for the particle velocity magnitude; (a) particle velocity, (b) particle diameter, and (c) particle temperature for the KHRT breakup models. Figure 10. Partial enlarged view of droplet convergence area.
The particle diameter and temperature data for the particle trajectories extracted from Fluent is presented in Figure 11. Due to the Rosin-Rammler distribution of the initial parent droplets, in a short period of 0.5 ms, the diameters of the droplets were spread over a wide range of 20~2000 μm, and their temperatures remained high. After that, the particle size sharply reduced to 20~100 μm because Figure 9. The particle trajectories for the particle velocity magnitude; (a) particle velocity, (b) particle diameter, and (c) particle temperature for the KHRT breakup models.
Processes 2020, 8, x FOR PEER REVIEW 10 of 19 Figure 9. The particle trajectories for the particle velocity magnitude; (a) particle velocity, (b) particle diameter, and (c) particle temperature for the KHRT breakup models. Figure 10. Partial enlarged view of droplet convergence area.
The particle diameter and temperature data for the particle trajectories extracted from Fluent is presented in Figure 11. Due to the Rosin-Rammler distribution of the initial parent droplets, in a short period of 0.5 ms, the diameters of the droplets were spread over a wide range of 20~2000 μm, and their temperatures remained high. After that, the particle size sharply reduced to 20~100 μm because of the breakup effect, and this was accompanied by a dramatic decrease in temperature. Finally, most droplets is gradually decreased to a negative value, and the diameter remains in the range of 20~100 μm, indicating that the droplets travel through the axis, and there is essentially no more breakup in the downstream area. In order to quantitatively evaluate the diameter reduction of the droplets during the atomizing process, we sampled droplets at the exit boundary and extracted their sizes. Detailed analysis is provided in the next section.

Effect of Atomizing Parameters
The KHRT breakup model was used to extract the droplet sizes after the breakup process to assess the effects of the gas-to-melt mass flow ratio (GMR), operating pressure P, operating temperature T, half-taper angle α, protrusion length h, and gas slit nozzle diameter D. The statistical distribution of the droplets obtained from the model exit region will obtain the d50, d84/d50 and SPAN = (d90 − d10)/d50 of the droplets to evaluate the breakup effect of different parameters. during the atomizing process, we sampled droplets at the exit boundary and extracted their sizes. Detailed analysis is provided in the next section.

Effect of Atomizing Parameters
The KHRT breakup model was used to extract the droplet sizes after the breakup process to assess the effects of the gas-to-melt mass flow ratio (GMR), operating pressure P, operating temperature T, half-taper angle α, protrusion length h, and gas slit nozzle diameter D. The statistical distribution of the droplets obtained from the model exit region will obtain the d50, d84/d50 and SPAN = (d90 − d10)/d50 of the droplets to evaluate the breakup effect of different parameters.

Effect of Atomizing Parameters
The KHRT breakup model was used to extract the droplet sizes after the breakup process to assess the effects of the gas-to-melt mass flow ratio (GMR), operating pressure P, operating temperature T, half-taper angle α, protrusion length h, and gas slit nozzle diameter D. The statistical distribution of the droplets obtained from the model exit region will obtain the d 50 , d 84 /d 50 and SPAN = (d 90 − d 10 )/d 50 of the droplets to evaluate the breakup effect of different parameters.

Effect of GMR
The gas-to-melt mass flow ratio describes the gas energy input during atomization. It is a measure of the economic efficiency of the atomization equipment [36]. Utilizing the predictive ability of the model, the effect of GMR on d 50 was assessed. Figure 13 shows the particle size distribution and cumulative mass fraction at a melt mass flow rate of 0.16 kg/s. Like the atomized powders, the simulation data presents a typical normal distribution. To demonstrate the effect of GMR on particle size distribution, we further sampled the d 50 , then calculated the values of d 84 /d 50 and SPAN at different GMRs. As shown in Figure 14, the d 84 /d 50 refers to the standard deviation, SPAN = (d 90 − d 10 )/d 50 , and is used as a measure of the width of the particle size distribution. From this, we can see that a smaller melt mass flow rate resulted in smaller values of d 50 , d 84 /d 50 and SPAN. This is expected, higher GMR means that more gas jet energy input interacts with the same melt. This will obviously decrease the particle mean diameter and create a narrower spread of the distribution. different GMRs. As shown in Figure 14, the d84/d50 refers to the standard deviation, SPAN = (d90 − d10)/d50, and is used as a measure of the width of the particle size distribution. From this, we can see that a smaller melt mass flow rate resulted in smaller values of d50, d84/d50 and SPAN. This is expected, higher GMR means that more gas jet energy input interacts with the same melt. This will obviously decrease the particle mean diameter and create a narrower spread of the distribution.   that a smaller melt mass flow rate resulted in smaller values of d50, d84/d50 and SPAN. This is expected, higher GMR means that more gas jet energy input interacts with the same melt. This will obviously decrease the particle mean diameter and create a narrower spread of the distribution.    Figure 15 illustrates the effect of operating pressure on the d 50 process data and measure of standard deviation. As we can see in this simulation, with the increase of the operating pressure from 1.2 MPa to 2.7 MPa, the median diameter and standard deviation decrease significantly. This means that increasing the atomization pressure not only helps to obtain a fine powder, but also decreases the spread of the particle sizes. However, when the atomization pressure is increased to~3 MPa, the median diameter and standard deviation start to maintain stability, a phenomenon that has also been observed in actual experiments [37]. Although high pressure promotes the full breakup of droplets, collision and fusion among broken droplets become stronger, which forms larger-sized particles and satellite powders; therefore, fine powder exhibits little change at a high gas pressure. Essentially, increasing the operating pressure has the same effects on particle size distribution as does the large GMR; they both enlarge the gas energy input during atomization. observed in actual experiments [37]. Although high pressure promotes the full breakup of droplets, collision and fusion among broken droplets become stronger, which forms larger-sized particles and satellite powders; therefore, fine powder exhibits little change at a high gas pressure. Essentially, increasing the operating pressure has the same effects on particle size distribution as does the large GMR; they both enlarge the gas energy input during atomization.  Figure 16 represents the relationship of the operating gas temperature with the d50 process data and the measure of standard deviation. A similar trend can be identified for increases in both the operating pressure and operating temperature. High-temperature gas brings about the formation of finer particles as compared to a cold gas temperature due to higher kinetic energy, and therefore contributes to a greater gas velocity. However, the standard deviation and SPAN exhibits no obvious change in different gas temperatures. From the data, it is evident that the standard deviation values fluctuate around 1.2, which is significantly lower than the standard deviation of classical atomized powder (usually 2.0) [36].  Figure 16 represents the relationship of the operating gas temperature with the d 50 process data and the measure of standard deviation. A similar trend can be identified for increases in both the operating pressure and operating temperature. High-temperature gas brings about the formation of finer particles as compared to a cold gas temperature due to higher kinetic energy, and therefore contributes to a greater gas velocity. However, the standard deviation and SPAN exhibits no obvious change in different gas temperatures. From the data, it is evident that the standard deviation values fluctuate around 1.2, which is significantly lower than the standard deviation of classical atomized powder (usually 2.0) [36].

Effect of the Half-Taper Angle
The half-taper angle α is the gas nozzle convergent angle towards the centerline. Figure 17 shows the effect of the half-taper angle on the d50 process data and the measure of standard deviation. From

Effect of the Half-Taper Angle
The half-taper angle α is the gas nozzle convergent angle towards the centerline. Figure 17 shows the effect of the half-taper angle on the d 50 process data and the measure of standard deviation. From the simulation results, we know that with the increase of the half-taper angle, the d 50 exhibits a linearly decreasing trend, and that a larger half-taper angle is beneficial for obtaining finer powders. However, the standard deviation and SPAN did not follow this reduction. Previous research has revealed that a larger half-taper angle also leads to a severe aspiration pressure turbulence of the recirculation zone, which may block the outflow of the molten metal melt [10].

Effect of the Half-Taper Angle
The half-taper angle α is the gas nozzle convergent angle towards the centerline. Figure 17 shows the effect of the half-taper angle on the d50 process data and the measure of standard deviation. From the simulation results, we know that with the increase of the half-taper angle, the d50 exhibits a linearly decreasing trend, and that a larger half-taper angle is beneficial for obtaining finer powders. However, the standard deviation and SPAN did not follow this reduction. Previous research has revealed that a larger half-taper angle also leads to a severe aspiration pressure turbulence of the recirculation zone, which may block the outflow of the molten metal melt [10].  Figure 18 shows the influence of the melt tube protrusion length on the normalized D 50 and standard deviation. We can clearly determine the effect of the protrusion length on the particle size distribution. When the protrusion length of the melt tube is 4~10 mm, with the increase of the protrusion length, the d 50 exhibits a tendency of unilateral decline, and the standard deviation and SPAN remains steady. However, when the protrusion length of the melt tube is 12 mm, both the d 50 and standard deviation increase sharply. From the results of the atomization flow field simulation study, we can obtain a lower aspiration pressure when the melt tube protrusion length is in the range of 6~10 mm. protrusion length, the d50 exhibits a tendency of unilateral decline, and the standard deviation and SPAN remains steady. However, when the protrusion length of the melt tube is 12 mm, both the d50 and standard deviation increase sharply. From the results of the atomization flow field simulation study, we can obtain a lower aspiration pressure when the melt tube protrusion length is in the range of 6~10 mm.

Effect of the Gas Slit Nozzle Diameter
The gas slit nozzle diameter refers to the diameter of the gas jet out of the nozzle exit. Figure 19 shows the influence of the gas slit nozzle diameter on the normalized d50 and standard deviation. In this research, gas slit nozzle diameters of 26.8 mm to 50.8 mm at intervals of 6 mm were simulated. It can be seen that as the gas slit nozzle diameter increases, so does the standard deviation and SPAN. However, the d50 value is reduced continuously as the gas slit nozzle diameter increases. This means that a larger gas slit nozzle diameter can result in finer powders, but a wider size range. However, it will also increase the flight distance of the gas, which can cause a reduction of gas kinetic energy.

Effect of the Gas Slit Nozzle Diameter
The gas slit nozzle diameter refers to the diameter of the gas jet out of the nozzle exit. Figure 19 shows the influence of the gas slit nozzle diameter on the normalized d 50 and standard deviation. In this research, gas slit nozzle diameters of 26.8 mm to 50.8 mm at intervals of 6 mm were simulated. It can be seen that as the gas slit nozzle diameter increases, so does the standard deviation and SPAN. However, the d 50 value is reduced continuously as the gas slit nozzle diameter increases. This means that a larger gas slit nozzle diameter can result in finer powders, but a wider size range. However, it will also increase the flight distance of the gas, which can cause a reduction of gas kinetic energy.  Figure 19. The effect of gas slit nozzle diameter on d50, d84/d50 and SPAN for the KHRT breakup model simulation.

Discussion
According to the above analysis, it is evident that the two-way coupling between melt droplets

Discussion
According to the above analysis, it is evident that the two-way coupling between melt droplets and gas flow has a huge impact on the gas flow field during the atomization breakup process; periodic reflected waves were replaced by a ring-wrapped gas flow. In the range of x = 4~8 mm, near the nozzle top corner where the RT instability dominates the breakup behavior, the released parent droplets will be broken up in a short period time (0.5 ms). The broken droplets then flow downstream with the gas flow to form an hourglass molten droplet cloud, and most of the child droplets fly away from the simulation area in 2 ms.
We then explored the effects of the atomizing process conditions and nozzle geometry parameter on the particle size distribution for the KHRT model. Generally speaking, we intended to obtain finer particles and a concentrated particle size distribution with a low standard deviation in an industrial manufacturing process. For the atomization structure parameters, the above numerical simulation data suggest that a larger half-taper angle α results in finer metal powders. However, gas flow field research shows that it also increased the aspiration pressure in the recirculation zone, which would block the outflow of the metal melt. Therefore, it is more appropriate to use the median value of 30 • . The protrusion length range of 8~10 mm of the delivery tube h was found to be better for obtaining finer powders. Additionally, although increasing the gas slit nozzle diameter D can reduce the average particle size, it will also simultaneously increase the standard deviation and enlarge the gas flow field; therefore, we tend to choose a smaller gas slit nozzle diameter, such as 32.8 mm. For the operating parameter, although increasing the operating gas temperature can effectively reduce the average particle size, it demands complex extra equipment and costs a significant amount of energy in industrial production. Therefore, it must be determined according to the actual situation. We also examined the influence of atomizing conditions on the particle size distribution, which can provide an effective basis for parameter selection in industrial production. Typically, as powder must be supplied with tailored distributions to meet the requirements of its intended use (e.g., laser cladding demands an average particle size in the range of 80~95 µm), it is essential to minimize this off-size powder inventory to ensure a powder producer remains competitive. Therefore, we can roughly determine the range of operating pressure to be 2 MPa and the melt mass flow rate to be 0.6 kg/s by the simulation of the operating conditions. This work provides a method to design a de Laval-type gas nozzle atomizer. The specified FeNiCrBSiNbC alloy secondary breakup process was simulated by the validated KHRT model. Research on the atomization structural parameters and atomizing conditions can be used for the design of the atomizing nozzles and for optimizing their parameter selection in industrial production. In addition, it will contribute to the understanding of the principle of the secondary breakup process in the de Laval-type nozzle. However, many problems remain to be solved, such as improving the reliability of the model and studying the entire process of atomization, rather than dividing it into primary and secondary breakup processes. At the same time, the predictive capability of the KHRT model need be proven by experiments, which is our main task in future work.

Conclusions
A de Laval-type gas outlet atomizer was designed and studied using CFD software. To assess the interaction between the gas jet and molten metal, we used the Euler-Lagrange method to describe the gas flow field and metal droplets. The validated KHRT model was used for the quantitative assessment of the particle size distribution. From this simulation, we can see that an hourglass cloud of child droplets was generated; a far more realistic particle track than Zeoli [14] was produced by this breakup model. The characteristics of the flow field with the particle load demonstrate that the two-way coupling was vital for approximating the impact of the particles on the behavior of the fluid flow. This simulation approach was used for the quantitative evaluation of the effects of altering the atomizing process condition and nozzle geometry on the particle size distribution of the produced powders. From this study, we have drawn the following conclusions: