# Numerical Investigation of Droplet Properties of a Liquid Jet in Supersonic Crossflow

Authors: Yu-Qi Wang; Feng Xiao; Sen Lin; Yao-Zhi Zhou
Journal: International Journal of Aerospace Engineering
 (2021)

---

## Abstract

The atomization process of a liquid jet in supersonic crossflow with a Mach number of 1.94 was investigated numerically under the Eulerian-Lagrangian scheme. The droplet stripping process was calculated by the KH (Kelvin-Helmholtz) breakup model, and the secondary breakup due to the acceleration of shed droplets was calculated by the combination of the KH breakup model and the RT (Rayleigh-Taylor) breakup model. In our research, the existing KH-RT model was modified by optimizing the empirical constants incorporated in this model. Moreover, it was also found that the modified KH-RT breakup model is applied better to turbulent inflow of a liquid jet than laminar inflow concluded from the comparisons with experimental results. To validate the modified breakup model, three-dimensional spatial distribution and downstream distribution profiles of droplet properties of the liquid spray in theMa=1.94 airflow were successfully predicted in our simulations. Eventually, abundant numerical cases under different operational conditions were launched to investigate the correlations of SMD (Sauter Mean Diameter) with the nozzle diameter as well as the airflow Mach number, and at the same time, modified multivariate power functions were developed to describe the correlations.


---

## Body

## 1. Introduction

Within the combustion chamber of the Scramjet engine, a transverse liquid jet is injected into the supersonic airflow at a certain flow rate. Once exposed to the incoming airflow with a high Mach number, the jet deforms immediately and breaks apart into small fragments rapidly under strong aerodynamic forces [1]. None of any clear boundaries can be utilized to divide the whole breakup process into several specific stages. However, generally, the process of the liquid column breaking apart into pieces of initially large droplets can be recognized as “the primary breakup,” while the process of initially large droplets further breaking up into small-sized droplets is usually recognized as “the secondary breakup.” Gorokhovski [2] pointed out that the primary breakup has a significant impact on the subsequent droplet formation and dispersion, and yet, it is of huge challenge to understand the whole physical mechanisms within for various complicated structures in the flow field as well as the large time and space spans which thus make it even harder to investigate the microstructures both experimentally and numerically. On the other hand, since our goal is to produce small droplets downstream with appropriate sizes to increase the evaporation and mixing rates of the liquid fuel, one of the most significantly practical motivations to investigate the atomization process is to determine the conditions in which the desired final fragment sizes can be acquired. As noted by Tryggvason [3], the highest airflow velocity does not always guarantee the smallest droplet diameter. Therefore, the secondary breakup process should also be clearly understood to meet the actual engineering demands.Plenty of experimental investigations have been launched regarding the liquid jet in supersonic crossflow. With the application of laser technology and computer technology, many advanced testing technologies have been developed, such as phase Doppler particle analyzer (PDA/PDPA/LDV), laser scattering technology/laser holography technology, laser-induced fluorescence technology (PLIF), and particle image velocity instrument (PIV). The liquid spray properties measured by different instuments turn out differently. Koh [4] found that the mass distribution of PDA is obviously smaller than that of optical imaging. Lin [5] contrasted the effects on the penetration height caused by different optical approaches and pointed out that PDA is more sensitive to thin liquid mist thus measuring a higher penetration height than other approaches such as high-speed photography and schlieren. Lin [6] studied the structures of water jets injected into a Ma=1.94 crossflow by utilizing a two-component PDPA and successfully discovered the S-type distribution profiles of droplet size in the downstream position. What is more, the correlations of the operational conditions and droplet properties can be investigated, and the whole physical mechanisms of the atomization process can be further understood by numerical calculations. Currently, numerical approaches to investigate the whole gas-liquid mixing and interacting process in supersonic conditions can be generally divided into two types. One is to capture the movement of gas-liquid interface based on the Eulerian scheme such as VOF (volume of fluid) and LS (level set) approaches. These approaches are at high accuracy but huge computational costs. The other is to track the position of single droplets by integration of time based on the Lagrangian scheme. In this paper, the lagrangian method is focused on because compared with the Eulerian scheme, it can remarkably reduce computation costs and combine itself with breakup models compatibly; although, the initial size and velocity distributions of large droplets need to be determined.In the beginning stage of the investigation, out of the restriction of experimental conditions, most of the research was focused on the dependence of deformation and breakup of droplets on those dimensionless numbers such as Reynolds number and Weber number [7–9]. At present, it is still extremely difficult to simulate the whole atomization process from continuous liquid column to uniform spray plume with one single physical model. Therefore, existing breakup models have been coupled and improved to calculate breakup processes of liquid jets to acquire good agreement with experimental results. Currently, the theories of surface waves induced by the Kelvin-Helmholtz (KH) and the Rayleigh-Taylor (RT) instabilities are considered to explain the essential mechanisms of the liquid breakup process mostly [10]. The KH-RT breakup model is based on the surface wave theories and takes both KH and RT instabilities into account. At present, it is considered as the most appropriate breakup model to truly reveal the breakup mechanisms of the liquid phase in the supersonic environment. Liu [11] studied the breakup parameters of the KH-RT model and pointed out that in the supersonic environment, parameters controlling the breakup time have quite limited impacts for the very short breakup time in supersonic airflows. Yang et al. (Yang, Zhu, Sun, and Chen, 2017) [12] improved the KH-RT breakup model by considering the compressible effects of the gaseous environment and found that the penetration height of the improved model agreed well but the spread and size distributions of the liquid spray were still different from the experiment results. Li et al. (Li, Wang, Sun, and Wang, 2017) [13] adopted the KH breakup model to simulate the droplet stripping process near the nozzle and coupled RT breakup and TAB (Taylor analogy breakup) models to simulate the secondary breakup of droplets. In his research, a LES code program of two-phase flow was carried out under high resolution grid, and the downstream atomization characteristic profiles obtained were in good agreement with the experimental results.The inflow turbulence of a liquid jet was widely investigated both experimentally and numerically. It was noted to control the primary breakup process of a liquid jet in supersonic crossflow and thus further control the downstream droplet properties. Mazallon et al. (Mazallon, Dai, and Faeth, 1999) [14] and Sallam et al. [15] (Sallam, Ng, Sankarakrishnan, Aalburg, and Lee, 2006) investigated the laminar liquid jet injected into a uniform gaseous crossflow by using pulsed shadowgraphy and pulsed holography. It was found that surface waves formed on the upstream side of the laminar jet column, and that the wavelength decreased with the increase of the gaseous Weber number which indicated that the surface waves in the primary breakup process originate from RT instability. Xiao et al. [16] (Xiao, Dianat, and McGuirk, 2013) also drew the same conclusion with a two-phase-flow large-eddy simulation. It was also concluded that it is liquid rather than gaseous turbulence that determines the initial liquid-jet instability and interface characteristics. Lee et al. [17] (Lee, Aalburg, Diez, Faeth, and Sallam, 2007) discovered that the SMD of droplets stripped from the surface of a turbulent liquid jet column was not influenced by the crossflow, and thus the liquid turbulence controls the primary breakup process.In this paper, the atomization process of a liquid jet in supersonic crossflow with a Mach number of 1.94 was investigated numerically using the Lagrangian method. In this method, large droplets with initial size and velocity distributions were injected into a supersonic crossflow in substitution for a continuous liquid jet. Then, the droplets broke apart into small fragments calculated with the modified KH-RT breakup model and eventually formed the spray plume. This paper is organized as follows: In Section2, relevant mathematical models are presented or modified including the dynamic equations, breakup models, and states of internal nozzle flow. In Sections 3 and 4, computational conditions are introduced, and our modified models are verified via comparing our numerical results with experiment. Furthermore, the correlation functions of SMD and two important parameters are summarized. Eventually, several important conclusions are presented in Section 5.


## 2. Mathematical Models

### 2.1. Eulerian-Lagrangian Governing Equations

#### 2.1.1. Gas-Phase Dynamic Equation

(1)du⟶dt=g⟶−∇Pρ+v∇2u⟶+u⟶−u⟶pτr+Fothers,(2)τr=ρpdp218μ24CdRe,(3)Re=ρdpu⟶−u⟶pμ,where ρ is the gas density, u⟶ is the gas velocity, P is the static pressure, v is the gas dynamic viscosity, F⟶ is the momentum source term of droplets, Cd is the drag coefficient, u⟶p is the droplet velocity, ρp is the droplet density, dp is the droplet diameter, Fother is other forces per unit gas mass, τr is the droplet relaxation time [18], and Re is the relative Reynolds number .


