Sage Journals: Discover world-class research

Abstract

A pure two-fluid model was used for investigating transverse liquid jet to a supersonic crossflow. The well-posedness problem of the droplet phase governing equations was solved by applying an equation of state in the kinetic theory. A k-ε-k_p turbulence model was used to simulate the turbulent compressible multiphase flow. Separation of boundary layer in front of the liquid jet was predicted with a separation shock induced. A bow shock was found to interact with the separation shock in the simulation result, and the adjustment of shock structure caused by the interaction described the whipping phenomena. The predicted penetration height showed good agreement with the empirical correlations. In addition, the turbulent kinetic energies of both the gas and droplet phases were presented for comparison, and effects of the jet-to-air momentum flux ratio and droplet diameter on the penetration height were also examined in this work.

Keywords

Liquid injection supersonic crossflow kinetic theory penetration height

Introduction

The description of liquid jet injection into a supersonic crossflow has become an important research area as the development of high-speed flight vehicles. Such a process is usually encountered in thrust vector control, external burning on projectiles, and particularly in the design of scramjet engines. For liquid fuel propulsion system, the combustion efficiency depends strongly on the breakup, atomization, and evaporation of liquid fuel jet.^1–5 However, the complexity of mechanism in these physical processes makes it difficult for detailed analysis from a theoretical point, especially in compressible regimes.

Extensive experimental investigations have been carried out previously to study the spray structures.^6–8 Kush and Schetz⁹ used water as the liquid jet normally injected into a supersonic gas crossflow and found the liquid surface layers appear to be different in size, shape, and depth with varying jet-to-air momentum flux ratio. Schetz et al.¹⁰ employed photographic observations to study the jet column breakup mechanism and found that the formation and behavior of the waves surrounding the jet column are insensitive to large variations in surface tension and viscosity. Lin et al.¹¹ investigated the structures of liquid jets injected vertically into a Mach 1.94 crossflow experimentally and used Phase Doppler Particle Analyzer (PDPA) measurement to establish distributions of droplet and spray plume properties. They found the normalized distribution profiles for droplet size and droplet velocity exhibit S and mirrored-S shapes, respectively, which are caused by the presence of the bottom floor. Perurena et al.¹² conducted an experimental investigation of a water jet into a crossing Mach 6 air flow, penetration heights for different shaped injectors were compared, and the liquid fracture was shown to be influenced by the bow shock–separation shock interaction. Based on the previous works, a jet plume model has been devised, and the schematic of the liquid jet into a high-speed crossflow is shown in Figure 1.

Figure 1.

Schematic of liquid jet injection into supersonic flow.

However, limited numerical investigations have been seen in the literatures, which is due to the complexity of the problem both in physical and theoretical aspects. There are two common approaches to investigate the two-phase flow problem. One is the Lagrangian approach, which takes the droplets as mass points, and traces the droplets to obtain their locations and velocities individually. The other one is the Eulerian approach, which is a two-fluid model, and the droplet phase plays a role as continuum fluid. Comparing with the former approach, the Eulerian approach solves the governing equations in the Eulerian frame, which could avoid massive cost of computation for the droplet trajectories. The Eulerian approach has been widely used in simulations of incompressible regimes.^13–17 Nevertheless, for modeling compressible multiphase flow, unlike the Navier–Stokes equations of the gas phase, the droplet phase governing equation system is found to be non-hyperbolic, which will lead to an ill-posed problem and produce non-physical results.^18–21 In order to obtain a well-posed equation system, the kinetic theory of granular flow^22–25 is introduced in the current study.

The kinetic theory of granular flow is on the basis of kinetic theory of non-uniform dense gas as described by Chapman and Cowling.²⁶ The main idea of the theory is to treat the particles in a dense state as continuum fluid and apply the analogy between particles in granular flow and molecules in dense gas flow. The velocity distribution of a single particle is assumed to follow the Maxwellian distribution, and a pseudo temperature for the particle phase is derived to describe the random motion of the particles, which is called granular temperature²⁷ with respect to the thermal temperature in gas phase. In particular, an equation of state is derived from the kinetic theory and modifies the characteristics of the governing equations of the particle phase. For the liquid injection in this article, the jet is supposed to be atomized into dense droplets plume. By assuming the breakup and atomization processes are completed quickly after the jet enters the supersonic flow, the kinetic theory is implemented consequently in the simulation of droplet phase.

The scope of this article is to construct a well-posed equation system for the droplet governing equations based on the kinetic theory and carry out the numerical simulation of transverse liquid jet to a supersonic crossflow. For turbulence modeling of the two phases, the k-ε model was applied to the gas phase with compressible correction,²⁸ while the one-equation k_p model was used to describe the turbulence in droplet phase.²⁴ Characteristics analysis was employed to investigate the well-posedness of the governing equation system, and the hyperbolic nature based on the kinetic theory was then demonstrated. As the authors know, it is the initial work to extent the kinetic theory with k_p turbulence model into the simulation of compressible multiphase flow. The penetration height of the liquid jet from the numerical results was compared with the empirical correlations derived from the experimental data, and good agreement was obtained. The predicted distribution of the turbulent kinetic energy of both phases indicated the shear layer is the main mixing area for the liquid jet and the supersonic crossflow air. In addition, the jet-to-air momentum flux ratio showed to be important in predicting the penetration height, while the influence of the droplet diameters on the numerical results was slight.

Methodology

Governing equations for gas phase

In this work, the Navier–Stokes equations in conservative form for the compressible gas phase were solved, and the expressions are presented as follows:

Gas phase continuity equation

