Abstract
The atomization process of a liquid jet in supersonic crossflow with a Mach number of 1.94 was investigated numerically under the Eulerian-Lagrangian scheme. The droplet stripping process was calculated by the KH (Kelvin-Helmholtz) breakup model, and the secondary breakup due to the acceleration of shed droplets was calculated by the combination of the KH breakup model and the RT (Rayleigh-Taylor) breakup model. In our research, the existing KH-RT model was modified by optimizing the empirical constants incorporated in this model. Moreover, it was also found that the modified KH-RT breakup model is applied better to turbulent inflow of a liquid jet than laminar inflow concluded from the comparisons with experimental results. To validate the modified breakup model, three-dimensional spatial distribution and downstream distribution profiles of droplet properties of the liquid spray in the 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.
1. Introduction
Within the combustion chamber of the Scramjet engine, a transverse liquid jet is injected into the supersonic airflow at a certain flow rate. Once exposed to the incoming airflow with a high Mach number, the jet deforms immediately and breaks apart into small fragments rapidly under strong aerodynamic forces [1]. None of any clear boundaries can be utilized to divide the whole breakup process into several specific stages. However, generally, the process of the liquid column breaking apart into pieces of initially large droplets can be recognized as “the primary breakup,” while the process of initially large droplets further breaking up into small-sized droplets is usually recognized as “the secondary breakup.” Gorokhovski [2] pointed out that the primary breakup has a significant impact on the subsequent droplet formation and dispersion, and yet, it is of huge challenge to understand the whole physical mechanisms within for various complicated structures in the flow field as well as the large time and space spans which thus make it even harder to investigate the microstructures both experimentally and numerically. On the other hand, since our goal is to produce small droplets downstream with appropriate sizes to increase the evaporation and mixing rates of the liquid fuel, one of the most significantly practical motivations to investigate the atomization process is to determine the conditions in which the desired final fragment sizes can be acquired. As noted by Tryggvason [3], the highest airflow velocity does not always guarantee the smallest droplet diameter. Therefore, the secondary breakup process should also be clearly understood to meet the actual engineering demands.
Plenty of experimental investigations have been launched regarding the liquid jet in supersonic crossflow. With the application of laser technology and computer technology, many advanced testing technologies have been developed, such as phase Doppler particle analyzer (PDA/PDPA/LDV), laser scattering technology/laser holography technology, laser-induced fluorescence technology (PLIF), and particle image velocity instrument (PIV). The liquid spray properties measured by different instuments turn out differently. Koh [4] found that the mass distribution of PDA is obviously smaller than that of optical imaging. Lin [5] contrasted the effects on the penetration height caused by different optical approaches and pointed out that PDA is more sensitive to thin liquid mist thus measuring a higher penetration height than other approaches such as high-speed photography and schlieren. Lin [6] studied the structures of water jets injected into a crossflow by utilizing a two-component PDPA and successfully discovered the S-type distribution profiles of droplet size in the downstream position. What is more, the correlations of the operational conditions and droplet properties can be investigated, and the whole physical mechanisms of the atomization process can be further understood by numerical calculations. Currently, numerical approaches to investigate the whole gas-liquid mixing and interacting process in supersonic conditions can be generally divided into two types. One is to capture the movement of gas-liquid interface based on the Eulerian scheme such as VOF (volume of fluid) and LS (level set) approaches. These approaches are at high accuracy but huge computational costs. The other is to track the position of single droplets by integration of time based on the Lagrangian scheme. In this paper, the lagrangian method is focused on because compared with the Eulerian scheme, it can remarkably reduce computation costs and combine itself with breakup models compatibly; although, the initial size and velocity distributions of large droplets need to be determined.
In the beginning stage of the investigation, out of the restriction of experimental conditions, most of the research was focused on the dependence of deformation and breakup of droplets on those dimensionless numbers such as Reynolds number and Weber number [7–9]. At present, it is still extremely difficult to simulate the whole atomization process from continuous liquid column to uniform spray plume with one single physical model. Therefore, existing breakup models have been coupled and improved to calculate breakup processes of liquid jets to acquire good agreement with experimental results. Currently, the theories of surface waves induced by the Kelvin-Helmholtz (KH) and the Rayleigh-Taylor (RT) instabilities are considered to explain the essential mechanisms of the liquid breakup process mostly [10]. The KH-RT breakup model is based on the surface wave theories and takes both KH and RT instabilities into account. At present, it is considered as the most appropriate breakup model to truly reveal the breakup mechanisms of the liquid phase in the supersonic environment. Liu [11] studied the breakup parameters of the KH-RT model and pointed out that in the supersonic environment, parameters controlling the breakup time have quite limited impacts for the very short breakup time in supersonic airflows. Yang et al. (Yang, Zhu, Sun, and Chen, 2017) [12] improved the KH-RT breakup model by considering the compressible effects of the gaseous environment and found that the penetration height of the improved model agreed well but the spread and size distributions of the liquid spray were still different from the experiment results. Li et al. (Li, Wang, Sun, and Wang, 2017) [13] adopted the KH breakup model to simulate the droplet stripping process near the nozzle and coupled RT breakup and TAB (Taylor analogy breakup) models to simulate the secondary breakup of droplets. In his research, a LES code program of two-phase flow was carried out under high resolution grid, and the downstream atomization characteristic profiles obtained were in good agreement with the experimental results.
The inflow turbulence of a liquid jet was widely investigated both experimentally and numerically. It was noted to control the primary breakup process of a liquid jet in supersonic crossflow and thus further control the downstream droplet properties. Mazallon et al. (Mazallon, Dai, and Faeth, 1999) [14] and Sallam et al. [15] (Sallam, Ng, Sankarakrishnan, Aalburg, and Lee, 2006) investigated the laminar liquid jet injected into a uniform gaseous crossflow by using pulsed shadowgraphy and pulsed holography. It was found that surface waves formed on the upstream side of the laminar jet column, and that the wavelength decreased with the increase of the gaseous Weber number which indicated that the surface waves in the primary breakup process originate from RT instability. Xiao et al. [16] (Xiao, Dianat, and McGuirk, 2013) also drew the same conclusion with a two-phase-flow large-eddy simulation. It was also concluded that it is liquid rather than gaseous turbulence that determines the initial liquid-jet instability and interface characteristics. Lee et al. [17] (Lee, Aalburg, Diez, Faeth, and Sallam, 2007) discovered that the SMD of droplets stripped from the surface of a turbulent liquid jet column was not influenced by the crossflow, and thus the liquid turbulence controls the primary breakup process.
In this paper, the atomization process of a liquid jet in supersonic crossflow with a Mach number of 1.94 was investigated numerically using the Lagrangian method. In this method, large droplets with initial size and velocity distributions were injected into a supersonic crossflow in substitution for a continuous liquid jet. Then, the droplets broke apart into small fragments calculated with the modified KH-RT breakup model and eventually formed the spray plume. This paper is organized as follows: In 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.
2. Mathematical Models
2.1. Eulerian-Lagrangian Governing Equations
2.1.1. Gas-Phase Dynamic Equation
where is the gas density, is the gas velocity, is the static pressure, is the gas dynamic viscosity, is the momentum source term of droplets, is the drag coefficient, is the droplet velocity, is the droplet density, is the droplet diameter, is other forces per unit gas mass, is the droplet relaxation time [18], and is the relative Reynolds number .
2.1.2. Liquid-Phase Dynamic Equation
where is the drag force per unit droplet mass, and 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 can be computed fromwhere includes accelerations due to all other forces except drag force, and and represent droplet velocity and gas velocity at the old location.
2.2. Turbulence Model
In this paper, the (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 is splitted into a mean and fluctuating components . For compressible flow, the density varies so widely that a mass-weighted average (called Farve average) is usually preferred . Any quantity may be splinted into mean and fluctuating components as
The average balance equations of continuous phase are written aswhere are specific heat capacity of constant volume and pressure, respectively. The turbulent kinetic energy is written aswhere is the molecular viscosity, and is thermal conductivity.
To provide closures for the unknown turbulent kinetic energy , Menter [19] proposed the model. The classical turbulence model has high prediction ability with high Reynolds number region in far field while the classical turbulence model has better prediction ability and more stable numerical properties with the near-wall region. The turbulence model combines the advantages of both and turbulence models switched by a mixing function . The transport equations of turbulence model are as follows:
The mixing function is defined bywhere is the vertical distance from the wall surface, is the positive part of the cross diffusion term, and is determined by
All the coefficients can be calculated in a uniform form bywhere is the corresponding coefficients of and turbulence models, respectively.
The eddy viscosity coefficient is defined bywhere is the eddy vector.
In Eq. (23), the first component is the eddy viscosity coefficient of the turbulence model, and the second component is acquired from the one-equation turbulence model based on the characteristics of shear stress in laminar boundary layers.
2.3. Modified KH-RT Breakup Model
The KH-RT breakup model is based on the linearized instability theory and takes both KH and RT breakup into account. It is often used to simulate high Weber number sprays. Within this model, a length-limited liquid core in the near nozzle region is assumed as shown in 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 corresponding to wave length is obtained as below [20].where , , , , and .
The radius of child droplet stripped from the cylindrical jet is
is the KH breakup time calculated bywhere is the KH breakup radius constant, and 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 belowwhere is the acceleration of the droplet, is the RT breakup radius constant, and is the RT breakup time constant.
When the KH-RT breakup model is adopted to track wave growth on the surface of droplets, it is often beneficial to adjust the constants to acquire optimal simulation results compared with the experiments. According to the research by Yang et al. (Yang, Le, He, and Chen, 2012) [23], the optimal values of breakup constants , , , and in supersonic flows are summarized below as 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 rather than time constant . It is also found that the RT radius constant has much less influence on the downstream SMD than the KH radius constant . 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 is probably because the present liquid core length is larger than expected.
(a)
(b)
Therefore, using the idea of controlling variables, the optimal values of breakup constants for supersonic simulation of two-phase interactions calculated by FLUENT software are shown in Table 2.
2.4. Inflow Turbulence of a Liquid Jet
In this paper, the DPM (discrete particle model) is used under the Eulerian-Laglangian frame via the fluid mechanic simulation software FLUENT. In the DPM model, the “Blob” model was applied by injecting large droplets into the free airflow in substitution for a continuous liquid jet. Then, the behavior and trajectory of droplets are calculated by breakup models and interaction with the crossflow. Although with the exact same orifice diameter specified, the inflow turbulence of a liquid jet is discovered to have significant effects over the primary breakup process by affecting the initial droplet sizes and velocities. Therefore, to accurately predict the spray characteristics, the correct turbulent state of the internal nozzle flow must be specified. As can be seem from 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 is the controllable upstream and downstream pressure, respectively.
(a)
(b)
(c)
In the rest of this section, the varied initial droplet size and velocity of different flow turbulence are focused on. And it is shown that with the same nozzle diameter and flow mass rate specified, the turbulent (single-phase) liquid flow fits better with the experiment results than the laminar (flipped) liquid flow.
2.4.1. Determination of the Initial Droplet Size Distribution
According to our investigation, the initial droplet diameter distribution is closely related to the nozzle state. To indicate the connection, a two-parameter RR (Rosin-Rammler) distribution is used to represent the droplet diameter distribution as below, characterized by the most probable droplet size and a spread parameter . The mass fraction of droplets with diameters greater than is calculated by
The first parameter required to specify the droplet size distribution is the most probable droplet size . For a single-phase nozzle flow, the correlation of Wu et al. (Wu, Tseng, and Faeth, 1992) [24] is applied to calculate SMD 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.where is the radial integral length scale at the jet exit based upon fully developed turbulent pipe flow, and 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 , where 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 jetwhere 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 . The values for the spread parameter are determined from past modeling experience and experimental observations. Table 3 lists the values of 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 in the central plane distinguished by the flow turbulence with different initial droplet diameter distributions. Therefore, the turbulent liquid jet fits better with the experimental results in the aspect of the droplet size distribution.
(a)
(b)
2.4.2. Determination of Initial Droplet Velocity
For a single-phase nozzle, the estimate of exit velocity 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 is the contraction coefficient by Nurick’s [27] fit, is the internal upstream pressure of the nozzle, is the internal downstream pressure of the nozzle, and 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. Figure 7 indicates the higher velocity of laminar flow compared with turbulent flow.
(a)
(b)
3. Computational Conditions
Lin et al. (Lin and Kennedy, 2002) [6] carried out the test research of a water liquid jet injected into supersonic crossflow with a Mach number of 1.94 and acquired abundant experimental data using a two-component phase Doppler particle analyzer (PDPA), which has been commonly simulated to verify new numerical models. In this section, the exact case mentioned above was simulated by FLUENT software to validate our physical models.
Considering the computational cost, the calculation domain was set to be a rectangular region near the injector with . 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 , the vicinity of where was locally refined as shown in Figure 8. The 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/m3, viscosity of kg/(m‧s), and surface tension of 0.072 N/m. The momentum flux ratio is set to be constant 7, and other detailed parameters of liquid and gas are given in Tables 4 and 5.
4. Results and Discussion
4.1. Model Validation
4.1.1. Spatial Distribution of the Liquid Spray
Spray penetration height is an important atomization characteristic to indicate the liquid-gas mixing effect. 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 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.
4.1.2. Centerline Profiles of Droplet Properties
Centerline distribution profiles of droplet properties in the free stream direction can be normalized by local penetration height to obtain universal curves. Figure 12 illustrates the normalized SMD distribution along the -axis at different 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 -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 grid points in the , , and directions. And the comparison shows considerable consistency in the SMD, absolute, and relative velocities distributions at the position of 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 are contrasted in Figure 14. It can be found that the simulated gas-phase vortex field is remarkably weaker than the actual vortex field because the droplet volumes are neglected in the Eulerian-Lagrangian scheme. As a result, in simulations, it is the aerodynamic shear force that simply controls the downstream SMD distribution. That is to say that the closer to the free airflow, the smaller the size of droplets is. However, the experimental results show that the maximum droplet size appears at both the top of the spray plume and the location due to the complicated liquid-gas interaction and complex gas-phase vortex field.
(a)
(b)
4.2. Grid Independence Verification
To investigate the reliability of current grids, a higher-resolution set of grid with a total number of 2093184 cells is used to calculate this problem. The top view of this refined set of grid is screenshot as 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.
4.3. Modeling of the SMD Distribution
Due to actual demands, the scale effects of droplet properties need to be explored by several simulated cases with different airflow Mach numbers and nozzle diameters. Although the S-shaped profile of the SMD distribution can inevitably increase the modeling difficulty, a new modeling approach is presented in this section to help settle the SMD quantification problem. To solve that, modified power functions are established which show good accuracy in predicting the SMD distribution. Dozens of numerical cases are launched to acquire statistically enough data for mathematical analysis with the aid of the multivariate nonlinear fitting method. Detailed operational conditions are listed below, and the expressions of the SMD distribution in the free stream direction as well as the errors are exhibited.
4.3.1. Operational Conditions
In this section, a total of 54 sets of operating conditions are designed, among which the nozzle diameter varies in 0.3, 0.5, 1, 2, and 3 mm, and the airflow Mach number varies in 2.1, 3, and 4. It is worth noting that since the flux momentum ratio and the airflow Weber number 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 and always constant. Detailed settings of parameters are shown in Table 6.
4.3.2. Modeling Instruments
In this section, the traditional power function for developing penetration height formulasis modified into a new-form power function to describe the SMD distribution in the free stream directionwhere is the boundary layer thickness, and is the nozzle diameter, Ma is the free stream Mach number, is the distance from the nozzle. , , , , and are parameters to be determined.
These unknown parameters are optimized by using the multivariate nonlinear fitting approach to correlate the local SMD at different positions with various nozzle diameters and airflow Mach numbers. Eventually, the formulas of SMD are given separately at different positions in Table 7, and the MSE (mean squared error) represents the fitting error. As can be seen, the fitting results are acceptably accurate within a certain error range. Figure 22 is the schematic diagram of the fitting results.
4.3.3. Result Analysis
The whole expression of SMD can be analyzed from two components. For a fixed position, the former component is controlled by the nozzle diameter and airflow Mach number while the latter component remains constant. Since gaseous Weber number and flux momentum ratio remain unchanged in different settings of operational conditions, the latter component should be regarded as inherent influence of and . Mathematically, it is quite reasonable because it was proven that has a significant effect on the droplet diameter. About the former component, it can be observed that SMD is roughly proportionate with and negatively proportionate with . Therefore, it can be concluded that apart from the component affected by the changeless parameters, the changeable component of the value of SMD is closely related to the nozzle diameter and airflow Mach number.
5. Conclusions
This research mainly focuses on the atomization process of a liquid jet injected into supersonic crossflows. The discrete phase model and modified KH-RT breakup model are adopted to investigate the inflow turbulence of a liquid jet which can significantly affect the distribution profiles of downstream droplet properties. Also, a new method to describe the SMD distribution is presented as our main innovative point. Abundant numerical data is acquired, and main conclusions are summarized below.(1)The inflow turbulence of a liquid jet has significant influence on the primary breakup, which is reflected on the droplet initial size and velocity. And a turbulent jet is found to fit better with experimental results than a laminar one(2)In the supersonic environment, it was discovered that the KH breakup has more control on the whole breakup process than RT breakup. The KH breakup size constant can significantly affect the downstream SMD distribution profile(3)The distribution of SMD can be well described with a modified power function presented in this paper. The mathematical results implied that gaseous Weber number and flux momentum ratio still have considerable effects on the final droplet size in the supersonic airflow, while the scale effect of nozzle diameter and airflow Mach number can also be significant. And their effects on the SMD distribution are fairly comparable
Data Availability
The datasets generated during and/or analyzed during the current study are available from the corresponding author on reasonable request.
Conflicts of Interest
The authors declare that they have no conflicts of interest.