#### 2.1.2. Liquid-Phase Dynamic Equation

(4)du⟶pdt=u⟶−u⟶pτr+g⟶ρp−ρρp+F⟶′where u⟶−u⟶p/τr is the drag force per unit droplet mass, and F⟶′ is an additional acceleration per unit droplet mass term, such as “virtual mass” force(5)F⟶′=Cvmρρpu⟶p∇u⟶−du⟶pdt.Integration of time in Equation (4) yields the velocity of the droplet at each point along the trajectory, with the trajectory itself predicted by(6)dxdt=up.The new locationxpn+1 can be computed from(7)xpn+1=xpn+Δtun+aτr+τr1−e−Δtτrupn−un−aτr,where a includes accelerations due to all other forces except drag force, and upn and nn represent droplet velocity and gas velocity at the old location.


### 2.2. Turbulence Model

In this paper, thek−ω−SST(Shear-stress Transport) model is used for its high prediction ability in the far-field and near-wall regions. In RANS simulations, the instantaneous quantity f is splitted into a mean f¯ and fluctuating f′ components f=f¯+f′. For compressible flow, the density ρ varies so widely that a mass-weighted average f~ (called Farve average) is usually preferred f~=ρf/ρ¯. Any quantity f may be splinted into mean and fluctuating components as(8)f=f~+f′withf~″=0.The average balance equations of continuous phase are written as(9)∂ρ¯∂t+∂ρ¯u~i∂xi=0,(10)∂ρ¯u~i∂t+∂ρ¯u~iu~j∂xj=−∂p¯δij∂xi+∂t¯ij−τij∂xj,(11)∂ρ¯E∂t+∂ρ¯Eu~j∂xj=−∂−qLj−qTj∂xj+∂u~jt¯ij+τij∂xj,(12)p¯=ρ¯RT~,where E=e~+uiui/2+k,H=h~+uiui/2+k,e~=cpT~,h~=cpT~cv′cp are specific heat capacity of constant volume and pressure, respectively. The turbulent kinetic energy k is written as(13)k=ρui′ui′¯/ρ,(14)t¯=2μSij−23∂u~k∂xkδij,(15)qLi=−μPrL∂h~∂xj,qTi=−KPrT∂h~∂xj,where μ is the molecular viscosity, and K is thermal conductivity.To provide closures for the unknown turbulent kinetic energyk, Menter [19] proposed the k−ω−SST model. The classical k−ε turbulence model has high prediction ability with high Reynolds number region in far field while the classical k−ω turbulence model has better prediction ability and more stable numerical properties with the near-wall region. The k−ω−SST turbulence model combines the advantages of both k−ε and k−ω turbulence models switched by a mixing function F. The transport equations of k−ω−SST turbulence model are as follows:(16)∂k∂t+Ui∂k∂xi=∂∂xiv+σkvt∂k∂xi+pk−Cuωk,(17)∂ω∂t+Ui∂ω∂xi=∂∂xiv+σωvt∂ω∂xi+γωkpk−βω2+1−F12σωω∂k∂xi∂ω∂xi.The mixing functionF1 is defined by(18)F1=tanharg14,(19)arg1=minmaxkCuwy,500vy2w,4ρσw2kCDkwy2,where y is the vertical distance from the wall surface, Dkw is the positive part of the cross diffusion term, and Dkw is determined by(20)CDkw=max2ρσw21w∂k∂xi∂w∂xi,10−20.All the coefficientsσk,σw,β,γ can be calculated in a uniform form by(21)ϕ=F1ϕ1+1−F1ϕ2,where ϕ1′ϕ2′ is the corresponding coefficients of k−ε and k−ω turbulence models, respectively.The eddy viscosity coefficientvt is defined by(22)vt=a1kmaxa1,ΩF2,where F2=tanharg22,Ω is the eddy vector.(23)arg2=max2kCμω,500vy2ω.In Eq. (23), the first component is the eddy viscosity coefficient of the k−ω turbulence model, and the second component is acquired from the one-equation turbulence model based on the characteristics of shear stress in laminar boundary layers.


### 2.3. Modified KH-RT Breakup Model

The KH-RT breakup model is based on the linearized instability theory and takes both KH and RT breakup into account. It is often used to simulate high Weber number sprays. Within this model, a length-limited liquid core in the near nozzle region is assumed as shown in Figure1. The droplet stripping process within the liquid core is considered by KH instability, and the acceleration of shed droplets was calculated by the competition of both KH and RT instabilities.Figure 1
Schematic diagram of the liquid core.For the KH breakup model, the propagation equations of unstable waves on the surface of cylindrical jet are calculated numerically, and the maximum growth rateΩKH corresponding to wave length ΛKH is obtained as below [20].(24)ΛKHr=9.021+0.45Oh0.51+0.4T0.71+0.87Weg1.670.6,(25)ΩKHρdr3σ0.5=0.34+0.38Weg1.51+Oh1+1.470.6,where Oh=We/Red, Wed=ρdrv2/σ, Weg=ρgrv2/σ, Red=rv/υd, and T=OhWeg.The radius of child droplet stripped from the cylindrical jet is(26)r′=B0ΛKH,(27)drdt=−r−rτKH,τKH is the KH breakup time calculated by(28)τKH=3.726B1rΩKHΛKH,where B0 is the KH breakup radius constant, and B1 is the KH breakup time constant.For the RT breakup model, the size of child droplet and breakup time depends on the fastest growing wave. The wave length of the fastest growing wave is calculated by [21, 22].(29)ΛRT=2πC13σadρd−ρg.The RT breakup child droplet radius related to the wave length of the fastest growing wave is calculated by(30)r″=0.5ΛRT.The RT breakup time is calculated as below(31)τRT=C2σ0.5ρd+ρg23adρd+ρg1.5,where αd is the acceleration of the droplet, C1 is the RT breakup radius constant, and C2 is the RT breakup time constant.When the KH-RT breakup model is adopted to track wave growth on the surface of droplets, it is often beneficial to adjust the constants to acquire optimal simulation results compared with the experiments. According to the research by Yang et al. (Yang, Le, He, and Chen, 2012) [23], the optimal values of breakup constants B0, B1, C1, and C2 in supersonic flows are summarized below as Table 1.Table 1
Optimal KH-RT breakup model constants from Yang’s research.
ConstantB0B1C1C2Optimal value0.151.730.51On the basis of Yang’s work, the optimal KH-RT breakup model constants were further investigated in this paper. Figure2 illustrates that the downstream SMD distribution is mainly controlled by the KH radius constant B0 rather than time constant B1. It is also found that the RT radius constant C1 has much less influence on the downstream SMD than the KH radius constant B0. It can be explained that in the KH-RT breakup model, the KH instability participates both the stripping and acceleration breakup process, and that the breakup time for a droplet in the supersonic environment is extremely short; thus, the impacts of the time constants can be neglected. Typically, the RT instability grows faster when droplet acceleration is high, and this effect dominates for high Weber number sprays. Therefore, the cause of the little control shown by C1 is probably because the present liquid core length is larger than expected.Figure 2
SMD distribution in the free-stream direction with varied KH breakup model constants: (a)B0 and (b) B1.
(a)(b)Therefore, using the idea of controlling variables, the optimal values of breakup constantsB0,B1,C1,C2 for supersonic simulation of two-phase interactions calculated by FLUENT software are shown in Table 2.Table 2
Improved results of optimal KH-RT breakup model constants.
ConstantB0B1C1C2Optimal value0.31.730.51


### 2.4. Inflow Turbulence of a Liquid Jet

In this paper, the DPM (discrete particle model) is used under the Eulerian-Laglangian frame via the fluid mechanic simulation software FLUENT. In the DPM model, the “Blob” model was applied by injecting large droplets into the free airflow in substitution for a continuous liquid jet. Then, the behavior and trajectory of droplets are calculated by breakup models and interaction with the crossflow. Although with the exact same orifice diameter specified, the inflow turbulence of a liquid jet is discovered to have significant effects over the primary breakup process by affecting the initial droplet sizes and velocities. Therefore, to accurately predict the spray characteristics, the correct turbulent state of the internal nozzle flow must be specified. As can be seem from Figure3, the internal nozzle flow can be divided into single-phase flow, cavitating flow, and flipped flow, respectively, with the intensity of turbulence decreasing in turn. And the turbulence within liquid jet is mainly determined by Reynolds number.(32)Rel=dρlμ2p1−p2ρ1,where p1,p2 is the controllable upstream and downstream pressure, respectively.Figure 3
States of internal nozzle flow: (a) single-phase, (b) cavitating, and (c) flipped.
(a)(b)(c)In the rest of this section, the varied initial droplet size and velocity of different flow turbulence are focused on. And it is shown that with the same nozzle diameter and flow mass rate specified, the turbulent (single-phase) liquid flow fits better with the experiment results than the laminar (flipped) liquid flow.