\frac{\partial}{\partial t} (α_{g} ρ_{g}) + \frac{\partial}{\partial x_{j}} (α_{g} ρ_{g} v_{g j}) = 0

(1)

Gas phase momentum equation

\frac{\partial}{\partial t} (α_{g} ρ_{g} v_{g i}) + \frac{\partial}{\partial x_{j}} (α_{g} ρ_{g} v_{g i} v_{g j}) = - α_{g} \frac{\partial p_{g}}{\partial x_{i}} + \frac{\partial τ_{i j}}{\partial x_{j}} + f_{g p i}

(2)

Gas phase total energy equation

\begin{array}{l} \frac{\partial}{\partial t} (α_{g} ρ_{g} E_{g}) + \frac{\partial}{\partial x_{j}} [α_{g} ρ_{g} v_{g j} (E_{g} + \frac{p_{g}}{α_{g} ρ_{g}})] \\ = \frac{\partial v_{g i} τ_{i j}}{\partial x_{j}} - \frac{\partial}{\partial x_{j}} (λ_{g} \frac{\partial T_{g}}{\partial x_{j}}) - Q_{g p} \end{array}

(3)

where the subscript g represents the gas phase, and α_g, ρ_g, and v_g are the volume fraction, density, and velocity of the gas, respectively. The total energy is given by

E_{g} = e_{g} + \frac{1}{2} v_{gi}^{2} + k

(4)

and k is the turbulent kinetic energy. Besides, p_g and T_g are the gas phase pressure and thermal temperature, which follow the equation of state

p_{g} = (γ - 1) (E_{g} - \frac{1}{2} v_{gi}^{2} - k)

(5)

and γ is the ratio of specific heats for the gas phase. The stress $τ_{ij}$ is given by

τ_{ij} = μ_{g} (\frac{\partial v_{gi}}{\partial x_{j}} + \frac{\partial v_{gj}}{\partial x_{i}}) - δ_{ij} \frac{2}{3} μ_{g} \frac{\partial v_{gk}}{\partial x_{k}}

(6)

where µ_g is the gas viscosity coefficient. In equations (2) and (3), f_gpi and Q_gp are the interaction source terms coupling the two phases, which will be discussed in section “Interaction modeling.”

Governing equations for droplet phase

The kinetic theory introduces the concept of granular temperature Θ used to describe the droplet random motion as molecules in the gas phase.^26,27 Based on the kinetic theory,²⁹ the governing equations of the droplet phase can be expressed in conserved form in the Eulerian frame as follows:

Droplet phase continuity equation

\frac{\partial}{\partial t} (a_{p} ρ_{p}) + \frac{\partial}{\partial x_{j}} (α_{p} ρ_{p} v_{pj}) = 0

(7)

Droplet phase momentum equation

\frac{\partial}{\partial t} (α_{p} ρ_{p} v_{p i}) + \frac{\partial}{\partial x_{j}} (α_{p} ρ_{p} v_{p i} v_{p j}) = - α_{p} \frac{\partial p_{g}}{\partial x_{i}} + \frac{\partial Π_{i, j}}{\partial x_{j}} - f_{g p i}

(8)

Droplet phase granular temperature equation

\begin{array}{l} \frac{\partial}{\partial t} (α_{p} ρ_{p} Θ) + \frac{\partial}{\partial x_{j}} (α_{p} ρ_{p} v_{p j} Θ) \\ = \frac{2}{3} Π_{i, j} \frac{\partial v_{p j}}{\partial x_{i}} + \frac{2}{3} \frac{\partial}{\partial x_{j}} (Γ_{Θ} \frac{\partial Θ}{\partial x_{j}}) - \frac{2}{3} ω \end{array}

(9)

Droplet phase energy equation

\begin{array}{l} \frac{\partial}{\partial t} (α_{p} ρ_{p} C_{p}^{p} T_{p}) + \frac{\partial}{\partial x_{j}} (α_{p} ρ_{p} v_{p j} C_{p}^{p} T_{p}) \\ = \frac{\partial}{\partial x_{j}} (C_{p}^{p} \frac{μ_{p}}{σ_{T}} \frac{\partial T_{p}}{\partial x_{j}}) + Q_{g p} \end{array}

(10)

where the subscript p represents the droplet phase, α_p, ρ_p, v_p, and T_p are the volume fraction, density, velocity and thermal temperature of the droplet, respectively, and $C_{p}^{p}$ is the droplet specific heat capacity. In addition, Π_i,j is the shear stress of droplet phase

\begin{array}{l} Π_{i, j} = - p_{p} δ_{i, j} + ζ_{p} δ_{i, j} \frac{\partial v_{p j}}{\partial x_{j}} + \\ μ_{p} [(\frac{\partial v_{p i}}{\partial x_{j}} + \frac{\partial v_{p j}}{\partial x_{i}}) - \frac{2}{3} δ_{i, j} \frac{\partial v_{p j}}{\partial x_{j}}] \end{array}

(11)

where µ_p and ζ_p are droplet phase viscosity and bulk viscosity, respectively. It should be noted that according to the kinetic theory, ζ_p can be expressed as

ζ_{p} = \frac{4}{3} α_{p}^{2} ρ_{p} d_{p} g_{0} (1 + β) \sqrt{\frac{Θ}{π}}

(12)

and µ_p is defined in dilute and dense phase conditions, respectively

Droplet phase dilute viscosity

μ_{p, dilute} = \frac{5 \sqrt{π}}{96} ρ_{p} d_{p} Θ^{1 / 2}

(13)

Droplet phase dense viscosity

