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

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 the Ma = 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.


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 investi-gate 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][8][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 drop-lets. 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  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 Section 2, 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.

Mathematical Models
International Journal of Aerospace Engineering 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, C d is the drag coefficient, u ! p is the droplet velocity, ρ p is the droplet density, d p is the droplet diameter, F other is other forces per unit gas mass, τ r is the droplet relaxation time [18], and R e is the relative Reynolds number .

Liquid-Phase Dynamic Equation.
d u 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 Integration of time in Equation (4) yields the velocity of the droplet at each point along the trajectory, with the trajectory itself predicted by The new location x n+1 p can be computed from where a includes accelerations due to all other forces except drag force, and u n p and n n represent droplet velocity and gas velocity at the old location.

Turbulence Model.
In this paper, the k − ω − 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 averagef (called Farve average) is usually preferred ðf = ρf / ρÞ. Any quantity f may be splinted into mean and fluctuating components as The average balance equations of continuous phase are written as  ðc v ′ c p are specific heat capacity of constant volume and pressure, respectively. The turbulent kinetic energy k is written as where μ is the molecular viscosity, and K is thermal conductivity.
To provide closures for the unknown turbulent kinetic energy k, 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: The mixing function F 1 is defined by where y is the vertical distance from the wall surface, D kw is the positive part of the cross diffusion term, and D kw is determined by International Journal of Aerospace Engineering All the coefficients σ k , σ w , β, γ can be calculated in a uniform form by where ϕ 1 ′ ϕ 2 ′ is the corresponding coefficients of k − ε and k − ω turbulence models, respectively. The eddy viscosity coefficient v t is defined by where F 2 = tanh ðarg 2 2 Þ, Ω is the eddy vector.   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.

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 Figure 1. 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.
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]. where We g p . The radius of child droplet stripped from the cylindrical jet is τ KH is the KH breakup time calculated by where B 0 is the KH breakup radius constant, and B 1 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].
The RT breakup child droplet radius related to the wave length of the fastest growing wave is calculated by The RT breakup time is calculated as below where α d is the acceleration of the droplet, C 1 is the RT breakup radius constant, and C 2 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      Table 1.
On the basis of Yang's work, the optimal KH-RT breakup model constants were further investigated in this paper. Figure 2 illustrates that the downstream SMD distribution is mainly controlled by the KH radius constant B 0 rather than time constant B 1 . It is also found that the RT radius constant C 1 has much less influence on the downstream SMD than the KH radius constant B 0 . 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 C 1 is probably because the present liquid core length is larger than expected.
Therefore, using the idea of controlling variables, the optimal values of breakup constants B 0 , B 1 , C 1 , C 2 for supersonic simulation of two-phase interactions calculated by FLUENT software are shown in Table 2.

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 Figure 3, 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.
where p 1 , p 2 is the controllable upstream and downstream pressure, respectively.
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

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 size d 0 and a spread parameter n. The mass fraction of droplets with diameters greater than d is calculated by The first parameter required to specify the droplet size distribution is the most probable droplet size d 0 . For a single-phase nozzle flow, the correlation of Wu et al. (Wu, Tseng, and Faeth, 1992) [24] is applied to calculate SMD ðd 32 Þ 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.
We = ρ l u 2 λ σ , ð35Þ 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 λ = d eff /8, where d eff 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 where C ct 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 parameter n. 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.
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. Figure 4 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 5 illustrates the SMD distribution at x = 50,100 mm in the central plane distinguished by the flow turbulence 9 International Journal of Aerospace Engineering 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.

Determination of Initial Droplet Velocity.
For a single-phase nozzle, the estimate of exit velocity u comes from the conservation of mass and the assumption of a uniform exit velocity: 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: where C c is the contraction coefficient by Nurick's [27] fit, p 1 is the internal upstream pressure of the nozzle, p 2 is the internal downstream pressure of the nozzle, and p v 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: 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.

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 with L x × L y × L z = 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 ðx 0 , y 0 , z 0 Þ = ð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.
Water was used as the simulated liquid in the research which has a density of 998 kg/m 3 , viscosity of 2:67 × 10 −3 kg/(m‧s), and surface tension of 0.072 N/m. The momentum flux ratio q = ρ l v 2 l /ρ ∞ v 2 ∞ is set to be constant 7, and other detailed parameters of liquid and gas are given in Tables 4  and 5.

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. Figure 9 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].
Cross-sectional distribution of the spray is another important parameter to evaluate the mixing characteristic of liquid and gas phases. Figure 10 shows the spray spread at the position of x = 50 mm 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].
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.  Figure 12 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.
To further investigate the difference between numerical and experimental results, the filled contours of the SMD distribution at the position of x = 50 mm 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

12
International Journal of Aerospace Engineering 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.

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 prob-lem. The top view of this refined set of grid is screenshot as Figure 15.
The results of the penetration, spread distribution and downstream properties curves are compared between two sets of grid as shown in Figures 16-21. It can be naturally concluded that despite the trivial differences, the grid independence is got basically verified. 13 International Journal of Aerospace Engineering 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.

Operational Conditions.
In this section, a total of 54 sets of operating conditions are designed, among which the nozzle diameter d 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.
is modified into a new-form power function to describe the SMD distribution in the free stream direction 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 different x 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

Data Availability
The datasets generated during and/or analyzed during the current study are available from the corresponding author on reasonable request.

16
International Journal of Aerospace Engineering