#### 2.4.1. Determination of the Initial Droplet Size Distribution

According to our investigation, the initial droplet diameter distribution is closely related to the nozzle state. To indicate the connection, a two-parameter RR (Rosin-Rammler) distribution is used to represent the droplet diameter distribution as below, characterized by the most probable droplet sized0 and a spread parameter n. The mass fraction of droplets with diameters greater than d is calculated by(33)Φd=e−d/d0n.The first parameter required to specify the droplet size distribution is the most probable droplet sized0. For a single-phase nozzle flow, the correlation of Wu et al. (Wu, Tseng, and Faeth, 1992) [24] is applied to calculate SMD d32 considering the initial drop size related to the turbulence quantities of the liquid jet. Snyder [25] gives the most general relationship between SMD and most probable diameter for a Rosin-Rammler distribution.(34)d32=133.0λWe−e0.74,(35)We=ρlu2λσ,(36)d0=1.2726d321−121/n,where λ=d/8,λ is the radial integral length scale at the jet exit based upon fully developed turbulent pipe flow, and We is the Weber number of the liquid jet.For a cavitating nozzle flow, the correlation of Wu can still be applied, and yet the length scale for a cavitating nozzle isλ=deff/8, where deff is the effective diameter of the exiting liquid jet according to Schmidt and Corradini [26].For the case of a flipped nozzle flow, the initial droplet diameter is set to the diameter of the liquid jet(37)d0=dCct,where Cct is a theoretical constant equal to 0.611, which comes from potential flow analysis of flipped nozzles.The second parameter required to specify the droplet size distribution is the spread parametern. The values for the spread parameter are determined from past modeling experience and experimental observations. Table 3 lists the values of n for three flow states. The larger the value of the spread parameter, the narrower the droplet size distribution.Table 3
Theoretical values of spread parametern for different liquid states.
StateSpread parameterSingle-phase3.5Cavitating1.5Flipped∞Having specified the most probable diameter and the spread parameter, the initial droplet size distribution can thus be determined. It should be noted that the actual size distribution may be a little different with the theoretical one for the limited number of droplets injected from the nozzle exit. Figure4 shows the initial droplet size distribution of three nozzle flow states. It can be seen that the size of injected droplets of turbulent internal flow tends to be smaller and more uniform while the laminar internal flow tends to be larger and more centralized. It has been experimentally figured out that the droplet size distribution in the near nozzle region has a tremendous effect on downstream droplet properties.Figure 4
Initial droplet size distribution of different flow states.Figure5 illustrates the SMD distribution at x=50,100mm in the central plane distinguished by the flow turbulence with different initial droplet diameter distributions. Therefore, the turbulent liquid jet fits better with the experimental results in the aspect of the droplet size distribution.Figure 5
Downstream SMD distributions with different flow state: (a) turbulent and (b) laminar.
(a)(b)


#### 2.4.2. Determination of Initial Droplet Velocity

For a single-phase nozzle, the estimate of exit velocityu comes from the conservation of mass and the assumption of a uniform exit velocity:(38)u=meffρlA.For a cavitating nozzle, an expression for a higher velocity over a reduced area derived by Schmidt and Corradini [26] is presented here instead of a uniform exit velocity:(39)u=2Ccp1−p2+1−2CcpvCc2ρlρ1−pv,where Cc is the contraction coefficient by Nurick’s [27] fit, p1 is the internal upstream pressure of the nozzle, p2 is the internal downstream pressure of the nozzle, and pv is the vapor pressure of the nozzle.For a flipped nozzle, the exit velocity is derived from the conservation of mass and the value of the reduced flow area:(40)u=meffρlCctA.The tremendous effect on the initial droplet velocity directly affects the penetration height of the liquid spray. According to our investigation, the penetration height of a turbulent liquid jet fits better with experimental results [6] than that of a laminar liquid jet in the same gaseous conditions as shown in Figure 6. It can be explained that the laminar liquid jet has a larger initial droplet velocity with the mass flow rate specified due to the flow contraction in the nozzle, which is inconsistent with the actual situation. Figure 7 indicates the higher velocity of laminar flow compared with turbulent flow.Figure 6
Comparison of spray penetration height with experimental results in the central plane: (a) turbulent and (b) laminar.
(a)(b)Figure 7
Calculation domain and meshing.


## 2.1. Eulerian-Lagrangian Governing Equations

### 2.1.1. Gas-Phase Dynamic Equation

(1)du⟶dt=g⟶−∇Pρ+v∇2u⟶+u⟶−u⟶pτr+Fothers,(2)τr=ρpdp218μ24CdRe,(3)Re=ρdpu⟶−u⟶pμ,where ρ is the gas density, u⟶ is the gas velocity, P is the static pressure, v is the gas dynamic viscosity, F⟶ is the momentum source term of droplets, Cd is the drag coefficient, u⟶p is the droplet velocity, ρp is the droplet density, dp is the droplet diameter, Fother is other forces per unit gas mass, τr is the droplet relaxation time [18], and Re is the relative Reynolds number .


### 2.1.2. Liquid-Phase Dynamic Equation

(4)du⟶pdt=u⟶−u⟶pτr+g⟶ρp−ρρp+F⟶′where u⟶−u⟶p/τr is the drag force per unit droplet mass, and F⟶′ is an additional acceleration per unit droplet mass term, such as “virtual mass” force(5)F⟶′=Cvmρρpu⟶p∇u⟶−du⟶pdt.Integration of time in Equation (4) yields the velocity of the droplet at each point along the trajectory, with the trajectory itself predicted by(6)dxdt=up.The new locationxpn+1 can be computed from(7)xpn+1=xpn+Δtun+aτr+τr1−e−Δtτrupn−un−aτr,where a includes accelerations due to all other forces except drag force, and upn and nn represent droplet velocity and gas velocity at the old location.


## 2.1.1. Gas-Phase Dynamic Equation

(1)du⟶dt=g⟶−∇Pρ+v∇2u⟶+u⟶−u⟶pτr+Fothers,(2)τr=ρpdp218μ24CdRe,(3)Re=ρdpu⟶−u⟶pμ,where ρ is the gas density, u⟶ is the gas velocity, P is the static pressure, v is the gas dynamic viscosity, F⟶ is the momentum source term of droplets, Cd is the drag coefficient, u⟶p is the droplet velocity, ρp is the droplet density, dp is the droplet diameter, Fother is other forces per unit gas mass, τr is the droplet relaxation time [18], and Re is the relative Reynolds number .


## 2.1.2. Liquid-Phase Dynamic Equation

(4)du⟶pdt=u⟶−u⟶pτr+g⟶ρp−ρρp+F⟶′where u⟶−u⟶p/τr is the drag force per unit droplet mass, and F⟶′ is an additional acceleration per unit droplet mass term, such as “virtual mass” force(5)F⟶′=Cvmρρpu⟶p∇u⟶−du⟶pdt.Integration of time in Equation (4) yields the velocity of the droplet at each point along the trajectory, with the trajectory itself predicted by(6)dxdt=up.The new locationxpn+1 can be computed from(7)xpn+1=xpn+Δtun+aτr+τr1−e−Δtτrupn−un−aτr,where a includes accelerations due to all other forces except drag force, and upn and nn represent droplet velocity and gas velocity at the old location.


## 2.2. Turbulence Model