\begin{array}{l} μ_{p, d e n s e} = \frac{10 ρ_{p} d_{p} \sqrt{π Θ}}{96 (1 + β) g_{0}} {[1 + \frac{4}{5} (1 + β) g_{0} α_{p}]}^{2} \\ + \frac{4}{5} α_{p}^{2} ρ_{p} d_{p} g_{0} (1 + β) \sqrt{\frac{Θ}{π}} \end{array}

(14)

The differences between the two expressions lie in the collisions and subsequent interactions among droplets, which contribute to the production of the droplet phase pressure. The droplet phase pressure can be obtained by an equation of state based on the kinetic theory as

p_{p} = α_{p} ρ_{p} [1 + 2 (1 + β) α_{p} g_{0}] Θ

(15)

with g₀ and β representing radial distribution function and the coefficient of restitution

g_{0} = \frac{3}{5} {[1 - {(\frac{α_{p}}{α_{p, \max}})}^{\frac{1}{3}}]}^{- 1}

(16)

Note that g₀ should reach infinite when the volume fraction α_p equals the package limit α_p,max, which is 0.6 in the present study. When the droplet phase is in dilute condition and g₀ then approaches 0, the “kinetic pressure” is obtained as described by Gidaspow²⁹

p_{p} = α_{p} ρ_{p} Θ

(17)

The coefficient of restitution β ranges from 0 to 1 depending on the collision regime, and set to be 0.9 in the present study. In equation (9), Γ_Θ is the transport coefficient of the granular temperature, which is similarly defined for the dilute and dense conditions respectively as

Droplet phase dilute granular conductivity

Γ_{Θ, dilute} = \frac{75 \sqrt{π}}{384} ρ_{p} d_{p} Θ^{1 / 2}

(18)

Droplet phase dense granular conductivity

\begin{array}{l} Γ_{Θ, d e n s e} = \frac{150 ρ_{p} d_{p} \sqrt{π Θ}}{384 (1 + β) g_{0}} {[1 + \frac{6}{5} (1 + β) g_{0} α_{p}]}^{2} \\ + 2 α_{p}^{2} ρ_{p} d_{p} g_{0} (1 + β) \sqrt{\frac{Θ}{π}} \end{array}

(19)

and ω is the dissipation term of the granular temperature

ω = 3 (1 - β^{2}) α_{p}^{2} ρ_{p} g_{0} Θ (\frac{4}{d_{p}} \sqrt{\frac{Θ}{π}} - \frac{\partial v_{pk}}{\partial x_{k}})

(20)

where d_p is the droplet diameter.

Interaction modeling

In this study, two-way coupling was established through the source terms in the governing equations between the gas phase and droplet phase, including the momentum exchange and heat transfer. The viscous drag force is caused by the velocity slip of the two phases and expressed in terms of the mechanical droplet response time τ_rp, which is also called droplet motion relaxation time

f_{gp} = \frac{α_{p} ρ_{p}}{τ_{rp}} (v_{gi} - v_{pi}) and τ_{rp} = \frac{d_{p}^{2} ρ_{pm}}{18 μ_{g}} \frac{1}{C_{D}} \frac{24}{R e_{p}}

(21)

where ρ_pm is the material density of the droplets, and C_D is the drag coefficient consisting of three separate correlations based on the range of the droplet Reynolds number Re_p

