Eulerian–Lagrangian method for liquid jet atomization in supersonic crossflow using statistical injection model

Abstract

An integrated Eulerian–Lagrangian method has been proposed to simulate liquid jet in supersonic flow. The carrier fluid surrounding the liquid droplets is described by multicomponent Navier–Stokes equations based on stationary (Eulerian) Cartesian grid, and the Lagrangian-discrete droplet method is employed to represent the behavior of atomized droplets in this work. A total of three classic breakup models, Taylor analogy breakup model, Reitz wave model, and Kelvin–Helmholtz/Rayleigh–Taylor hybrid model, are discussed under supersonic conditions, and model predictions are compared to experimental data through multiple perspectives quantitatively. More accurate predictions of liquid penetration, as well as droplet size distribution, can be achieved for specific conditions with Kelvin–Helmholtz/Rayleigh–Taylor hybrid model comparing to other breakup models. Additionally, a statistical injection model has been introduced to depict droplet dispersion in the near-nozzle region. The probability density function is utilized to predict the size and position of the injected droplet, and the four commonly used probability density functions, uniform, chi-squared, Nukiyama–Tanasawa, and Rosin–Rammler, are analyzed. The influence of statistical representation of injection condition is discussed. The simulation results indicate that the random components in the velocity and droplet size of the injected droplets have a significant impact on the structure of the liquid jet in high-speed flow.

Keywords

Numerical simulation Eulerian–Lagrangian method multiphase flow high-speed flow probability density function atomization

Introduction

A liquid jet flow injected into high-speed gaseous crossflows is one of the most widely used techniques for fuel–air atomization enhancement and mixing augmentation in propulsion systems such as the internal combustion engine, the liquid rocket engine, and the scramjet.¹ The understanding of procedures of fuel spray is a great challenge in liquid injection for practical applications. The procedure of these propulsion systems involves a number of complex, closely coupled physical and chemical processes, which include transient three-dimensional liquid atomization mixing with multicomponent gases and fuel combustion. The residence time of liquid fuel is millisecond ordered between the engine capture and exit through the nozzle in air-breathing vehicles, and therefore, the liquid breakup and fuel mixing fleeting become a key issue in high-speed crossflow.² Thus, an accurate prediction of liquid fuel spray process is necessary for better producing reliable engine performances.

Several studies on atomization of liquid jet have been presented in the literatures.^3–11 Experimental investigations of nonturbulent round liquid jets injected into gaseous crossflows conducted by Mazallon et al.¹² show that liquid surface deformation and breakup properties were not strongly affected by the jet-to-air momentum flux ratio at relatively large liquid/gas density ratios. Lin and his research group considered both pure-liquid and aerated-liquid jets in supersonic crossflows, and their experimental studies focused on information about velocity, liquid volume flux, and drop size distributions after breakup.^13–15 Early numerical studies of various liquid spray breakup models used to simulate liquid atomization were reviewed by Bravo and Kweon.¹⁶ Madabhushi¹⁷ developed a breakup model to simulate spray structures of water jets injected into subsonic crossflows, accounting for the effect of jet bending, aerodynamic forces, droplet distortion, and dispersion. Im et al.¹⁸ performed a Reynolds-averaged Navier–Stokes (RANS) study of a Lagrangian-based supersonic crossflow and proposed a modified breakup model to simulate the catastrophic regime in a supersonic flow based on an existing Kelvin–Helmholtz/Rayleigh–Taylor (K-H/R-T) hybrid model. Detailed numerical studies of turbulent liquid jets injected into a subsonic gaseous crossflow were conducted by Herrmann,¹⁹ indicating the breakup of the liquid jet via two main simultaneous atomization mechanisms. Further to density ratio, he found that the resulting drop sizes were a function of the chosen density ratio. Xiao et al.^1,20 investigated the primary breakup mechanism of liquid jet in an air crossflow using a two-phase flow large eddy simulation (LES) and found that Rayleigh–Taylor instability plays an important role in the primary breakup of liquid jet in a supersonic crossflow, which is consistent with the atomization mechanism for a liquid jet in a subsonic incompressible crossflow. As reviewed by Huang²¹ and Huang and Yan,²² several aspects, including the jet-to-crossflow pressure ratio, the geometric configuration of the injection port, the number of injection ports, and the injection angle, could affect the transverse injection flow field structure.

In this work, a marker-based Eulerian–Lagrangian method for solving the atomization of a water jet in supersonic crossflow is outlined. The carrier fluid surrounding the liquid droplets is described by the multicomponent Navier–Stokes equations, and the behavior of the atomized droplets is tracked by Lagrangian markers. Effects of pressure gradient caused by shock wave and turbulence on the spray droplets are considered in the proposed method. Statistical descriptions of initial droplets, through the joint of probability density function (PDF) of the droplet velocity and diameter, are introduced to make the liquid jets closer to real situations. The statistical descriptions of injection conditions, which emphasize the effect of droplet sizes, velocities, and locations, are the unique feature of this approach. A total of three types of classical breakup models, namely, Taylor analogy breakup (TAB) model,²³ Reitz wave model,²⁴ and K-H/R-T hybrid model,²⁵ are compared with experimental data available from Lin et al.¹⁴ for modeling validation as well as a case without breakup model to prove the necessity of using breakup models. A total of four types of PDFs, which are uniform, chi-squared, Nukiyama–Tanasawa, and Rosin–Rammler, respectively, are analyzed, and the influence on calculation results with different statistical representation of injection conditions is discussed. Numerical results are compared with those of experimental investigations by spray penetration height and droplet size distribution. The most reasonable injection model in this article is found out by comparing the numerical results.