In this paper, thek−ω−SST(Shear-stress Transport) model is used for its high prediction ability in the far-field and near-wall regions. In RANS simulations, the instantaneous quantity f is splitted into a mean f¯ and fluctuating f′ components f=f¯+f′. For compressible flow, the density ρ varies so widely that a mass-weighted average f~ (called Farve average) is usually preferred f~=ρf/ρ¯. Any quantity f may be splinted into mean and fluctuating components as(8)f=f~+f′withf~″=0.The average balance equations of continuous phase are written as(9)∂ρ¯∂t+∂ρ¯u~i∂xi=0,(10)∂ρ¯u~i∂t+∂ρ¯u~iu~j∂xj=−∂p¯δij∂xi+∂t¯ij−τij∂xj,(11)∂ρ¯E∂t+∂ρ¯Eu~j∂xj=−∂−qLj−qTj∂xj+∂u~jt¯ij+τij∂xj,(12)p¯=ρ¯RT~,where E=e~+uiui/2+k,H=h~+uiui/2+k,e~=cpT~,h~=cpT~cv′cp are specific heat capacity of constant volume and pressure, respectively. The turbulent kinetic energy k is written as(13)k=ρui′ui′¯/ρ,(14)t¯=2μSij−23∂u~k∂xkδij,(15)qLi=−μPrL∂h~∂xj,qTi=−KPrT∂h~∂xj,where μ is the molecular viscosity, and K is thermal conductivity.To provide closures for the unknown turbulent kinetic energyk, Menter [19] proposed the k−ω−SST model. The classical k−ε turbulence model has high prediction ability with high Reynolds number region in far field while the classical k−ω turbulence model has better prediction ability and more stable numerical properties with the near-wall region. The k−ω−SST turbulence model combines the advantages of both k−ε and k−ω turbulence models switched by a mixing function F. The transport equations of k−ω−SST turbulence model are as follows:(16)∂k∂t+Ui∂k∂xi=∂∂xiv+σkvt∂k∂xi+pk−Cuωk,(17)∂ω∂t+Ui∂ω∂xi=∂∂xiv+σωvt∂ω∂xi+γωkpk−βω2+1−F12σωω∂k∂xi∂ω∂xi.The mixing functionF1 is defined by(18)F1=tanharg14,(19)arg1=minmaxkCuwy,500vy2w,4ρσw2kCDkwy2,where y is the vertical distance from the wall surface, Dkw is the positive part of the cross diffusion term, and Dkw is determined by(20)CDkw=max2ρσw21w∂k∂xi∂w∂xi,10−20.All the coefficientsσk,σw,β,γ can be calculated in a uniform form by(21)ϕ=F1ϕ1+1−F1ϕ2,where ϕ1′ϕ2′ is the corresponding coefficients of k−ε and k−ω turbulence models, respectively.The eddy viscosity coefficientvt is defined by(22)vt=a1kmaxa1,ΩF2,where F2=tanharg22,Ω is the eddy vector.(23)arg2=max2kCμω,500vy2ω.In Eq. (23), the first component is the eddy viscosity coefficient of the k−ω turbulence model, and the second component is acquired from the one-equation turbulence model based on the characteristics of shear stress in laminar boundary layers.


## 2.3. Modified KH-RT Breakup Model

The KH-RT breakup model is based on the linearized instability theory and takes both KH and RT breakup into account. It is often used to simulate high Weber number sprays. Within this model, a length-limited liquid core in the near nozzle region is assumed as shown in Figure1. The droplet stripping process within the liquid core is considered by KH instability, and the acceleration of shed droplets was calculated by the competition of both KH and RT instabilities.Figure 1
Schematic diagram of the liquid core.For the KH breakup model, the propagation equations of unstable waves on the surface of cylindrical jet are calculated numerically, and the maximum growth rateΩKH corresponding to wave length ΛKH is obtained as below [20].(24)ΛKHr=9.021+0.45Oh0.51+0.4T0.71+0.87Weg1.670.6,(25)ΩKHρdr3σ0.5=0.34+0.38Weg1.51+Oh1+1.470.6,where Oh=We/Red, Wed=ρdrv2/σ, Weg=ρgrv2/σ, Red=rv/υd, and T=OhWeg.The radius of child droplet stripped from the cylindrical jet is(26)r′=B0ΛKH,(27)drdt=−r−rτKH,τKH is the KH breakup time calculated by(28)τKH=3.726B1rΩKHΛKH,where B0 is the KH breakup radius constant, and B1 is the KH breakup time constant.For the RT breakup model, the size of child droplet and breakup time depends on the fastest growing wave. The wave length of the fastest growing wave is calculated by [21, 22].(29)ΛRT=2πC13σadρd−ρg.The RT breakup child droplet radius related to the wave length of the fastest growing wave is calculated by(30)r″=0.5ΛRT.The RT breakup time is calculated as below(31)τRT=C2σ0.5ρd+ρg23adρd+ρg1.5,where αd is the acceleration of the droplet, C1 is the RT breakup radius constant, and C2 is the RT breakup time constant.When the KH-RT breakup model is adopted to track wave growth on the surface of droplets, it is often beneficial to adjust the constants to acquire optimal simulation results compared with the experiments. According to the research by Yang et al. (Yang, Le, He, and Chen, 2012) [23], the optimal values of breakup constants B0, B1, C1, and C2 in supersonic flows are summarized below as Table 1.Table 1
Optimal KH-RT breakup model constants from Yang’s research.
ConstantB0B1C1C2Optimal value0.151.730.51On the basis of Yang’s work, the optimal KH-RT breakup model constants were further investigated in this paper. Figure2 illustrates that the downstream SMD distribution is mainly controlled by the KH radius constant B0 rather than time constant B1. It is also found that the RT radius constant C1 has much less influence on the downstream SMD than the KH radius constant B0. It can be explained that in the KH-RT breakup model, the KH instability participates both the stripping and acceleration breakup process, and that the breakup time for a droplet in the supersonic environment is extremely short; thus, the impacts of the time constants can be neglected. Typically, the RT instability grows faster when droplet acceleration is high, and this effect dominates for high Weber number sprays. Therefore, the cause of the little control shown by C1 is probably because the present liquid core length is larger than expected.Figure 2
SMD distribution in the free-stream direction with varied KH breakup model constants: (a)B0 and (b) B1.
(a)(b)Therefore, using the idea of controlling variables, the optimal values of breakup constantsB0,B1,C1,C2 for supersonic simulation of two-phase interactions calculated by FLUENT software are shown in Table 2.Table 2
Improved results of optimal KH-RT breakup model constants.
ConstantB0B1C1C2Optimal value0.31.730.51


## 2.4. Inflow Turbulence of a Liquid Jet

In this paper, the DPM (discrete particle model) is used under the Eulerian-Laglangian frame via the fluid mechanic simulation software FLUENT. In the DPM model, the “Blob” model was applied by injecting large droplets into the free airflow in substitution for a continuous liquid jet. Then, the behavior and trajectory of droplets are calculated by breakup models and interaction with the crossflow. Although with the exact same orifice diameter specified, the inflow turbulence of a liquid jet is discovered to have significant effects over the primary breakup process by affecting the initial droplet sizes and velocities. Therefore, to accurately predict the spray characteristics, the correct turbulent state of the internal nozzle flow must be specified. As can be seem from Figure3, the internal nozzle flow can be divided into single-phase flow, cavitating flow, and flipped flow, respectively, with the intensity of turbulence decreasing in turn. And the turbulence within liquid jet is mainly determined by Reynolds number.(32)Rel=dρlμ2p1−p2ρ1,where p1,p2 is the controllable upstream and downstream pressure, respectively.Figure 3
States of internal nozzle flow: (a) single-phase, (b) cavitating, and (c) flipped.
(a)(b)(c)In the rest of this section, the varied initial droplet size and velocity of different flow turbulence are focused on. And it is shown that with the same nozzle diameter and flow mass rate specified, the turbulent (single-phase) liquid flow fits better with the experiment results than the laminar (flipped) liquid flow.


### 2.4.1. Determination of the Initial Droplet Size Distribution