C_{D} = {\begin{matrix} 24 / R e_{p} & R e_{p} < 1 \\ 24 / R e_{p} (1 + \frac{1}{6} {Re}_{p}^{2 / 3}) & 1 \leq R e_{p} \leq 1000 and R e_{p} = \frac{ρ_{g} | v_{g} - v_{p} | d_{p}}{μ_{g}} \\ 0.44 & R e_{p} > 1000 \end{matrix}

(22)

In equation (3), Q_gp represents the heat transfer between the two phases per unit volume, and the expression for Q_gp is³⁰

Q_{gp} = α_{p} ρ_{p} π d_{p} N u_{p} λ_{g} (T_{g} - T_{p})

(23)

where the correlation of droplet Nusselt number is given by

N u_{p} = 2.0 + 0.6 {Re}_{p}^{0.5} P r^{0.33}

(24)

and Pr is the Prandtl number $\Pr = μ_{g} C_{p}^{g} / λ_{g}$ .

Turbulence modeling

Turbulence of the two phases was considered in the current simulation. The two equation k-ε turbulence model with Wilcox compressibility correction was used for the gas phase.²⁸ In order to incorporate the turbulence of the droplet phase, the one-equation k_p turbulence model was employed, while k_p denotes the turbulent kinetic energy of the droplet phase, and the concept is consistent with the turbulent kinetic energy of the gas phase. Therefore, the combined k-ε-k_p turbulence model²⁴ was presented for the compressible multiphase flow in this article.

Droplet phase turbulent kinetic energy

\begin{array}{l} \frac{\partial}{\partial t} (α_{p} ρ_{p} k_{p}) + \frac{\partial}{\partial x_{j}} (α_{p} ρ_{p} v_{p j} k_{p}) \\ = \frac{\partial}{\partial x_{j}} (\frac{μ_{p}}{σ_{p}} \frac{\partial k_{p}}{\partial x_{j}}) + G_{p k} - α_{p} ρ_{p} ε_{p} \end{array}

(25)

where the production and dissipation term in the turbulent kinetic energy equation are given by

\begin{array}{l} G_{p k} = μ_{p} (\frac{\partial v_{p i}}{\partial x_{j}} + \frac{\partial v_{p j}}{\partial x_{i}}) \frac{\partial v_{p i}}{\partial x_{j}} and \\ ε_{p} = - \frac{2}{τ_{r p}} (c_{p} \sqrt{k k_{p}} - k_{p}) \end{array}

(26)

and the turbulent viscosity is defined as

μ_{p} = α_{p} ρ_{p} c_{μ p} k_{p}^{0.5} k^{1.5} / ε

(27)

The empirical constants for the two-phase governing equations above are listed in Table 1 It should be mentioned that the empirical constants were chosen in the range of common values as presented in the previous studies,^15,24,25 where the constants have been examined with the experimental data.

Table 1.

Empirical constants in the governing equations.

σ_T	σ_p	β	c_p	c_µp
1.0	0.35	0.9	0.75	0.01

Characteristics analysis

As in the compressible gas phase, the equation of state in the kinetic theory connected the granular temperature with pressure in the droplet phase. Therefore, the continuity equation, momentum equation, and granular temperature equation of the droplet phase were combined and solved simultaneously in the numerical simulation. Since the gas phase governing equations have already been demonstrated as a hyperbolic, thus well-posed system,³¹ characteristics analysis was employed to investigate the well-posedness of the droplet phase governing equations.

To be concise, the hyperbolic characteristic in x direction of the governing equations was chosen to be demonstrated for representative in the present study, and the same analysis could be easily extended into the other two dimensions. Assuming the conservative variables of the droplet phase are included in W_p, and the convection terms of the governing equations are included in F_p, the expressions are then given as

W_{p} = {[\begin{matrix} α_{p} ρ_{p} & α_{p} ρ_{p} u_{p} & α_{p} ρ_{p} v_{p} & α_{p} ρ_{p} w_{p} & α_{p} ρ_{p} Θ \end{matrix}]}^{T}

(28)

and

F_{p} = {[\begin{matrix} α_{p} ρ_{p} u_{p} & α_{p} ρ_{p} u_{p} u_{p} + p_{p} & α_{p} ρ_{p} v_{p} u_{p} & α_{p} ρ_{p} w_{p} u_{p} & α_{p} ρ_{p} Θ u_{p} \end{matrix}]}^{T}

(29)

Note that the pressure in the droplet phase shear stress in equation (8) is displaced to the convection term, and in this way, the characteristics analysis for the droplet phase governing equations can be carried out subsequently. By applying the chain rule to the derivative of the convection term F_p with respect to the conservation variable W_p

\frac{\partial F_{p}}{\partial x} = \frac{\partial F_{p}}{\partial W_{p}} \frac{\partial W_{p}}{\partial x}

(30)

the terms on the left side of the equations (7)–(9) can be expressed in form of the Burgers’ equation

\frac{\partial W_{p}}{\partial t} + A (W_{p}) \frac{\partial W_{p}}{\partial x} = 0

(31)

and the coefficient matrix A(W_p) is the Jacobian matrix of the equation system

A (W_{p}) = \frac{\partial F_{p}}{\partial W_{p}} = [\begin{matrix} 0 & 1 & 0 & 0 & 0 \\ - u_{p}^{2} & 2 u_{p} & 0 & 0 & R_{p} \\ - u_{p} v_{p} & v_{p} & u_{p} & 0 & 0 \\ - u_{p} w_{p} & w_{p} & 0 & u_{p} & 0 \\ - Θ u_{p} & Θ & 0 & 0 & u_{p} \end{matrix}]

(32)

where $R_{p} = 1 + 2 (1 + β) α_{p} g_{0}$ , and u_p, v_p, w_p are the velocities in three coordinate directions. The major reason of characteristics analysis is that the choice of numerical method and the stability of the equation system are intimately connected with the nature of the partial differential equations, which has been described by the Lax Equivalence Theorem.³² This theorem states that the stability is both a necessary and a sufficient condition for convergence if a finite difference approximation to a well-posed linear initial-value problem is consistent. Von Neumann analysis also predicts that no stable numerical scheme can be found at least for the linear system, since the partial differential equations are not hyperbolic or ill-posed.²⁹

Based on the concept of diagonalization, a hyperbolic system is often defined as a system with real eigenvalues and diagonalizable coefficient matrix.³¹ Therefore, diagonalization of the Jacobian matrix was tried to be conducted, and five real eigenvalues were obtained then

λ_{1} = u_{p} - a_{p}, λ_{2} = u_{p} + a_{p}, λ_{3} = λ_{4} = λ_{5} = u_{p}

(33)

where $a_{p} = \sqrt{R_{p} Θ}$ is called as “sound speed” of the droplet phase, and it is a function of the granular temperature. It is interesting to note that the expression of the “sound speed” of the droplet phase is similar with the sound speed in the gas phase, which is also proportional to the square root of temperature. Through the above analysis, the equation system of droplet phase is demonstrated a hyperbolic system, and a well-posed problem for droplet phase based on kinetic theory can be constructed.

The governing equations of the gas and droplet phases were solved by finite differential method. The inviscid flux for the gas phase was computed using Roe-type Riemann solver. For the droplet phase, because of the application of kinetic theory, hyperbolic nature of the governing equation system was obtained, and the same numerical solver can be applied to the inviscid flux as in the gas phase. The face-states on the computational grid were reconstructed by Monotone Upstream-Centred Schemes for Conservation Law (MUSCL) method. To avoid un-realistic solutions, a simple modified Harten formulation was applied to satisfy the entropy condition for the Roe scheme.³³ For the viscid flux, second-order accurate central difference scheme was used. A two-step Runge–Kutta method with second-order precision was applied to integrate the equations temporally. To ensure the temporal convergence, the time step is set to be 1.0 × 10⁻⁸ s. Considering the Courant-Friedrichs-Lewy (CFL) condition, it is a quite small integration time step for explicit time marching in transient simulations.

Results and discussion

System description

The test domain is based on the experimental facility in the Virginia Tech supersonic blow down wind tunnel.³⁴ The experiment was conducted under the inflow Mach number 3.0 configuration. The simulations in the present study were carried out under two-dimensional condition, and the computational domain was 0.2 m × 0.18 m as shown in Figure 2. The supersonic air flow entered into the test section in the direction parallel to the flat plate, and water was used in the experiment as the transverse liquid jet, which was injected from an orifice located 80 mm downstream of the leading edge with its diameter d₀ of 2.0 × 10⁻³ m. A grid system of 251 × 121 was used with refinement near the injection orifice and the flat plate wall, and the size of the minimum grid cell is 5 × 10⁻⁴ m as shown in Figure 3.

Figure 2.

Schematic of the liquid injection into the supersonic crossflow.

Figure 3.

Grid set 251 × 121.

In the Eulerian approach, the uniform droplet size is used since only one droplet size can be calculated by one set of droplet phase governing equations. According to the liquid spray experiment by Less and Schetz,³⁵ the diameter of the droplets after atomization process mostly varies from 7 × 10⁻⁶ m to 100 × 10⁻⁶ m. In this article, the breakup and atomization of the liquid jet were assumed to be completed quickly in the supersonic crossflow.¹¹ The breakup and atomization processes were ignored and the uniform droplet diameter according to the liquid spray experiment was used in the present simulations, which was assumed to be 30 × 10⁻⁶ m. The effect of the droplet diameter on the numerical results will be discussed later in section “Effect of the momentum flux ratio and droplet diameter on penetration height.” The inflow stagnation pressure in the experiment maintained at 4.5 × 10⁵ ± 2% Pa and the stagnation temperature was close to the ambient air. The primary parameter used in injection investigations is the jet-to-air momentum flux ratio

q = ρ_{j} v_{j}^{2} / ρ_{\infty} v_{\infty}^{2}

(34)

where $ρ_{j}$ and $ρ_{\infty}$ are the density of the injected water and the inflow air, respectively. In the present study, the computation was conducted under the experimental condition q = 6, and the corresponding mass flow rate of the injected water was maintained at 0.09 kg/s as in the experiment. The inlet boundary was treated as a supersonic inflow condition, where the variables in the computation are given directly. The variables on the outlet and upside boundary were obtained through extrapolation. The bottom wall was set to be a no-slip boundary condition.

Penetration height

Since the penetration height of fuel jet injection into the combustor has great effect on the diffusive mixing of fuel and freestream air, which is needed for the preconditions of chemical reaction, many studies have been devoted to formulate empirical correlations for the liquid jet’s penetration height. The expressions for the penetration height are commonly described as functions of the momentum flux ratio, spatial distance from the orifice, and orifice diameter, such as^12,35

\frac{h}{d_{0}} = 1.15 {(q)}^{0.5} \ln (1 + 6 x / d_{0})

(35)

\frac{h}{d_{0}} = 3.75 {(q)}^{0.414} {(\frac{x}{d_{0}})}^{0.239}

(36)

\frac{h}{d_{0}} = 4.73 {(q)}^{0.3} {(\frac{x}{d_{0}})}^{0.3}

(37)

\frac{h}{d_{0}} = 3.94 {(q)}^{0.47} {(\frac{x}{d_{0}})}^{0.21}

(38)

Different definitions of the penetration height in experiments can be found in the previous investigations.^{8,11,12,34,35} In this article, the penetration height was obtained by the liquid jet plume boundary based on critical density of 0.1 kg/m³ of the droplet phase. Since the density of droplet phase can be obtained through the droplet phase continuity equation (7), the critical density of the liquid plume boundary is then easily recognized in the numerical results. Instead of tracking the droplets individually in Lagrangian approach, for the Eulerian approach, the number density of the droplets in each grid cell can be represented by the density of the droplet phase. Therefore, the statistical counting process for droplets in each grid cell interacted with the gas phase is unnecessary, avoiding the computation cost increase and the resulted statistical convergence problems. It is very important to demonstrate the grid independence for the numerical simulation.³⁶ Simulations were carried out on three different sets of grid systems, namely 231 × 101, 251 × 121, and 271 × 141. The comparison of the numerical results with the empirical correlations for the penetration height is shown in Figure 4(a) obtained with the grid set 251 × 121. It is noted that the numerical result by the current two-fluid model is well matched with the empirical correlations. The injected liquid column maintains its direction at first. After about 8d₀ transverse distance from the orifice, the jet starts to incline, and the curves of penetration height turn to the inflow direction, with the progressive increasing of jet penetration into the main flow slowing down. In Figure 4(b), the penetration heights resulted from the three grid sets are compared for grid independency analysis. It can be observed that the penetration height predicted is a little higher with coarse grid. However, the overall results do not change appreciably with the grid refinement. The grid set of 251 × 121 is then employed in the following numerical simulations.

Figure 4.

Penetration height of the numerical results: (a) comparison with empirical correlations and (b) grid independency.

Boundary layer separation

In the present study, time evolutions of the transverse liquid jet to the supersonic flow were obtained. Figure 5 shows the time evolutions of the Mach number contours from 0.4 × 10⁻³ s to 6.4 × 10⁻³ s. It is interesting to note that the boundary layer detachment, together with the bow shock and separation shock, is well described using the two-fluid model. Due to the boundary layer on the flat plate, the supersonic freestream is compressed when it enters into the computational domain, resulting in the existence of the leading edge shock near the inlet boundary. When the injected water penetrates into the supersonic flow, strong drag force acts on the high-speed air due to the high velocity slip between the two phases, and adverse pressure gradient is formed in front of the injection column, which results in the boundary layer separation and gives rise to the separation shock. Note that, initially, the separation region is small, and the separation shock and bow shock are weak. As the water injected into the main flow, the interaction force of the two phases becomes stronger, and alterations in pressure distribution near the injection are induced. The separation point is gradually moving upstream and makes the separation shock inclined, creating a further change in bow shock angle. The continuous adjustment of the shock structure is called whipping phenomena¹² and has been presented numerically by the current simulation. In addition, it can be observed that the leading shock behaves weak compared with the bow shock and separation shock, and its interaction with the bow shock has little effects on the flow structures. Actually, in the present study, when the boundary layer separation develops further toward the inlet as shown in Figure 5(a)–(e), the leading shock wave was covered up by the induced separation shock wave. With the development of time, as shown in Figure 5(f)–(i), the difference of the time evolutions of the flow structures in the Mach number contour can be found very small, indicating the numerical solution is converged.

Figure 5.

Time evolutions of the Mach number contours: (a) 0.4 × 10⁻³ s, (b) 0.8 × 10⁻³ s, (c) 1.6 × 10⁻³ s, (d) 2.4 × 10⁻³ s, (e) 3.2 × 10⁻³ s, (f) 4.0 × 10⁻³ s, (g) 4.8 × 10⁻³ s, (h) 5.6 × 10⁻³ s, and (i) 6.4 × 10⁻³ s.

As discussed above, separation of the boundary layer can be seen in front of the jet column. For a close inspection, the vector contour of the gas phase velocity is shown in Figure 6. Due to the separation, a circulation zone is formed, and the velocity of the gas phase decreases to subsonic state. The height of the circulation zone is observed to be about 8d₀, the same as the inclination point of the jet penetration curves shown in Figure 6, indicating that the reduced velocity of the gas phase caused by the boundary layer separation could benefit the penetration height.

Figure 6.

Vector contours of the gas phase velocity.

Since the liquid fuel droplets can be captured by the recirculation zone, the pre-mixing of the fuel with air may occur in this area. Furthermore, the low-speed flow in the recirculation zone provides the mechanism for stabilizing the combustion flame, as in the scramjet combustor where liquid fuel is injected into the supersonic crossflow. Figure 7 shows the droplet phase density distributions in different streamwise locations of the recirculation zone, namely x₁ = 0.065 m, x₂ = 0.070 m, and x₃ = 0.075 m, respectively. The present simulation obtains the droplet phase entraining into the recirculation zone in front of the liquid jet. The predicted droplet phase density increases when the streamwise location approaches the jet, and the captured droplets show to concentrate near the bottom wall. Although the captured droplets are limited, the existence of the pre-mixing process plays an important role in holding the flame in supersonic freestream condition.

Figure 7.

Droplet phase density distributions in the recirculation zone.

It can be expected that the liquid injection in the supersonic flow should result in great changes in pressure distribution in the flow field. Figure 8 shows the pressure profiles in different y locations along the direction of the supersonic inflow. The adverse pressure gradient begins with the separation shock and maintains high pressure in the circulation zone as exhibited on the wall at y = 0 m, which reveals the reason for the boundary layer detachment. Moreover, the separation shock gives rise to the first increase of pressure at y = 0.03 m. After that, the pressure increases more rapidly in the shock interaction point, resulting in the highest value in the pressure field. For pressure distribution at y = 0.1 m, the original flow field developed by the inflow condition is not disturbed until the bow shock appears above the jet plume, leading to the peak value in the pressure profile. It can be noted that the shock waves in the flow field may cause strong resistance for the inflow air.

Figure 8.

Pressure profiles along the inflow direction in different vertical locations.

In numerical simulations, the resistance of the liquid jet to the supersonic inflow air is exerted by the drag force in the source terms of the momentum governing equations (2) and (8). As a result, the two-way coupling computation will cause a great amount of total pressure loss for the supersonic crossflow. To evaluate the influence of the drag force, the total pressure recovery ratio, defined as the ratio between the total pressure p₀ in the internal flow field and the initial inflow total pressure p_{0, total}, has been employed in this study. Figure 9 shows the contours of the total pressure recovery ratio. It can be found that there is a large total pressure drop when the supersonic inflow air encounters the liquid jet. In addition, the induced bow shock and separation shock also cause the total pressure loss. Figure 10 shows the distributions of the total pressure recovery ratio along y direction in three different streamwise locations, namely, x₁ = 0.01 m, x₂ = 0.10 m, and x₃ = 0.18 m. In the x₁ location, there is nearly no total pressure loss. However, for locations x₂ and x₃ behind the liquid jet, the total pressure recovery ratio decreases to 0.05 in the lower part of the spray plume. The contributions of the bow shock to the total pressure loss can be also recognized easily in the distributions with a sudden drop identified in the profiles.

Figure 9.

Contours of total pressure recovery ratio.

Figure 10.

Distributions of total pressure recovery ratio in different streamwise locations.

Turbulent kinetic energies for two phases

Figure 11 shows the turbulent kinetic energies for the air flow and liquid jet. It can be seen that the distribution of the turbulent kinetic energies of the two phases is different. The intensity of the gas turbulence seems to be strong in the circulation zone in front of the jet column and also in the shear layer between the liquid and air stream. This phenomenon can be explained by the large velocity gradient of the gas phase, which appears in these regions and results in strong turbulence there. However, for the droplet phase, the initial part of the jet from the orifice is not disturbed by the gas phase evidently. As the jet penetrates further, in the region where the jet starts to incline to the supersonic inflow direction, the turbulent kinetic energy of the droplet phase is increasing. Note that the shear layer is the main area for the momentum exchange between the two phases, and the velocity of the droplet phase experiences great changes, which may account for the intensified turbulence of the droplet phase there. Therefore, the shear layer is found to possess great significance in the mixing of the fuel jet with the supersonic air flow.

Figure 11.

Turbulent kinetic energies of the gas phase and droplet phase: (a) gas phase and (b) droplet phase.

Effect of the momentum flux ratio and droplet diameter on penetration height

In the transverse gas injection regimes, there are many crucial parameters that have important impacts on the transverse flow field, such as the jet-to-momentum flux ratio, the geometric configuration of the injection orifice, the number of injection orifice, and the injection angle. Among these parameters, the jet-to-crossflow momentum flux ratio implies the energy of the transverse jet in comparison with the crossflow. For the gas injection, the jet-to-crossflow momentum flux ratio shows to be a key parameter which influences the transverse flow field, including the shapes of the bow shock wave, the separation point of boundary layer, and the Ω-shaped vortices.³⁷ For the liquid injection regime, the jet-to-crossflow momentum flux ratio shows its significance by the penetration height of the liquid jet as the empirical correlations (35)–(38) have revealed, which is the most concerned parameter in liquid jet experiments. In Figure 12, the predicted penetration heights of different momentum flux ratio q varying from 3 to 6 are compared. The inflow condition in different cases remains the same, while the jet velocity is modified with respect to q. It is shown that the penetration heights are increasing with the momentum flux ratio. The initial states of the injection out of the orifice are similarly perpendicular. However, the jet with larger injection velocity gains more inertia and is supposed to penetrate deeper into the main flow, leading to the changes of jet inclination points. After the jets inclination, all of the penetration height curves are almost parallel to each other propagating downstream.

Figure 12.

Droplet penetration height of different jet-to-air momentum flux ratios.

In this article, the breakup and atomization of the liquid jet were assumed to be completed quickly, and uniform diameter of droplets was used in the current simulation. To examine the effect of droplet size on the penetration height, in this section, different diameters of droplet phases were tested. In Figure 13, the predicted penetration heights of the liquid jet with droplet diameter varying from 10 × 10⁻⁶ m to 40 × 10⁻⁶ m were obtained under the condition of momentum flux ratio q = 6. It is shown that the larger droplets seem to penetrate a little higher than the smaller ones because of the inertia effect. However, the drag force acted on droplets also increases as the droplets become larger, and the difference in penetration heights appears to be slight as the droplet diameter further increasing. Therefore, it can be the seen that the uniform droplet diameter assumption has small influence on the numerical results.

Figure 13.

Droplet penetration height of different droplet diameters.

Conclusions

A liquid jet into supersonic crossflow has been investigated numerically using the Eulerian approach. The kinetic theory of granular flow was introduced to simulate the droplet phase, and hyperbolic nature of the governing equation system for the droplet phase was demonstrated through characteristics analysis. Therefore, a well-posed problem was established for compressible multiphase flow based on the kinetic theory. The comparison between the penetration height predicted by the numerical work with the empirical correlations derived from experiments was conducted and showed good agreement. The boundary layer separation induced by the adverse pressure gradient in front of the liquid injection was observed with droplet phase captured in the recirculation zone. The whipping phenomenon caused by the bow shock and separation shock interaction was also well described numerically. In addition, a great amount of total pressure loss was predicted when the supersonic freestream enters the spray plume. The present simulation shows that the k-ε-k_p turbulence model is adaptable in compressible multiphase flow and gives a good description of the turbulent flow field of liquid jet injection. The shear layer between the liquid jet plume and the inflow air is found to be an important area for the momentum exchange and mixing process between the two phases. The effects of jet-to-air momentum flux ratio and droplet diameter on the penetration height were examined in the present study, and the liquid jet shows to penetrate deeper into the supersonic flow with larger momentum flux ratio, while the influence of the droplet diameter is slight.

Footnotes

Appendix 1

Academic Editor: Moran Wang

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) disclosed receipt of the following financial support for the research, authorship, and/or publication of this article: This research was supported by the National Natural Science Foundation of China under Grant Nos 51176099 and 51390493.