Methodology

Gas-phase models

The effect of mass exchange and energy transportation between vapor and gas phase at the overall process of evaporation is significant. Thus, the calculation of the heat through latent heat contribution and the partial enthalpy of each component is based on real-gas mixture assumption, which is performed with multicomponent three-dimensional Navier–Stokes equations. The vector forms of the governing equations are given by⁹

\frac{\partial Q}{\partial t} + \frac{\partial E}{\partial x} + \frac{\partial F}{\partial y} + \frac{\partial G}{\partial z} = \frac{\partial E_{v}}{\partial x} + \frac{\partial F_{v}}{\partial y} + \frac{\partial G_{v}}{\partial z} + S_{l}

(1)

where the Q is the vector of conservative variables; E , F , and G are the vectors of convective fluxes; and E _v, F _v, and G _v are the vectors of viscous fluxes as follows²⁶

Q = {[\begin{matrix} ρ_{i} & ρ u & ρ v & ρ w & e_{t} \end{matrix}]}^{T}

(2)

E = {[\begin{matrix} ρ_{i} u & ρ u^{2} + p & ρ uv & ρ uw & (e_{t} + p) u \end{matrix}]}^{T}

(3)

F = {[\begin{matrix} ρ_{i} v & ρ vu & ρ v^{2} + p & ρ vw & (e_{t} + p) v \end{matrix}]}^{T}

(4)

G = {[\begin{matrix} ρ_{i} w & ρ wu & ρ wv & ρ w^{2} + p & (e_{t} + p) w \end{matrix}]}^{T}

(5)

E_{v} = {[\begin{matrix} ρ D_{im} \frac{\partial Y_{i}}{\partial x} & τ_{xx} & τ_{xy} & τ_{xz} & τ_{xx} u + τ_{xy} v + τ_{xz} w \end{matrix}]}^{T}

(6)

F_{v} = {[\begin{matrix} ρ D_{im} \frac{\partial Y_{i}}{\partial y} & τ_{yx} & τ_{yy} & τ_{yz} & τ_{yx} u + τ_{yy} v + τ_{yz} w \end{matrix}]}^{T}

(7)

G_{v} = {[\begin{matrix} ρ D_{im} \frac{\partial Y_{i}}{\partial z} & τ_{zx} & τ_{zy} & τ_{zz} & τ_{zx} u + τ_{zy} v + τ_{zz} w \end{matrix}]}^{T}

(8)

with

e_{t} = \sum_{i} ρ_{i} h_{i} - p + \frac{1}{2} ρ (u^{2} + v^{2} + w^{2})

(9)

\sum_{i = 1}^{N_{sp}} Y_{i} = 1

(10)

where e_t denotes the total energy per unit mass, h_i is the enthalpy of species i, Y_i denotes the mass fraction, D_im is the effective binary diffusivity of species i, $τ$ is the viscous stress tensor, and N_sp is the number of species. Together with Dalton’s law²⁶

p = \sum_{i = 1}^{Nsp} ρ_{i} \frac{R_{u}}{M W_{i}} T

(11)

where R_u is the universal gas constant and W_i is the molecular weight of species i.

In general, to determine an equilibrium temperature T of a gas mixture, N_sp+ 1 variables are needed where the explicit form of temperature cannot be found because the gas is in thermodynamical equilibrium, and all species possess the same temperature; thus, the temperature has to be calculated iteratively²⁷ where the Newton-iterative method is invoked.

Substituting equation (11) into equation (9)

\begin{matrix} F (T) = \sum_{i} ρ_{i} h_{i} - \sum_{i = 1}^{Nsp} ρ_{i} \frac{R_{u}}{M W_{i}} T \\ + \frac{1}{2} ρ (u^{2} + v^{2} + w^{2}) - e_{t} = 0 \end{matrix}

(12)