According to our investigation, the initial droplet diameter distribution is closely related to the nozzle state. To indicate the connection, a two-parameter RR (Rosin-Rammler) distribution is used to represent the droplet diameter distribution as below, characterized by the most probable droplet sized0 and a spread parameter n. The mass fraction of droplets with diameters greater than d is calculated by(33)Φd=e−d/d0n.The first parameter required to specify the droplet size distribution is the most probable droplet sized0. For a single-phase nozzle flow, the correlation of Wu et al. (Wu, Tseng, and Faeth, 1992) [24] is applied to calculate SMD d32 considering the initial drop size related to the turbulence quantities of the liquid jet. Snyder [25] gives the most general relationship between SMD and most probable diameter for a Rosin-Rammler distribution.(34)d32=133.0λWe−e0.74,(35)We=ρlu2λσ,(36)d0=1.2726d321−121/n,where λ=d/8,λ is the radial integral length scale at the jet exit based upon fully developed turbulent pipe flow, and We is the Weber number of the liquid jet.For a cavitating nozzle flow, the correlation of Wu can still be applied, and yet the length scale for a cavitating nozzle isλ=deff/8, where deff is the effective diameter of the exiting liquid jet according to Schmidt and Corradini [26].For the case of a flipped nozzle flow, the initial droplet diameter is set to the diameter of the liquid jet(37)d0=dCct,where Cct is a theoretical constant equal to 0.611, which comes from potential flow analysis of flipped nozzles.The second parameter required to specify the droplet size distribution is the spread parametern. The values for the spread parameter are determined from past modeling experience and experimental observations. Table 3 lists the values of n for three flow states. The larger the value of the spread parameter, the narrower the droplet size distribution.Table 3
Theoretical values of spread parametern for different liquid states.
StateSpread parameterSingle-phase3.5Cavitating1.5Flipped∞Having specified the most probable diameter and the spread parameter, the initial droplet size distribution can thus be determined. It should be noted that the actual size distribution may be a little different with the theoretical one for the limited number of droplets injected from the nozzle exit. Figure4 shows the initial droplet size distribution of three nozzle flow states. It can be seen that the size of injected droplets of turbulent internal flow tends to be smaller and more uniform while the laminar internal flow tends to be larger and more centralized. It has been experimentally figured out that the droplet size distribution in the near nozzle region has a tremendous effect on downstream droplet properties.Figure 4
Initial droplet size distribution of different flow states.Figure5 illustrates the SMD distribution at x=50,100mm in the central plane distinguished by the flow turbulence with different initial droplet diameter distributions. Therefore, the turbulent liquid jet fits better with the experimental results in the aspect of the droplet size distribution.Figure 5
Downstream SMD distributions with different flow state: (a) turbulent and (b) laminar.
(a)(b)


### 2.4.2. Determination of Initial Droplet Velocity

For a single-phase nozzle, the estimate of exit velocityu comes from the conservation of mass and the assumption of a uniform exit velocity:(38)u=meffρlA.For a cavitating nozzle, an expression for a higher velocity over a reduced area derived by Schmidt and Corradini [26] is presented here instead of a uniform exit velocity:(39)u=2Ccp1−p2+1−2CcpvCc2ρlρ1−pv,where Cc is the contraction coefficient by Nurick’s [27] fit, p1 is the internal upstream pressure of the nozzle, p2 is the internal downstream pressure of the nozzle, and pv is the vapor pressure of the nozzle.For a flipped nozzle, the exit velocity is derived from the conservation of mass and the value of the reduced flow area:(40)u=meffρlCctA.The tremendous effect on the initial droplet velocity directly affects the penetration height of the liquid spray. According to our investigation, the penetration height of a turbulent liquid jet fits better with experimental results [6] than that of a laminar liquid jet in the same gaseous conditions as shown in Figure 6. It can be explained that the laminar liquid jet has a larger initial droplet velocity with the mass flow rate specified due to the flow contraction in the nozzle, which is inconsistent with the actual situation. Figure 7 indicates the higher velocity of laminar flow compared with turbulent flow.Figure 6
Comparison of spray penetration height with experimental results in the central plane: (a) turbulent and (b) laminar.
(a)(b)Figure 7
Calculation domain and meshing.


## 2.4.1. Determination of the Initial Droplet Size Distribution

According to our investigation, the initial droplet diameter distribution is closely related to the nozzle state. To indicate the connection, a two-parameter RR (Rosin-Rammler) distribution is used to represent the droplet diameter distribution as below, characterized by the most probable droplet sized0 and a spread parameter n. The mass fraction of droplets with diameters greater than d is calculated by(33)Φd=e−d/d0n.The first parameter required to specify the droplet size distribution is the most probable droplet sized0. For a single-phase nozzle flow, the correlation of Wu et al. (Wu, Tseng, and Faeth, 1992) [24] is applied to calculate SMD d32 considering the initial drop size related to the turbulence quantities of the liquid jet. Snyder [25] gives the most general relationship between SMD and most probable diameter for a Rosin-Rammler distribution.(34)d32=133.0λWe−e0.74,(35)We=ρlu2λσ,(36)d0=1.2726d321−121/n,where λ=d/8,λ is the radial integral length scale at the jet exit based upon fully developed turbulent pipe flow, and We is the Weber number of the liquid jet.For a cavitating nozzle flow, the correlation of Wu can still be applied, and yet the length scale for a cavitating nozzle isλ=deff/8, where deff is the effective diameter of the exiting liquid jet according to Schmidt and Corradini [26].For the case of a flipped nozzle flow, the initial droplet diameter is set to the diameter of the liquid jet(37)d0=dCct,where Cct is a theoretical constant equal to 0.611, which comes from potential flow analysis of flipped nozzles.The second parameter required to specify the droplet size distribution is the spread parametern. The values for the spread parameter are determined from past modeling experience and experimental observations. Table 3 lists the values of n for three flow states. The larger the value of the spread parameter, the narrower the droplet size distribution.Table 3
Theoretical values of spread parametern for different liquid states.
StateSpread parameterSingle-phase3.5Cavitating1.5Flipped∞Having specified the most probable diameter and the spread parameter, the initial droplet size distribution can thus be determined. It should be noted that the actual size distribution may be a little different with the theoretical one for the limited number of droplets injected from the nozzle exit. Figure4 shows the initial droplet size distribution of three nozzle flow states. It can be seen that the size of injected droplets of turbulent internal flow tends to be smaller and more uniform while the laminar internal flow tends to be larger and more centralized. It has been experimentally figured out that the droplet size distribution in the near nozzle region has a tremendous effect on downstream droplet properties.Figure 4
Initial droplet size distribution of different flow states.Figure5 illustrates the SMD distribution at x=50,100mm in the central plane distinguished by the flow turbulence with different initial droplet diameter distributions. Therefore, the turbulent liquid jet fits better with the experimental results in the aspect of the droplet size distribution.Figure 5
Downstream SMD distributions with different flow state: (a) turbulent and (b) laminar.
(a)(b)


## 2.4.2. Determination of Initial Droplet Velocity

For a single-phase nozzle, the estimate of exit velocityu comes from the conservation of mass and the assumption of a uniform exit velocity:(38)u=meffρlA.For a cavitating nozzle, an expression for a higher velocity over a reduced area derived by Schmidt and Corradini [26] is presented here instead of a uniform exit velocity:(39)u=2Ccp1−p2+1−2CcpvCc2ρlρ1−pv,where Cc is the contraction coefficient by Nurick’s [27] fit, p1 is the internal upstream pressure of the nozzle, p2 is the internal downstream pressure of the nozzle, and pv is the vapor pressure of the nozzle.For a flipped nozzle, the exit velocity is derived from the conservation of mass and the value of the reduced flow area:(40)u=meffρlCctA.The tremendous effect on the initial droplet velocity directly affects the penetration height of the liquid spray. According to our investigation, the penetration height of a turbulent liquid jet fits better with experimental results [6] than that of a laminar liquid jet in the same gaseous conditions as shown in Figure 6. It can be explained that the laminar liquid jet has a larger initial droplet velocity with the mass flow rate specified due to the flow contraction in the nozzle, which is inconsistent with the actual situation. Figure 7 indicates the higher velocity of laminar flow compared with turbulent flow.Figure 6
Comparison of spray penetration height with experimental results in the central plane: (a) turbulent and (b) laminar.
(a)(b)Figure 7
Calculation domain and meshing.


## 3. Computational Conditions