References

Rachner

Becker

Hassa

. Modelling of the atomization of a plain liquid fuel jet in crossflow at gas turbine conditions. Aerosp Sci Technol 2002; 6: 495–506.

Mathur

Gruber

Jackson

. Supersonic combustion experiments with a cavity-based fuel injector. J Propul Power 2001; 17: 1305–1312.

Vinogradov

Kobigskij

Petrov

MD.

Experimental investigation of kerosene fuel combustion in supersonic flow. J Propul Power 1995; 11: 130–134.

Yang

. A comprehensive study of combustion oscillations in a hydrocarbon-fueled scramjet engine. In: Proceedings of the 45th AIAA aerospace sciences meeting and exhibit, Reno, NV, 8–11 January 2007, paper no. 2007-836, pp.1–12. Reston, VA: AIAA.

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, paper no. 2005-732, pp.1–9. Reston, VA: AIAA.

Sallam

Aalburg

Faeth

. 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, paper no. 2004-970, pp.1–10. Reston, VA: AIAA.

Lin

Kennedy

. Spray structures of aerated-liquid jets in subsonic crossflows. In: Proceedings of the 39th AIAA aerospace sciences meeting and exhibit, Reno, NV, 8–11 January 2001, paper no. 2001-330, pp.1–18. Reston, VA: AIAA.