T^{(n + 1)} = T^{(n)} - \frac{F (T^{(n)})}{F' (T^{(n)})}

(13)

\begin{matrix} T^{(n + 1)} = \\ \frac{T^{(n)} \sum_{i} ρ_{i} C p_{i} (T^{(n)}) - \sum_{i} ρ_{i} h_{i} (T^{(n)}) - \frac{1}{2} ρ (u^{2} + v^{2} + w^{2}) + e_{t}}{\sum_{i} ρ_{i} C p_{i} (T^{(n)}) - \sum_{i = 1}^{Nsp} ρ_{i} \frac{R_{u}}{M W_{i}}} \end{matrix}

(14)

where the superscripts n and n+ 1 denote the previous time level and new time level, respectively.

The mass, momentum, and energy coupling between phases are implemented as a source term S _l in the governing equations

S_{l} = {[\begin{matrix} {\overset{\cdot}{ρ}}_{i} & F_{x} & F_{y} & F_{z} & \overset{\cdot}{Q} \end{matrix}]}^{T}

(15)

The production rate of species i due to evaporating spray ${\overset{\cdot}{ρ}}_{i}$ is given by

{\overset{\cdot}{ρ}}_{i} = - \frac{1}{Ω Δ t} \sum N_{n} ρ_{d} \frac{4}{3} π [{(r^{n + 1})}^{3} - {(r^{n})}^{3}]

(16)

The momentum source terms F_x, F_y, and F_z are given by

\begin{matrix} F_{x} = - \frac{1}{Ω Δ t} \sum N_{n} ρ_{d} \frac{4}{3} π [{(r^{n + 1})}^{3} u_{d}^{n + 1} - {(r^{n})}^{3} u_{d}^{n}] \\ F_{y} = - \frac{1}{Ω Δ t} \sum N_{n} ρ_{d} \frac{4}{3} π [{(r^{n + 1})}^{3} v_{d}^{n + 1} - {(r^{n})}^{3} v_{d}^{n}] \\ F_{z} = - \frac{1}{Ω Δ t} \sum N_{n} ρ_{d} \frac{4}{3} π [{(r^{n + 1})}^{3} w_{d}^{n + 1} - {(r^{n})}^{3} w_{d}^{n}] \end{matrix}

(17)

The work done by droplets to the gas phase $\overset{\cdot}{Q}$ is given by

\begin{matrix} \overset{\cdot}{Q} = - \frac{1}{Ω Δ t} \sum N_{n} ρ_{d} \frac{4}{3} π \\ {\begin{matrix} {(r^{n + 1})}^{3} e (T_{d}^{n + 1}) - {(r^{n})}^{3} e (T_{d}^{n}) \\ + {(r^{n + 1})}^{3} (u_{d}^{n + 1} - u_{d}^{n}) (u_{d}^{n + 1} - u^{n}) \\ + \frac{1}{2} [{(r^{n + 1})}^{3} - {(r^{n})}^{3}] {| u_{d}^{n + 1} - u_{d}^{n} |}^{2} \end{matrix}} \end{matrix}

(18)

where $Ω$ is the volume of a computational cell, $Δ t$ is the global physical time step, N_n is the number of droplets in the control volume, and r is the radius of droplets. The terms with the subscript d correspond to values of parcel, the superscript n is previous time level, and superscript n+ 1 is new time level.

Spray models

The Lagrangian-discrete droplet method is employed to represent the behavior of atomized droplets. This class of method tracks the spray droplets through Lagrangian markers moving on an Eulerian computational grid using a simple advection equation¹⁶

\frac{d x_{d}}{dt} = v_{d}

(19)

where subscript d denotes the computational parcel, which represents a number of droplets with the same size, position, temperature, and velocity.

The force acting on the droplet includes aerodynamic drag and gravitational force in general. While in high-speed flows, the pressure is rapidly increased within the shock front when considering a droplet pass through a shock wave, and the velocity of the droplet is accelerated due to the nonnegligible pressure gradients.²⁷ Therefore, the droplet acceleration is given by

\frac{d v_{d}}{dt} = \frac{3}{8} \frac{ρ}{ρ_{d}} \frac{| u + u' - v_{d} |}{r} (u + u' - v_{d}) C_{d} + g - \frac{\nabla p}{ρ_{d}}

(20)

where g is standard gravity and $\nabla p$ is the pressure gradient of the gas in the computational cell in which parcel d is located. The turbulence effect on the spray droplets is modeled by adding to a fluctuating velocity u′ to the gas velocity u as follows^28,29

u' = \sqrt{\frac{4 k}{3}} sign (2 R - 1) er f^{- 1} (| 2 R - 1 |)

(21)

where R is a random number ranged from 0 to 1, which satisfies a normal distribution, and k is the specific turbulent kinetic energy of the gas in the computational cell calculated using shear stress transport (SST) $k - ω$ turbulence model. The transport equations for k and $ω$ are given by²⁷

\frac{\partial ρ k}{\partial t} + \frac{\partial ρ u_{j} k}{\partial x_{j}} = - ρ τ_{ij} \frac{\partial u_{i}}{\partial x_{j}} - β^{*} ρ ω k + \frac{\partial}{\partial x_{j}} [(μ_{l} + σ_{k} μ_{t}) \frac{\partial k}{\partial x_{j}}]

(22)

\begin{matrix} \frac{\partial ρ ω}{\partial t} + \frac{\partial ρ u_{j} ω}{\partial x_{j}} = - \frac{α}{ν_{T}} ρ τ_{ij} \frac{\partial u_{i}}{\partial x_{j}} - β ρ ω^{2} \\ + \frac{\partial}{\partial x_{j}} [(μ_{l} + σ_{ω} μ_{t}) \frac{\partial ω}{\partial x_{j}}] \\ + 2 ρ (1 - F_{1}) σ_{ω 2} \frac{1}{ω} \frac{\partial k}{\partial x_{j}} \frac{\partial ω}{\partial x_{j}} \end{matrix}

(23)

For better describing the drag acting on the droplets in the supersonic regime, the high Mach number drag law is used which is similar to the spherical law with compressibility correction.^30,31 The C_d in equation (20) is the drag coefficient that can be expressed as

\begin{matrix} C_{d} = 2 + (C_{d 0} - 2) e^{- 3.07 \sqrt{γ (T_{c})} g ({Re}_{d}) M_{d} / M_{d} {Re}_{d} {Re}_{d}} \\ + \frac{h (M_{d})}{\sqrt{γ (T_{c})} M_{d}} e^{- {Re}_{d} / 2 M_{d}} \end{matrix}

(24)

where $γ (T_{c})$ is the specific heat ratio of the gas mixture and is defined by

γ (T) = \frac{C_{p} (T)}{C_{v} (T)}

(25)

For each species, the specific heat as a function of temperature is given in the form of least squares coefficients. The general form of this equation is as follows²⁷

\frac{C_{p}^{0}}{R_{0}} = \sum a_{i} T^{q_{i}}

(26)

In equation (24), g(Re_d) and h(M_d) are the two functions

g ({Re}_{d}) = \frac{1 + {Re}_{d} (12.278 + 0.548 {Re}_{d})}{1 + 11.278 {Re}_{d}}

(27)

and

h (M_{d}) = \frac{5.6}{1 + M_{d}} + 1.7 \sqrt{\frac{T_{d}}{T_{c}}}

(28)

where T_d is the droplet temperature and T_c is the temperature of the gas. C_d₀ is the drag coefficient for a Mach number approach to zero, namely, the spherical drag law is given by

C_{d 0} = {\begin{matrix} \frac{24}{{Re}_{d}} (1 + 1 / 6 {Re}_{d}^{2 / 3}) & {Re}_{d} < 1000 \\ 0.424 & {Re}_{d} \geq 1000 \end{matrix}

(29)

and the liquid Reynolds number Re_d in equation (29) is given by

{Re}_{d} = \frac{2 ρ | u + u' - v_{d} | r}{μ_{d}}

(30)

Injection model

An essential component of the theoretical model is the manner in which the effect of the initial droplet size distribution and the injection position on droplet dispersion is represented. In previous studies, an average droplet size and distribution have been used to simulate injected droplets.^1,17,18 However, strong turbulence structures will develop inside the injector at high-speed flow condition, the average droplet size and distribution are no longer reasonable.³² Furthermore, the droplet sizes of injected parcels are strongly dependent on atomizer configurations. The spray of the liquid jet is significantly influenced by droplet size and distribution. To accurately depict the physical characteristic of injected droplets, a statistical injection model has been created to predict the velocity, size, and position of droplets in the section of injector.²⁹ The spray injector is set as a point in space with specified droplet mass flow rate, droplet distributions, sizes, velocities, temperatures, and breakup parameters.

The diameter of droplets in the parcel is randomly chosen from a specific probability density distribution

\int_{- \infty}^{ζ_{1}} f (x) dx = F (ζ_{1})

(31)

The investigation of the influence on droplet size of atomizer design indicates that the radius of the droplet with different atomizers can usually be characterized by the PDF.³³ In this work, four types of commonly used PDFs f(x) are available: uniform, chi-squared, Nukiyama–Tanasawa, and Rosin–Rammler.

The position vector (x, y, z) is given by

\begin{matrix} x = x_{0} + d \cos ϑ n_{1 x} + d \sin ϑ n_{2 x} \\ y = y_{0} + d \cos ϑ n_{1 y} + d \sin ϑ n_{2 y} \\ z = z_{0} + d \cos ϑ n_{1 z} + d \sin ϑ n_{2 z} \end{matrix}

(32)

where d and $ϑ$ are, respectively, the radial distance from the injector axis and the angle, satisfying the uniform distribution in (0, d_inj) and (0, $2 π$ ), where d_inj is the diameter of the injector. The vector (x₀, y₀, z₀) specifies the centroid of the injector. The vectors (n_1x, n_1y, n_1z)(n_1x, n_1y, n_1z) and (n_2x, n_2y, n_2z) represent two orthogonal directions. These two vectors form a plane perpendicular to the incident direction of injected droplets.

The droplet injection velocity is randomly chosen from a Gaussian distribution with standard deviation $u_{inj}$ , where u _inj is the initial droplet velocity. By generating a series of random variables, the initial injection conditions (droplet sizes, velocities, and locations) for all injected parcels can be established.

Model validation

To validate the simulation method, the structure of water jets injected into a Mach 1.94 crossflow is numerically studied as the baseline case, where well-defined boundary conditions and experimental data are available in Lin et al.’s¹⁴ work. The free-stream conditions are defined by the pressure, temperature, and density which are 29 kPa, 304.1K and 0.331 kg/m³, respectively. The temperature and the radius of injected parcels are set to constant values of 298.15K and $5 \times 10^{- 5} m$ . The jet-to-air momentum flux ratio q₀, which is defined as $q_{0} = ρ_{liq} v_{liq}^{2} / ρ_{\infty} v_{\infty}^{2}$ , is derived from the physical diameter of the nozzle orifice. According to the experimental conditions, a constant q₀ = 10 is used to imply the constant liquid flow rate. Based on experimental data, a calibration process, including mesh density and spray submodels, is carried out to determine the best parameters to properly reproduce the reference case which is performed by the Eulerian–Lagrangian method with Jameson scheme and SST turbulence model in this article. Geometrical parameters in numerical calculation are depicted in Figure 1.

Figure 1.

Geometrical parameters of the test case. All length dimensions are in millimeter.

Grid sensitive test

The cell size is determined in the calibration process because the Lagrangian submodels are highly dependent on it. The rectangular duct is discretized by structured, hexahedral mesh elements (Figure 2). In this calibration study, three different meshes are evaluated to model the reference case condition where a coarse grid with $0.5 \times 1.0 \times 1.0 - mm$ cells, a refined grid with $0.5 \times 0.8 \times 0.8 - mm$ cells, and a refined grid with $0.5 \times 0.65 \times 0.65 - mm$ cells are generated. The mesh density is increased in the vicinity of the fuel injection orifice, which is located at the lower wall. These grid resolutions are typical for viscous flow simulations for y⁺ < 10.

Figure 2.

Computational grid.

In Figure 3 (left), spray penetration heights predicted by the proposed method using the three different meshes are compared with each other. Very close penetration heights are obtained with two refined grids, while the numerical result of the coarse mesh is relatively lower. In addition, the gas density distribution along the x-axis at the center plane of y-axis and z = 0 mm on the lower wall is also compared and analyzed for better illustrating the model predictions in Figure 3 (right). Gas density distributions show that predictions of the coarse mesh have a higher density in the near-injector region and further downstream region. It has to be noted that with higher spatial resolution, the liquid parcels contained in one computational cell becomes less, and as a result, the momentum transfer in the near-injector region from the fuel to the air in the refined grid is slower. Therefore, a higher penetration is reached compared to the coarse case. When the most refined mesh is used, the axial gas density in the near-injector region of the jet is slightly lower where liquid droplets are transported farther from the injector nozzle because of the smaller external aerodynamic drag on the parcels, which lead to a higher liquid penetration height of the specific simulation. Similar trends have been obtained by Desantes et al.’s³⁴ work. To save the computational cost reasonably, the mesh with the grid resolution of $0.5 \times 0.8 \times 0.8 m m^{3}$ is used since the difference between two refined meshes is small enough (within 3%).

Figure 3.

Comparisons of spray penetration heights and x-axial gas density distributions in the near-injector region with three different meshes.

Breakup model

Spray breakup process plays the role in mixing the fuel with air which increases the droplet surface area for rapid evaporation and combustion. It significantly affects the droplet size, penetration height, flow field structure, especially in supersonic crossflows. The performance of various breakup models in high-speed crossflows is mainly studied in this section. A further investigation and comparison is carried out for three classical breakup models such as the TAB model, the Reitz wave model, and the K-H/R-T hybrid model for better understanding of the physical mechanisms. Detailed ambient conditions and the breakup models of all cases are presented in Table 1.

Table 1.

Conditions for spray atomization in supersonic crossflows.

No.	Incoming air conditions		Injection conditions		Breakup model
No.	Pressure (kPa)	Temperature (K)	q ₀	Temperature (K)	Breakup model
1	29.00	304.1	10.0	298.15	No breakup
2					TAB model
3					Reitz wave model
4					K-H/R-T hybrid model

K-H/R-T: Kelvin–Helmholtz/Rayleigh–Taylor.

Figure 4 shows penetration heights of atomizing liquid jets when the physical flow time t > 0.2 ms where the liquid jet can be considered as stable state. The penetration height is compared with the experimental correlation function, which is developed by Liu,³⁵ namely

\frac{h}{d} = 3.14 q_{0}^{0.29} {(\frac{x}{d})}^{0.38}

(33)

where $h$ is the penetration height and $d$ is the nozzle orifice diameter.

Figure 4.

Comparisons of spray penetration heights: (a) no breakup, (b) TAB model, (c) Reitz wave model, and (d) K-H/R-T hybrid model.

Results show in Figure 4(a) depict that the case without breakup model significantly overpredicted the penetration height. This can be explained that the sizes and types of droplets remain unchanged when breakup contributions are ignored which results in a higher momentum and a slower momentum transfer from the fuel to the air for every droplet is obtained. Therefore, liquid droplets are transported farther away from the injector nozzle where the spray penetration overprediction occurred. The cases with different breakup models are shown in Figure 4(b)–(d). The case with Reitz wave model predicts a higher penetration height which outclasses that of K-H/R-T model which has a good match with the experimental data. This can be explained that in the near-injector region, where the fluid acceleration has a different direction to the density gradient, the R-T instabilities dominate. The K-H model is applied immediately after the injection region to provide the aerodynamic instabilities that will begin to grow on the droplet surface, which shear off smaller secondary droplets from the parent droplet surface compared to Reitz wave model and result in earlier breakup and hence lead to a shorter liquid penetration height. Additionally, it has to be noted that with the TAB breakup model, only one oscillation mode corresponding to the lowest order harmonic in the near-injector region can be represented which is the most important oscillation mode near the injector. But other modes dominate immediately once the droplets are transported farther from the injector, where the Weber number becomes larger in high-speed flow.¹⁶ This explains that droplet appears at the bottom of the flow field in Figure 4(b), where Weber number is small and the TAB model results in earlier breakup compared to K-H/R-T model. As previously mentioned, the high Mach number drag law is used for accurately predicting the drag acting on droplets in the supersonic regime. For better understanding the influences of drag coefficients on calculation results, the simulation with the high Mach number drag law and the spherical drag law³¹ is compared as shown in Figure 4(a)–(d). The predicted values of the simulation with high Mach number drag law are slightly higher than that of spherical drag law in terms of the penetration height profiles. This can be justified because of the fact that in the near-injector region corresponding to the higher relative velocity between the droplet and gas, a larger drag coefficient is predicted by the high Mach number drag law and hence results in earlier breakup, and the liquid penetration becomes shorter.

Figure 5 shows the Sauter mean diameter (SMD) profiles along the z-axis at the center plane of y-axis when x = 50 mm from injector origin. The droplet size predicted by the TAB model is the smallest, while the K-H/R-T model obtains the most accurate droplet size compared with the experimental data although there are still distinguishable differences on SMD distribution trends. So, the K-H/R-T model is the most reasonable one among these three breakup models.

Figure 5.

SMD distribution of droplets.

Injection model

The proposed injection model is under the assumption where the droplet dispersion is dependent on the droplet size and initial injection conditions. The injection condition of each parcel is established by PDF where the state of the parcel is characterized by parameters such as position, velocity, and radius. Effects of the initial droplet size distribution and the injection position on droplet dispersion are studied in this section. A total of four types of PDFs are adopted (uniform, chi-squared, Nukiyama–Tanasawa, and Rosin–Rammler). In previous statements, simulation results of the K-H/R-T hybrid model are more accurate than that of the remaining two models. Therefore, the K-H/R-T hybrid model is used in combination of the injection model for further description on procedures of the fuel spray where the mesh used for this case remains unchanged. Table 2 illustrates representative equations of the PDFs.

Table 2.

Liquid spray atomization modeling table.

Probability density function	Representative equations
Uniform distribution	$f (x) = 1 / (b - a)$
Chi-squared distribution	$f (x) = (x^{(k / 2) - 1} e^{- x / 2}) / (2^{k / 2} Γ (k / 2))$
Nukiyama–Tanasawa distribution	$f (x) = k x^{p} e^{(- {(x / m)}^{n})}$
Rosin–Rammler distribution	$f (x) = 1 - e^{(- {(x / m)}^{n})}$

Thousands of dispersed droplets have been released from the center of the inlet section of the injector with the spray evolution in time and space. The number of droplets in each size interval $Δ D$ is counted, recorded, and divided by the total number of droplets in the sample. Results from two representative positions are discussed in the following section: (a) the near-nozzle regime corresponds to the liquid column breakup that occurs and (b) the atomization regime where droplets are much smaller than the initial droplet size at 400 mm from the injector. The droplet radius distribution at 2.0 ms after the start of the injection is shown in Figure 6 in the form of a histogram (bar chart), where the liquid penetration is deemed to be a steady state. The ordinate corresponding to each size interval is known as the number frequency, f_n(D_i). The droplet size is identified by the discrete number frequency distribution. The sum of each bar is unity since the number in each size category has been divided by the total number in the sample. The distribution has been normalized, that is³¹

\sum_{i = 1}^{N} f_{n} (D_{i}) = 1

(34)

where N is the total number of intervals.

Figure 6.

Droplet size distribution: discrete number frequency distribution and discrete cumulative number distribution: (a) uniform distribution, (b) chi-squared distribution, (c) Nukiyama–Tanasawa distribution, and (d) Rosin–Rammler distribution.

The cumulative distribution is another commonly used method to quantify droplet size, which is the sum of the frequency distribution as also shown in Figure 6 and associated with size D_k³¹

F_{n} (D_{k}) = \sum_{i = 1}^{D_{k}} f_{n} (D_{i})

(35)

Results shown in Figure 6 refer to all the droplets contained in a disk area near the injector which is perpendicular to the x–y plane, where the thickness is equal to the minimum cell size and has a radius of 0.5 mm. And the cylindrical centerline of the disk coincides with the center of the injector. The black bars represent the number frequency distribution characterized by PDF scaled by the vertical coordinate on the left, while the gray area shows the cumulative number distribution corresponding to the number frequency distribution scaled by the vertical coordinate on the right. Since the effect of the internal turbulent channel flow dominates at the region near the nozzle exit, which is characterized by the injection model, the profile of droplet size distribution satisfies the specific probability density distribution.

One of the most important and difficult measurements in the high-speed crossflow is the cross-sectional droplet distribution along the flow direction, because it is directly related to the mixing process and eventually the combustion performance.¹⁸ Results shown in Figure 7 refer to all droplets contained in a slice with thickness of 0.5 mm and located at 400 mm from the injector normal to the x-axis. Different droplet size distributions are observed depending on the initial injection conditions. Strong aerodynamic forces acting on the droplet surface distort the liquid jet and promote the liquid atomization, and the droplet diameters are ranging from 20 to 40 µm, which are in accordance with experimental result.¹⁶

Figure 7.

Droplet size distribution at 400 mm from the injector nozzle: (a) uniform distribution, (b) chi-squared distribution, (c) Nukiyama–Tanasawa distribution, and (d) Rosin–Rammler distribution.

The number of droplets contained in Figure 7 and the mean radius of them are shown in Table 3. Quantitatively, the number of droplet prediction is the highest and the mean radius prediction is the smallest in the case of Nukiyama–Tanasawa distribution which can be explained that the percentage of small droplets injected compared to the other three cases is the highest. The smaller droplets, which are ripped out of the surface of both the liquid jet and fragment drops, are immediately transported with the gas flow due to their small inertia and highly accelerated by the co-flowing air which will reach the free-stream velocity on very short path lengths. In contrast, bigger drops with high inertia approach the free-stream velocity toward the end of the duct. With respect to spray spread, bigger drops are able to penetrate deeper into the duct because of their higher momentum which results in the higher penetration predicted by Nukiyama–Tanasawa distribution.

Table 3.

Comparison of droplet parameters for baseline conditions in the slice.

Probability density function	Number of droplets	Diameters( $μ m$ )
Uniform distribution	8315	25.5
Chi-squared distribution	20,846	28.7
Nukiyama–Tanasawa distribution	35,615	24.4
Rosin–Rammler distribution	17,160	28.1

Figure 8 shows the numerical prediction of two-dimensional droplet distributions with different droplet size PDFs. It is clear that the case with Rosin–Rammler distribution shows a wider and higher droplet distribution. The droplet size of the newly generated parcel with Rosin–Rammler distribution is larger than other three cases as shown in Figure 6(d). Thus, the parcels are less decelerated, which augment the penetration and create a spray where the liquid is more spread in the z-axial direction.

Figure 8.

Comparisons of spray penetration heights: (a) uniform distribution, (b) chi-squared distribution, (c) Nukiyama–Tanasawa distribution, and (d) Rosin–Rammler distribution.

Figure 9 shows the SMD profiles along the z-axis at the center plane of y-axis when x = 50 mm from injector origin. The Nukiyama–Tanasawa distribution more accurately predicted the size of child droplets, and the same behavior is shown also in the uniform distribution case and chi-squared distribution case. However, with the Rosin–Rammler distribution, the droplet size distribution is monotonously increasing which outclasses that of the experimental data. But calculation results are not able to catch the finer spray structure presented in experimental measures which is probably due to the nature of numerical method where time-averaged fields are obtained.

Figure 9.

SMD distribution on varying probability density function.

As presented previously, in Figures 6 –9 and Table 3, numerical results with four different droplet size PDFs are analyzed. Modifying the statistical description can affect the numerical results of spray procedure of free water injected into a supersonic crossflow, which are the penetration heights, SMD distributions, and so on. The primary atomization process and the exchange of momentum are greatly different, and different penetration heights and SMD distributions are achieved with different PDFs. According to comparisons of penetration heights and SMD distributions among four PDFs in this section, the prediction with the Nukiyama–Tanasawa distribution is closely related to the experimental value, showing a better spatial behavior reproducing the near-nozzle and far-nozzle flow structures.

Conclusion

An integrated Eulerian–Lagrangian method is proposed in this work which mainly emphasizes influences on numerical results when different breakup models and different statistical injection models are employed. The spray atomization of a water jet injecting into a supersonic crossflow is numerically studied by the method proposed in this article.

On one hand, detailed studies include the performance of three classic breakup models in the supersonic regime, TAB model, Reitz wave model, and K-H/R-T hybrid model, for breakup model evaluation. By comparing with the TAB model and the Reitz wave model, the K-H/R-T hybrid model shows the greatest accuracy on spray penetration height and droplet size distribution prediction, where the K-H/R-T hybrid model is selected for further study. On the other hand, studies of influences on the numerical results illustrated by the number frequency and the cumulative fraction with different PDFs are carried out. The droplet dispersion is also investigated with various droplet sizes and initial injection conditions. The droplet dispersion in the vicinity of the injector can be finely described which is related to initial injection conditions. Better reproduce of spray characteristics such as initial droplet size distribution and jet structure which will result in significant impacts on improving the engine combustion process can be obtained with reasonable PDF. Furthermore, the cross-sectional drop size distribution in the near-nozzle and far-nozzle regions is also discussed for better understanding the influence on the spray structure with different initial injection conditions. Numerical results show that simulations with the Nukiyama–Tanasawa PDF have the best match with the experimental data in terms of the local droplet size distribution and penetration heights.

Above all, the integrated Eulerian–Lagrangian method has a great ability to describe the entire spray process including the breakup that occurs near the nozzle exit in the primary breakup region.

Footnotes

Handling Editor: Takahiro Tsukahara

Declaration of conflicting interests

The author(s) declared no potential conflicts of interest with respect to the research, authorship, and/or publication of this article.

Funding

The author(s) received no financial support for the research, authorship, and/or publication of this article.

References

Xiao

Dianat

McGuirk

. Large eddy simulation of liquid-jet primary breakup in air crossflow. AIAA J 2013; 51: 2878–2893.

Segal

. The scramjet engine: processes and characteristics. New York: Cambridge University Press, 2011.

Zhao

et al . An experimental study of kerosene combustion in a supersonic model combustor using effervescent atomization. P Combust Inst 2005; 30: 2859–2866.

Brinckman

Hosangadi

Ahuja

et al . A CFD methodology for liquid jet breakup and vaporization predictions of compressible flows. In: Proceedings of the 46th AIAA aerospace sciences meeting and exhibit, Reno, NV, 7–10 January 2008. Reston, VA: AIAA.

Pai

Bermejo-Moreno

Desjardins

. Role of Weber number in the primary breakup of liquid jets in crossflow. College Park, MD: American Physical Society, 2009.

Pitsch

Desjardins

. Detailed numerical investigation of turbulent atomization of liquid jets. Atomization Spray 2010; 20: 311–336.

Khare

Chen

et al . Phenomenology of secondary breakup of Newtonian liquid droplets. In: Proceedings of the 50th AIAA aerospace sciences meeting, Nashville, TN, 9–12 January 2012. Reston, VA: AIAA.

Pastor

Garcia-Oliver

Pastor

et al . One-dimensional diesel spray modeling of multicomponent fuels. Atomization Spray 2015; 25: 485–517.

Fan

Wang

. A marker-based Eulerian-Lagrangian method for multiphase flow with supersonic combustion applications. Int J Mod Phys Conf Ser 2016; 13: 42.

10.

Kirkendal

Fuller

et al . Breakup processes of liquid jets in subsonic crossflows. J Propul Power 1997; 13: 67–73.

11.

Reitz

. Modeling atomization processes in high-pressure vaporizing sprays. Atomization Spray Technol 1987; 3: 309–337.

12.

Mazallon

Dai

Faeth

. Primary breakup of nonturbulent round liquid jets in gas crossflows. Atomization Spray 1999; 9: 291–311.

13.

Sallam

Aalburg

Faeth

et al . Breakup of aerated-liquid jets in supersonic crossflows. In: Proceedings of the 42nd AIAA aerospace sciences meeting and exhibit, Reno, NV, 5–8 January 2004. Reston, VA: AIAA.

14.

Lin

Kennedy

Jackson

. Structures of water jets in a Mach 1.94 supersonic crossflow. In: Proceedings of the 42nd AIAA aerospace sciences meeting and exhibit, Reno, NV, 5–8 January 2004. Reston, VA: AIAA.

15.

Kyoung

Lin

Lai

. Spray atomization of liquid jet in supersonic cross flows. In: Proceedings of the 43rd AIAA aerospace sciences meeting and exhibit, Reno, NV, 10–13 January 2005. Reston, VA: AIAA.

16.

Bravo

Kweon

. A review on liquid spray models for diesel engine computational analysis. Technical report no. 6932, May 2014. Adelphi, MD: Army Research Laboratory.

17.

Madabhushi

. A model for numerical simulation of breakup of a liquid jet in crossflow. Atomization Spray 2003; 13: 413–424.

18.

Lin

Lai

. Breakup modeling of a liquid jet in cross flow. Int J Automot Technol 2011; 12: 489–496.

19.

Herrmann

. Detailed numerical simulations of the primary atomization of a turbulent liquid jet in crossflow. J Eng Gas Turb Power 2010; 132: 61506–61510.

20.

Xiao

Wang

Sun

et al . Large eddy simulation of liquid jet primary breakup in supersonic air crossflow. Int J Multiphas Flow 2016; 87: 229–240.

21.

Huang

. Transverse jet in supersonic crossflows. Aerosp Sci Technol 2016; 50: 183–195.

22.

Huang

Yan

. Progress in research on mixing techniques for transverse injection flow fields in supersonic crossflows. J Zhejiang Univ Sci Appl Phys Eng 2013; 14: 554–564.

23.

O’Rourke

Amsden

. The TAB method for numerical calculation of spray droplet breakup. In: Proceedings of the international fuels and lubricants meeting and exposition, Toronto, ON, Canada, 2 November 1987. Warrendale, PA, USA: SAE Technical Paper.

24.

Reitz

. Computer modeling of sprays. Pittsburgh, PA: Spray Technology Short Course, 1996.

25.

Chang

. Hydrodynamics of liquid jet sprays. PhD Thesis, University of Wisconsin–Madison, Madison, WI, 1991.

26.

Blazek

. Computational fluid dynamics: principles and applications. Amsterdam: Elsevier, 2005.

27.

Liu

. Numerical methods for complex combustion flowfields. PhD Thesis, Nanjing University of Aeronautics and Astronautics, Nanjing, China, 2009.

28.

Gosman

Loannides

. Aspects of computer simulation of liquid-fueled combustors. J Energy 2012; 7: 482–490.

29.

Gawne

Liu

Bao

. Modelling of plasma–particle two-phase flow using statistical techniques. Surf Coat Technol 2005; 191: 242–254.

30.

Crowe

Babcock

Willoughby

. Drag coefficient for particles in rarefied, low mach–number flows. Proc Int Symp Two Phase Syst 1972; 6: 419–431.

31.

Crowe

Sommerfeld

Tsuji

. Multiphase flows with droplets and particles. 2nd ed.Boca Raton, FL: CRC Press; Taylor & Francis, 2012.

32.

Xiao

Dianat

Mcguirk

. LES of turbulent liquid jet primary breakup in turbulent coaxial air flow. Int J Multiphas Flow 2014; 60: 103–118.

33.

Wang

. Modeling and numerical simulations of internal combustion process of liquid rocket engines. Beijing, China: National Defense Industry Press, 2012.

34.

Desantes

Garcia-Oliver

Pastor

. A comparison of diesel sprays CFD modelling approaches: DDM vs ε-Y Eulerian atomization model. Atomization Spray 2016; 26: 713–737.

35.

Liu

. Numerical and experimental investigation of fuel spray in supersonic cross flow. PhD Thesis, Beihang University, Beijing, China, 2009.