Lin et al. (Lin and Kennedy, 2002) [6] carried out the test research of a water liquid jet injected into supersonic crossflow with a Mach number of 1.94 and acquired abundant experimental data using a two-component phase Doppler particle analyzer (PDPA), which has been commonly simulated to verify new numerical models. In this section, the exact case mentioned above was simulated by FLUENT software to validate our physical models.Considering the computational cost, the calculation domain was set to be a rectangular region near the injector withLx×Ly×Lz=200mm×40mm×40mm. The calculation domain was meshed by structured grids with a total number of 761904 cells as shown in Figure 7. The position of the injector was designed to be x0,y0,z0=50mm,0,0, the vicinity of where was locally refined as shown in Figure 8. The k−ω−SST turbulence model is applied to calculate the turbulence of supersonic gas due to its high prediction ability in the far-field and near-wall regions. Default values of the coefficients of this model embedded in FLUENT software are applied.Figure 8
Vertical view of encrypted mesh near the nozzle.Water was used as the simulated liquid in the research which has a density of 998 kg/m3, viscosity of 2.67×10−3 kg/(m‧s), and surface tension of 0.072 N/m. The momentum flux ratio q=ρlvl2/ρ∞v∞2 is set to be constant 7, and other detailed parameters of liquid and gas are given in Tables 4 and 5.Table 4
Air inlet conditions.
SpeciesTotal pressure/kPaTotal temperature/KFree stream velocity/m‧s-1Air79%N2+21%O2206533678.13Table 5
Water injection conditions.
SpeciesTemperature/KFlow rate/(g/s)Velocity(m/s)Liquid water298.156.41532.73


## 4. Results and Discussion

### 4.1. Model Validation

#### 4.1.1. Spatial Distribution of the Liquid Spray

Spray penetration height is an important atomization characteristic to indicate the liquid-gas mixing effect. Figure9 illustrates the numerical result compared with the experimental result. The black dots represent the averaged spray droplets of 100 instantaneous moments, and the red dashed line represents the experimental correlation function developed by Lin et al. (Lin and Kennedy, 2002) [6].Figure 9
Model verification in spray penetration height.Cross-sectional distribution of the spray is another important parameter to evaluate the mixing characteristic of liquid and gas phases. Figure10 shows the spray spread at the position of x=50mm compared with the experimental result. The black dots represent the averaged spray droplets of 100 instantaneous moments, and the red dashed line is the experimental correlation function as seen in Figure 11 developed by Wu [24].Figure 10
Model verification in spray cross-sectional spread.Figure 11
Schematic diagram of Wu’s empirical formulas.It is worth noting that although Lin also showed their experimental results of the spray cross-sectional distribution, however, the explicit empirical correlations were not summarized in their research as done by Wu. Therefore, in this part, our numerical result is compared with Wu’s experimental correlation function. As can be seen, the spray foot is not observed as expected, because the thin liquid spray in the near-wall region calculated has much weaker obstruction than the experiment. In the Eulerian-Lagrangian methods, it is assumed that the liquid phase is sufficiently dilute that droplet-droplet interactions, and the effects of the droplet volume fraction on the gas phase are negligible. The result can be better improved by increasing the droplet density of the liquid spray, which is beneficial to enhance the entrainment effect of gas-phase vortices on the droplets but increase the computational cost.The consistency with the experimental results of the spatial distributions of the liquid spray firmly proves that present physical models are reliable.


#### 4.1.2. Centerline Profiles of Droplet Properties

Centerline distribution profiles of droplet properties in the free stream direction can be normalized by local penetration height to obtain universal curves. Figure12 illustrates the normalized SMD distribution along the y-axis at different x positions in the central plane. It shows that the numerical calculation can basically reproduce the S-shaped profile observed in the experiment even though the increasing trend of SMD with the y-axis height in the upper periphery of the liquid spray is not successfully simulated, which can be attributed to the weaker liquid-gas interaction than reality. Apart from the experiment results, our results of modified models is also compared with Li’s simulation results [13] which were obtained under the same numerical conditions but refined calculation region of 408×201×201 grid points in the x, y, and z directions. And the comparison shows considerable consistency in the SMD, absolute, and relative velocities distributions at the position of x=50mm as shown in Figure 13. All these comparisons validate the reliability of our modified models.Figure 12
Model verification in the SMD distribution at different locations in the central plane.Figure 13
Comparisons of SMD, absolute, and relative velocities distributions with the simulation results of Li [13] at x=50mm.To further investigate the difference between numerical and experimental results, the filled contours of the SMD distribution at the position ofx=50mm are contrasted in Figure 14. It can be found that the simulated gas-phase vortex field is remarkably weaker than the actual vortex field because the droplet volumes are neglected in the Eulerian-Lagrangian scheme. As a result, in simulations, it is the aerodynamic shear force that simply controls the downstream SMD distribution. That is to say that the closer to the free airflow, the smaller the size of droplets is. However, the experimental results show that the maximum droplet size appears at both the top of the spray plume and the y/h=0.3 location due to the complicated liquid-gas interaction and complex gas-phase vortex field.Figure 14
Comparison of SMD contour with the experimental result atx=50mm: (a) numerical result and (b) experimental result.
(a)(b)


### 4.2. Grid Independence Verification

To investigate the reliability of current grids, a higher-resolution set of grid with a total number of 2093184 cells is used to calculate this problem. The top view of this refined set of grid is screenshot as Figure15.Figure 15
Vertical viewing of the refined mesh.The results of the penetration, spread distribution and downstream properties curves are compared between two sets of grid as shown in Figures16–21. It can be naturally concluded that despite the trivial differences, the grid independence is got basically verified.Figure 16
Spray penetration att=0.5ms with high grid resolution.Figure 17
Spray penetration att=0.5ms with low grid resolution.Figure 18
Spray spread att=0.5ms with high grid resolution.Figure 19
Spray spread att=0.5ms with low grid resolution.Figure 20
SMD comparison atx=50mm between two sets of grids.Figure 21
Velocity of droplet comparison atx=50mm between two sets of grids.


### 4.3. Modeling of the SMD Distribution

Due to actual demands, the scale effects of droplet properties need to be explored by several simulated cases with different airflow Mach numbers and nozzle diameters. Although the S-shaped profile of the SMD distribution can inevitably increase the modeling difficulty, a new modeling approach is presented in this section to help settle the SMD quantification problem. To solve that, modified power functions are established which show good accuracy in predicting the SMD distribution. Dozens of numerical cases are launched to acquire statistically enough data for mathematical analysis with the aid of the multivariate nonlinear fitting method. Detailed operational conditions are listed below, and the expressions of the SMD distribution in the free stream direction as well as the errors are exhibited.


#### 4.3.1. Operational Conditions

In this section, a total of 54 sets of operating conditions are designed, among which the nozzle diameterd varies in 0.3, 0.5, 1, 2, and 3 mm, and the airflow Mach number Ma varies in 2.1, 3, and 4. It is worth noting that since the flux momentum ratio q and the airflow Weber number We have a significant effect on droplet properties, when considering the scale effects of the nozzle diameter and the airflow Mach number, the injecting velocity of the liquid jet and the density of the compressible air are controlled accordingly to keep q and We always constant. Detailed settings of parameters are shown in Table 6.Table 6
Operational conditions with different nozzle diameters and airflow Ma.
No.d (mm)meff· (kg/s)x/dGas crossflowMa=1.94T0=300KP0=29kPa1-70.30.0011910~708-140.50.0047510~7015-2110.0190010~7022-2820.0760010~7029-3330.1710010~50No.MaP0 (pa)x/dWater jetρ=998kg/m3T=300K33-402.197.4310~7040-47347.7410~7047-54426.8510~70


#### 4.3.2. Modeling Instruments

In this section, the traditional power function for developing penetration height formulas(41)y=Ax1ax2bis modified into a new-form power function to describe the SMD distribution in the free stream direction(42)SMD/δ=Ax/dad/δbMa+B,where δ is the boundary layer thickness, and d is the nozzle diameter, Ma is the free stream Mach number, x is the distance from the nozzle. A, B, a, b, and c are parameters to be determined.These unknown parameters are optimized by using the multivariate nonlinear fitting approach to correlate the local SMD at differentx positions with various nozzle diameters and airflow Mach numbers. Eventually, the formulas of SMD are given separately at different y positions in Table 7, and the MSE (mean squared error) represents the fitting error. As can be seen, the fitting results are acceptably accurate within a certain error range. Figure 22 is the schematic diagram of the fitting results.Table 7
Distribution of SMD and the corresponding MSE at the position ofy/d=3,5,and7.
Expressions of the distribution of SMDMSE (mm)y/d=3SMD/δ=0.0023x/d−0.09d/δ1.27Ma−0.98+0.03920.53y/d=5SMD/δ=0.0050x/d−0.09d/δ1.04Ma−1.01+0.03870.97y/d=7SMD/δ=0.0037x/d−0.09d/δ1.08Ma−1.01+0.03880.84Figure 22
Schematic diagram of multivariate nonlinear fitting.


#### 4.3.3. Result Analysis