Lin

Kennedy

Jackson

. Structures of aerated-liquid jets in high-speed crossflows. In: Proceedings of the 32nd AIAA fluid dynamics conference and exhibit, St. Louis, MO, 24–26 June 2002, paper no. 2002-3178, pp.1–12. Reston, VA: AIAA.

Kush

Jr Schetz

JA.

Liquid jet injection into a supersonic flow. AIAA J 1973; 11: 1223–1224.

10.

Schetz

Kush

Joshi

PB.

Wave phenomena in liquid jet breakup in a supersonic crossflow. AIAA J 1980; 18: 774–778.

11.

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, paper no. 2004-971, pp.1–19. Reston, VA: AIAA.

12.

Perurena

Asma

Theunissen

. Experimental investigation of liquid jet injection into Mach 6 hypersonic crossflow. Exp Fluids 2009; 46: 403–417.

13.

Balachadar

Eaton

JK.

Turbulent dispersed multiphase flow. Annu Rev Fluid Mech 2010; 42: 111–133.

14.

Enwald

Peirano

Almstedt

AE.

Eulerian two-phase flow theory applied to fluidization. Int J Multiphas Flow 1996; 22: 21–66.

15.

Zheng

Wan

Qian

. Numerical simulation of the gas–particle turbulent flow in riser reactor based on k–ε–k_p–ε_p–Θ two-fluid model. Chem Eng Sci 2001; 56: 6813–6822.