The whole expression of SMD can be analyzed from two components. For a fixedx position, the former component is controlled by the nozzle diameter and airflow Mach number while the latter component remains constant. Since gaseous Weber number We and flux momentum ratio q remain unchanged in different settings of operational conditions, the latter component should be regarded as inherent influence of We and q. Mathematically, it is quite reasonable because it was proven that We has a significant effect on the droplet diameter. About the former component, it can be observed that SMD is roughly proportionate with d and negatively proportionate with Ma. Therefore, it can be concluded that apart from the component affected by the changeless parameters, the changeable component of the value of SMD is closely related to the nozzle diameter and airflow Mach number.


## 4.1. Model Validation

### 4.1.1. Spatial Distribution of the Liquid Spray

Spray penetration height is an important atomization characteristic to indicate the liquid-gas mixing effect. Figure9 illustrates the numerical result compared with the experimental result. The black dots represent the averaged spray droplets of 100 instantaneous moments, and the red dashed line represents the experimental correlation function developed by Lin et al. (Lin and Kennedy, 2002) [6].Figure 9
Model verification in spray penetration height.Cross-sectional distribution of the spray is another important parameter to evaluate the mixing characteristic of liquid and gas phases. Figure10 shows the spray spread at the position of x=50mm compared with the experimental result. The black dots represent the averaged spray droplets of 100 instantaneous moments, and the red dashed line is the experimental correlation function as seen in Figure 11 developed by Wu [24].Figure 10
Model verification in spray cross-sectional spread.Figure 11
Schematic diagram of Wu’s empirical formulas.It is worth noting that although Lin also showed their experimental results of the spray cross-sectional distribution, however, the explicit empirical correlations were not summarized in their research as done by Wu. Therefore, in this part, our numerical result is compared with Wu’s experimental correlation function. As can be seen, the spray foot is not observed as expected, because the thin liquid spray in the near-wall region calculated has much weaker obstruction than the experiment. In the Eulerian-Lagrangian methods, it is assumed that the liquid phase is sufficiently dilute that droplet-droplet interactions, and the effects of the droplet volume fraction on the gas phase are negligible. The result can be better improved by increasing the droplet density of the liquid spray, which is beneficial to enhance the entrainment effect of gas-phase vortices on the droplets but increase the computational cost.The consistency with the experimental results of the spatial distributions of the liquid spray firmly proves that present physical models are reliable.


### 4.1.2. Centerline Profiles of Droplet Properties

Centerline distribution profiles of droplet properties in the free stream direction can be normalized by local penetration height to obtain universal curves. Figure12 illustrates the normalized SMD distribution along the y-axis at different x positions in the central plane. It shows that the numerical calculation can basically reproduce the S-shaped profile observed in the experiment even though the increasing trend of SMD with the y-axis height in the upper periphery of the liquid spray is not successfully simulated, which can be attributed to the weaker liquid-gas interaction than reality. Apart from the experiment results, our results of modified models is also compared with Li’s simulation results [13] which were obtained under the same numerical conditions but refined calculation region of 408×201×201 grid points in the x, y, and z directions. And the comparison shows considerable consistency in the SMD, absolute, and relative velocities distributions at the position of x=50mm as shown in Figure 13. All these comparisons validate the reliability of our modified models.Figure 12
Model verification in the SMD distribution at different locations in the central plane.Figure 13
Comparisons of SMD, absolute, and relative velocities distributions with the simulation results of Li [13] at x=50mm.To further investigate the difference between numerical and experimental results, the filled contours of the SMD distribution at the position ofx=50mm are contrasted in Figure 14. It can be found that the simulated gas-phase vortex field is remarkably weaker than the actual vortex field because the droplet volumes are neglected in the Eulerian-Lagrangian scheme. As a result, in simulations, it is the aerodynamic shear force that simply controls the downstream SMD distribution. That is to say that the closer to the free airflow, the smaller the size of droplets is. However, the experimental results show that the maximum droplet size appears at both the top of the spray plume and the y/h=0.3 location due to the complicated liquid-gas interaction and complex gas-phase vortex field.Figure 14
Comparison of SMD contour with the experimental result atx=50mm: (a) numerical result and (b) experimental result.
(a)(b)


## 4.1.1. Spatial Distribution of the Liquid Spray

Spray penetration height is an important atomization characteristic to indicate the liquid-gas mixing effect. Figure9 illustrates the numerical result compared with the experimental result. The black dots represent the averaged spray droplets of 100 instantaneous moments, and the red dashed line represents the experimental correlation function developed by Lin et al. (Lin and Kennedy, 2002) [6].Figure 9
Model verification in spray penetration height.Cross-sectional distribution of the spray is another important parameter to evaluate the mixing characteristic of liquid and gas phases. Figure10 shows the spray spread at the position of x=50mm compared with the experimental result. The black dots represent the averaged spray droplets of 100 instantaneous moments, and the red dashed line is the experimental correlation function as seen in Figure 11 developed by Wu [24].Figure 10
Model verification in spray cross-sectional spread.Figure 11
Schematic diagram of Wu’s empirical formulas.It is worth noting that although Lin also showed their experimental results of the spray cross-sectional distribution, however, the explicit empirical correlations were not summarized in their research as done by Wu. Therefore, in this part, our numerical result is compared with Wu’s experimental correlation function. As can be seen, the spray foot is not observed as expected, because the thin liquid spray in the near-wall region calculated has much weaker obstruction than the experiment. In the Eulerian-Lagrangian methods, it is assumed that the liquid phase is sufficiently dilute that droplet-droplet interactions, and the effects of the droplet volume fraction on the gas phase are negligible. The result can be better improved by increasing the droplet density of the liquid spray, which is beneficial to enhance the entrainment effect of gas-phase vortices on the droplets but increase the computational cost.The consistency with the experimental results of the spatial distributions of the liquid spray firmly proves that present physical models are reliable.


## 4.1.2. Centerline Profiles of Droplet Properties

Centerline distribution profiles of droplet properties in the free stream direction can be normalized by local penetration height to obtain universal curves. Figure12 illustrates the normalized SMD distribution along the y-axis at different x positions in the central plane. It shows that the numerical calculation can basically reproduce the S-shaped profile observed in the experiment even though the increasing trend of SMD with the y-axis height in the upper periphery of the liquid spray is not successfully simulated, which can be attributed to the weaker liquid-gas interaction than reality. Apart from the experiment results, our results of modified models is also compared with Li’s simulation results [13] which were obtained under the same numerical conditions but refined calculation region of 408×201×201 grid points in the x, y, and z directions. And the comparison shows considerable consistency in the SMD, absolute, and relative velocities distributions at the position of x=50mm as shown in Figure 13. All these comparisons validate the reliability of our modified models.Figure 12
Model verification in the SMD distribution at different locations in the central plane.Figure 13
Comparisons of SMD, absolute, and relative velocities distributions with the simulation results of Li [13] at x=50mm.To further investigate the difference between numerical and experimental results, the filled contours of the SMD distribution at the position ofx=50mm are contrasted in Figure 14. It can be found that the simulated gas-phase vortex field is remarkably weaker than the actual vortex field because the droplet volumes are neglected in the Eulerian-Lagrangian scheme. As a result, in simulations, it is the aerodynamic shear force that simply controls the downstream SMD distribution. That is to say that the closer to the free airflow, the smaller the size of droplets is. However, the experimental results show that the maximum droplet size appears at both the top of the spray plume and the y/h=0.3 location due to the complicated liquid-gas interaction and complex gas-phase vortex field.Figure 14
Comparison of SMD contour with the experimental result atx=50mm: (a) numerical result and (b) experimental result.
(a)(b)


## 4.2. Grid Independence Verification

To investigate the reliability of current grids, a higher-resolution set of grid with a total number of 2093184 cells is used to calculate this problem. The top view of this refined set of grid is screenshot as Figure15.Figure 15
Vertical viewing of the refined mesh.The results of the penetration, spread distribution and downstream properties curves are compared between two sets of grid as shown in Figures16–21. It can be naturally concluded that despite the trivial differences, the grid independence is got basically verified.Figure 16
Spray penetration att=0.5ms with high grid resolution.Figure 17
Spray penetration att=0.5ms with low grid resolution.Figure 18
Spray spread att=0.5ms with high grid resolution.Figure 19
Spray spread att=0.5ms with low grid resolution.Figure 20
SMD comparison atx=50mm between two sets of grids.Figure 21
Velocity of droplet comparison atx=50mm between two sets of grids.


## 4.3. Modeling of the SMD Distribution