16.

Abbassi

Namah

Characterization of the speed of a two-phase interface in a porous medium. Math Probl Eng 2005; 6: 641–661.

17.

Zhao

Zhou

. DSMC prediction of particle behavior in gas-particle two-phase impinging streams. Math Probl Eng 2013; 2013: 254082.

18.

Slater

Young

JB.

The calculation of inertial particle transport in dilute gas-particle flows. Int J Multiphas Flow 2001; 27: 61–87.

19.

Toumi

Kumbaro

An approximate linearized Riemann solver for a two-fluid model. J Comput Phys 1996; 124: 286–300.

20.

Sachdev

Groth

CPT

Gottlieb

JJ.

A parallel solution-adaptive scheme for multi-phase core flows in solid propellant rocket motors. Int J Comput Fluid D 2005; 19: 159–177.

21.

Sauel

Daniel

Loraud

JC.

Two-phase flows: second-order schemes and boundary conditions. AIAA J 1994; 32: 1214–1221.

22.

Gidaspow

Bouillard

. Hydrodynamic simulation of gas–solid flow in a riser using kinetic theory of granular flow. Chem Eng J 2003; 95: 1–13.

23.

Lun

CKK

. Kinetic theory for granular flow of dense, slightly, inelastic, slightly rough spheres. J Fluid Mech 1991; 233: 539–559.

24.

Chan

Guo

Lau

KS.

Numerical modeling of gas–particle flow using a comprehensive kinetic theory with turbulence modulation. Powder Technol 2005; 150: 42–55.

25.

Cheng

Guo

Wei

. Modeling the hydrodynamics of downer reactors based on kinetic theory. Chem Eng Sci 1999; 54: 2019–2027.

26.

Chapman

Cowling

TG.

The mathematical theory of non-uniform gases. 3rd ed.Cambridge: Cambridge University Press, 1970.

27.

Jenkins

Savage

SB.

A theory for the rapid flow of identical, smooth, nearly elastic, spherical particles. J Fluid Mech 1983; 130: 187–202.

28.

Liu

Wang

Guo

. Effects of inflow Mach number and step height on supersonic flows over a backward-facing step. Adv Mech Eng 2013; 2013: 147916.

29.

Gidaspow

Multiphase flow and fluidization: continuum and kinetic theory descriptions. San Diego, CA: Academic Press, 1994.

30.

Ling

Haselbacher

Balachandar

Transient phenomena in one-dimensional compressible gas-particle flows. Shock Waves 2009; 19: 67–81.

31.

Toro

EF.

Riemann solvers and numerical methods for fluid dynamics: a practical introduction. 3rd ed.New York: Springer, 2009.

32.

Lax

Richtmyer

RD.

Survey of the stability of linear finite difference equations. Commun Pur Appl Math 1956; 9: 267–293.

33.

Kermani

Plett

. Modified entropy correction formula for the Roe scheme. In: Proceedings of the 39th aerospace sciences meeting and exhibit, Reno, NV, 8–11 January 2001, AIAA paper 2001-0083. Reston, VA: AIAA.

34.

Thomas

Schetz

JA.

Distributions across the plume of transverse liquid and slurry jets in supersonic airflow. AIAA J 1985; 23: 1892–1901.

35.

Less

Schetz

JA.

Transient behavior of liquid jets injected normal to a high-velocity gas stream. AIAA J 1986; 24: 1979–1986.

36.

Huang

Liu

. Influences of the turbulence model and the slot width on the transverse slot injection flow in supersonic flows. Acta Astronaut 2012; 73: 1–9.

37.

Huang

Yan

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

Numerical simulations of transverse liquid jet to a supersonic crossflow using a pure two-fluid model

Abstract

Keywords

Introduction

Methodology

Governing equations for gas phase

Governing equations for droplet phase

Interaction modeling

Turbulence modeling

Characteristics analysis

Results and discussion

System description

Penetration height

Boundary layer separation

Turbulent kinetic energies for two phases

Effect of the momentum flux ratio and droplet diameter on penetration height

Conclusions

Footnotes

Appendix 1

Declaration of conflicting interests

Funding

References