Due to actual demands, the scale effects of droplet properties need to be explored by several simulated cases with different airflow Mach numbers and nozzle diameters. Although the S-shaped profile of the SMD distribution can inevitably increase the modeling difficulty, a new modeling approach is presented in this section to help settle the SMD quantification problem. To solve that, modified power functions are established which show good accuracy in predicting the SMD distribution. Dozens of numerical cases are launched to acquire statistically enough data for mathematical analysis with the aid of the multivariate nonlinear fitting method. Detailed operational conditions are listed below, and the expressions of the SMD distribution in the free stream direction as well as the errors are exhibited.


### 4.3.1. Operational Conditions

In this section, a total of 54 sets of operating conditions are designed, among which the nozzle diameterd varies in 0.3, 0.5, 1, 2, and 3 mm, and the airflow Mach number Ma varies in 2.1, 3, and 4. It is worth noting that since the flux momentum ratio q and the airflow Weber number We have a significant effect on droplet properties, when considering the scale effects of the nozzle diameter and the airflow Mach number, the injecting velocity of the liquid jet and the density of the compressible air are controlled accordingly to keep q and We always constant. Detailed settings of parameters are shown in Table 6.Table 6
Operational conditions with different nozzle diameters and airflow Ma.
No.d (mm)meff· (kg/s)x/dGas crossflowMa=1.94T0=300KP0=29kPa1-70.30.0011910~708-140.50.0047510~7015-2110.0190010~7022-2820.0760010~7029-3330.1710010~50No.MaP0 (pa)x/dWater jetρ=998kg/m3T=300K33-402.197.4310~7040-47347.7410~7047-54426.8510~70


### 4.3.2. Modeling Instruments

In this section, the traditional power function for developing penetration height formulas(41)y=Ax1ax2bis modified into a new-form power function to describe the SMD distribution in the free stream direction(42)SMD/δ=Ax/dad/δbMa+B,where δ is the boundary layer thickness, and d is the nozzle diameter, Ma is the free stream Mach number, x is the distance from the nozzle. A, B, a, b, and c are parameters to be determined.These unknown parameters are optimized by using the multivariate nonlinear fitting approach to correlate the local SMD at differentx positions with various nozzle diameters and airflow Mach numbers. Eventually, the formulas of SMD are given separately at different y positions in Table 7, and the MSE (mean squared error) represents the fitting error. As can be seen, the fitting results are acceptably accurate within a certain error range. Figure 22 is the schematic diagram of the fitting results.Table 7
Distribution of SMD and the corresponding MSE at the position ofy/d=3,5,and7.
Expressions of the distribution of SMDMSE (mm)y/d=3SMD/δ=0.0023x/d−0.09d/δ1.27Ma−0.98+0.03920.53y/d=5SMD/δ=0.0050x/d−0.09d/δ1.04Ma−1.01+0.03870.97y/d=7SMD/δ=0.0037x/d−0.09d/δ1.08Ma−1.01+0.03880.84Figure 22
Schematic diagram of multivariate nonlinear fitting.


### 4.3.3. Result Analysis

The whole expression of SMD can be analyzed from two components. For a fixedx position, the former component is controlled by the nozzle diameter and airflow Mach number while the latter component remains constant. Since gaseous Weber number We and flux momentum ratio q remain unchanged in different settings of operational conditions, the latter component should be regarded as inherent influence of We and q. Mathematically, it is quite reasonable because it was proven that We has a significant effect on the droplet diameter. About the former component, it can be observed that SMD is roughly proportionate with d and negatively proportionate with Ma. Therefore, it can be concluded that apart from the component affected by the changeless parameters, the changeable component of the value of SMD is closely related to the nozzle diameter and airflow Mach number.


## 4.3.1. Operational Conditions

In this section, a total of 54 sets of operating conditions are designed, among which the nozzle diameterd varies in 0.3, 0.5, 1, 2, and 3 mm, and the airflow Mach number Ma varies in 2.1, 3, and 4. It is worth noting that since the flux momentum ratio q and the airflow Weber number We have a significant effect on droplet properties, when considering the scale effects of the nozzle diameter and the airflow Mach number, the injecting velocity of the liquid jet and the density of the compressible air are controlled accordingly to keep q and We always constant. Detailed settings of parameters are shown in Table 6.Table 6
Operational conditions with different nozzle diameters and airflow Ma.
No.d (mm)meff· (kg/s)x/dGas crossflowMa=1.94T0=300KP0=29kPa1-70.30.0011910~708-140.50.0047510~7015-2110.0190010~7022-2820.0760010~7029-3330.1710010~50No.MaP0 (pa)x/dWater jetρ=998kg/m3T=300K33-402.197.4310~7040-47347.7410~7047-54426.8510~70


## 4.3.2. Modeling Instruments

In this section, the traditional power function for developing penetration height formulas(41)y=Ax1ax2bis modified into a new-form power function to describe the SMD distribution in the free stream direction(42)SMD/δ=Ax/dad/δbMa+B,where δ is the boundary layer thickness, and d is the nozzle diameter, Ma is the free stream Mach number, x is the distance from the nozzle. A, B, a, b, and c are parameters to be determined.These unknown parameters are optimized by using the multivariate nonlinear fitting approach to correlate the local SMD at differentx positions with various nozzle diameters and airflow Mach numbers. Eventually, the formulas of SMD are given separately at different y positions in Table 7, and the MSE (mean squared error) represents the fitting error. As can be seen, the fitting results are acceptably accurate within a certain error range. Figure 22 is the schematic diagram of the fitting results.Table 7
Distribution of SMD and the corresponding MSE at the position ofy/d=3,5,and7.
Expressions of the distribution of SMDMSE (mm)y/d=3SMD/δ=0.0023x/d−0.09d/δ1.27Ma−0.98+0.03920.53y/d=5SMD/δ=0.0050x/d−0.09d/δ1.04Ma−1.01+0.03870.97y/d=7SMD/δ=0.0037x/d−0.09d/δ1.08Ma−1.01+0.03880.84Figure 22
Schematic diagram of multivariate nonlinear fitting.


## 4.3.3. Result Analysis

The whole expression of SMD can be analyzed from two components. For a fixedx position, the former component is controlled by the nozzle diameter and airflow Mach number while the latter component remains constant. Since gaseous Weber number We and flux momentum ratio q remain unchanged in different settings of operational conditions, the latter component should be regarded as inherent influence of We and q. Mathematically, it is quite reasonable because it was proven that We has a significant effect on the droplet diameter. About the former component, it can be observed that SMD is roughly proportionate with d and negatively proportionate with Ma. Therefore, it can be concluded that apart from the component affected by the changeless parameters, the changeable component of the value of SMD is closely related to the nozzle diameter and airflow Mach number.


## 5. Conclusions

This research mainly focuses on the atomization process of a liquid jet injected into supersonic crossflows. The discrete phase model and modified KH-RT breakup model are adopted to investigate the inflow turbulence of a liquid jet which can significantly affect the distribution profiles of downstream droplet properties. Also, a new method to describe the SMD distribution is presented as our main innovative point. Abundant numerical data is acquired, and main conclusions are summarized below.(1)
The inflow turbulence of a liquid jet has significant influence on the primary breakup, which is reflected on the droplet initial size and velocity. And a turbulent jet is found to fit better with experimental results than a laminar one(2)
In the supersonic environment, it was discovered that the KH breakup has more control on the whole breakup process than RT breakup. The KH breakup size constantB0 can significantly affect the downstream SMD distribution profile(3)
The distribution of SMD can be well described with a modified power function presented in this paper. The mathematical results implied that gaseous Weber numberWe and flux momentum ratio q still have considerable effects on the final droplet size in the supersonic airflow, while the scale effect of nozzle diameter and airflow Mach number can also be significant. And their effects on the SMD distribution are fairly comparable(43)z/d2a2+y/d−m2cy/d−b+b−m2=1,z/d=ky/d−ys/d2+zs/d,axd,q,y=1.33q0.21xd0.33−0.76q0.33xd0.23γ−0.5γ−0.5γ−0.51.1,bxd,q,y=3.34q0.36xd0.25−0.83q0.2xd0.32γ−0.5γ−0.5γ−0.51.1,mxd,q,y=1.93q0.48xd0.15−0.24q0.25xd0.35γ−0.5γ−0.5γ−0.51.1,c=0.25xd−1.7q−510−2.



---

*Source: 8828015-2021-07-09.xml*