Atomization of Shear Coaxial Liquid Jets

The instability and subsequent atomization of a viscous liquid jet emanated into a high-pressure gaseous surrounding is studied both computationally and experimentally. Liquid water issued into nitrogen gas at elevated pressures is used to simulate the flow conditions in a coaxial shear injector element relevant to liquid propellant rocket engines. The theoretical analysis is based on a simplified mathematical formulation of the continuity and momentum equations in their conservative form. Numerical solutions of the governing equations subject to appropriate initial and boundary conditions are obtained via a robust finite difference scheme. The computations yield real-time evolution and subsequent breakup characteristics of the liquid jet. The experimental investigation utilizes a digital imaging technique to measure resultant drop sizes. Data were collected for liquid Reynolds number between 2,500 and 25,000, aerodynamic Weber number range of 50-500 and ambient gas pressures from 150 to 1200 psia. Comparison of the model predictions and experimental data for drop sizes at gas pressures of 150 and 300 psia reveal satisfactory agreement particularly for lower values of investigated Weber number. The present model is intended as a component of a practical tool to facilitate design and optimization of coaxial shear atomizers.


Introduction
The breakup of a liquid jet in a surrounding gas has many practical applications including shear co-axial injectors of liquid-propellant rocket engines.In such an injector a liquid oxidizer is injected into a gaseous fuel surrounding.The aerodynamic interactions between the liquid surface and gaseous environment bring about the instability and disintegration of the liquid into small drops.These drops disperse and evaporate, forming a mixture with the fuel that can be efficiently burned to produce power.
The focus of the present effort is geared towards the development of a computational/analytical model for the fuel injection and atomization processes that take place in shear coaxial injectors used in liquid propellant rocket engines.The model will effect predictions of atomization characteristics of the round liquid jet produced by coaxial injectors, e.g., drop size, breakup time and length as functions of injector flow conditions.A parallel experimental study is conducted to acquire relevant high-quality data for the purpose of validating the model.Many of the existing models of jet atomization are linearized (e.g., Rayleigh, 1878, 1879, 1882, 1892, Weber, 1931, Tomotika, 1935, Sterling & Sleicher, 1975, Reitz & Bracco, 1982, Lin & Lian, 1990, and Lin & Ibrahim, 1990) and, hence, incapable of providing accurate predictions of resultant spray parameters since atomization is a highly nonlinear phenomenon.In addition, many of these models are riddled with empiricism, which limits their applications to the range of experimental data employed in deriving their empirical formulae.Therefore, it is imperative to advance a jet atomization model that accounts for nonlinear effects and free of empirical relations.
Attempts at formulating nonlinear models of liquid jet atomization have exploited varied mathematical and numerical approaches.Wang(1968), Yuen (1968), Nayfeh (1970), Nayfeh & Hassan (1971), Kakutani et al. (1974), Lafrance (1975), Taub (1976), and Chaudhary & Redekopp (1980) used the method of strained coordinates.Direct numerical solutions of Navier-Stokes equations in their axisymmetric form were obtained by Fromm (1984) and Shokoohi and Elrod (1987) for the viscous jet.Bogy (1978aBogy ( , 1978bBogy ( , 1979aBogy ( , 1979bBogy ( , 1979c) ) employed the Cosserat theory developed by Green & Naghdi (1974) and Green (1976).A weakly nonlinear instability analysis was advanced by Ibrahim & Lin (1992).Ashgriz & Mashayek (1995) adapted the Galerkin finite element method of Keunings (1986).Mansour &Lundgren (1990), andHilbing et al. (1995), Spangler et al. (1995), and Hilbing & Heister (1996, 1998) pioneered the application of boundary-element method to jet atomization studies.A common disadvantage of multidimensional methods is the need of very fine mesh to resolve the liquid-gas interface.This requires a very expensive computer effort, which prohibits practical simulations.Lee (1974) and later Pimbley & Lee (1977) developed a nonlinear-direct-simulation technique that proved to be a simple and practical approach to investigate the nonlinear instability of a liquid jet.Tropey (1989) and Schulkes (1993) utilized variations of Lee's numerical-simulation methodology in their analyses of the instability of inviscid liquid jets.By examining Weber's (1931) linear instability analysis of a viscous liquid jet, Sellens (1992) modified Lee's basic formulation to include viscous terms.He argued that viscous forces become more paramount at smaller jet radius.Since most atomizers are of small scale, the effects of liquid viscosity should be included in numerical simulations of the atomization process to enhance their accuracy.Sellens (1992), however, reported numerical instabilities in his solution, particularly at smaller viscous forces, perhaps due to the central-differencing based finite difference scheme he applied to obtain solutions of the governing equations.He also employed a linearized form of the surface tension term in his formulation.Yi & Reitz (2002) was the first to introduce the aerodynamic effects in the one-dimensional nonlinear modeling approach.Their work also included the effects of viscosity and the exact form of the surface tension term as given by Levich (1962).They employed a semi-explicit time-marching technique to obtain numerical solutions to their governing equations.
In the present work, Yi & Reitz's (2002) approach is adopted in the numerical simulations of a viscous water jet issued in air.However, a novel derivation method and a robust numerical solution technique of the governing equations are presented.A rigorous Total Variation Diminishing (TVD) with flux splitting and characteristic decomposition finite-difference scheme is used to ensure the accuracy and stability of the numerical simulations.Therefore, numerical instabilities reported by other investigators, e.g., Shokoohi and Elrod (1987) and Sellens (1992), in obtaining solutions of the governing equations are eliminated.These instabilities may have originated from the application of central-differencing based finite difference schemes particularly at smaller viscous forces.The present study focuses on predictions and validations of the atomization model by itself and not as a component in a global spray model such as that undertaken by Yi & Reitz (2002).In addition, the independent effects of viscous and aerodynamic forces on atomization are examined.

Theoretical Formulation
The instability of an infinitely long cylindrical liquid jet subject to an initial varicose disturbance is considered.The liquid is assumed to be viscous and incompressible.The axial velocity and the pressure are assumed to be constant over the cross section of the jet and dependent only on axial coordinate z and time t.These assumptions are consistent with the study of the case of long waves (Levich, 1962).Since the surrounding conditions and velocity distribution at each cross section within the jet is uniform, the jet surface will be axisymmetric during the wave growth.In a cylindrical coordinate system moving at the unperturbed (basic) axial jet velocity, U 0 , relative to the gas, the equations of motions may be written in their conservative form as; Continuity: ( ) ( ) Axial Momentum: where u and v are the axial and radial velocities of perturbation in the z and r directions, respectively, t is time, p is pressure, and , ρ ν l l are the respective liquid density and kinematic viscosity.Since the interface is a material surface, the radial velocity component is given by; where h (z,t) is the radius of the disturbed jet.By multiplying each of Eqs. ( 1) and ( 2) by the radial coordinate, r, integrating from 0 to h, the radius of the perturbed jet, and substituting for the radial velocity component from Eq.
(3), a set of unsteady one-dimensional equations is derived; ( ) It should be noted that Eqs. ( 4) and ( 5) above are identical to the ones derived by Sellens (1992) and adopted by Li & Reitz (2002) except for the viscous term.Sellens (1992) derived his equations by considering mass and axial momentum balances on a disc-shaped element of the liquid jet.Since both Sellens ' (1992) and the present equations are approximate due to simplifying assumptions used, discrepancies in the form of some terms may be expected.
The normal stresses due to liquid pressure and surface tension are balanced by the surrounding gas pressure, therefore, where μ l is the liquid dynamic viscosity, p g is the perturbation pressure in the gas, and p σ is the perturbation pressure due to surface tension given by Levich (1962) as; where σ is the surface tension.It should be noted that the viscous term in Eq. ( 6) vanishes by virtue of Eq. ( 3) and the assumption that the axial velocity is independent of the radial distance.Following Yi & Reitz (2002), the gas pressure is taken from the linearized gas equation of motion at the liquid-gas interface (Reitz & Bracco, 1982); where p ∞ is the static pressure of the ambient gas, ρ g is the gas density, k is the disturbance wave number, ( , ) ( , ) z t h z t a η = − , is the amplitude of the disturbance, with a being the undisturbed jet radius.K 0 and K 1 are the modified Bessel functions of the second kind of zeroth and first order, respectively.
Substitution from Eqs. ( 7) and ( 8) into Eq.( 6), carrying out the differentiation of the pressure term with respect to z, and rewriting both Eqs. ( 4) and ( 5) in dimensionless from yields; ( ) where S is the source matrix with two elements, namely, S(1,1) = 0, and S(2,1) given by; ( ) All lengths are normalized by the jet radius, a, velocity is nondimensionalized by the initial velocity, U 0 , and time is made dimensionless by multiplying by the initial velocity and dividing by the jet radius.K is the dimensionless wave number, K = ka, g / ρ ρ ρ l = is the gas to liquid density ratio, We U a / ρ σ Weber number, and Re U a / ν 0 l = is liquid Reynolds number.The terms on the right hand side of Eq. ( 11) correspond to gas, capillary, and viscous forces, respectively.The treatment of the gas, capillary and viscous terms as source terms is to avoid the numerical instability that may originate from these terms as was observed by previous investigators (Shokoohi & Elrod, 1987, Sellens, 1992).
The initial conditions correspond to a varicose disturbance imposed on the unperturbed interface, so that; 0 where η 0 is the dimensionless amplitude of the initial disturbance.Since it is not practical to simulate whole length of the jet, the computational domain is taken to be equal to one wavelength λ = 2π/k and a symmetric boundary condition is applied at the right and left ends of this domain.
The dimensionless continuity and axial momentum equations, given by Eqs. ( 8) and ( 9), may be written in matrix form as: where The present numerical simulations employ a Total Variation Diminishing (TVD) scheme with flux-splitting with characteristic decomposition to obtain solutions of the system of equations given by ( 10).The TVD upwind method uses two characteristic speeds of the convective eigenvector to determine the upwind direction.Therefore, the flux term can be calculated with the following formula; where A is the eigen vector of the system and i is index of the ith grid cell.Accordingly, the finite difference solutions of the dependent variables in the matrix E over the control volume can be obtained through; where ΔZ and ΔT are the spatial and temporal steps, respectively, and superscript m represents index of time marching.For the present computations, the time step was set as 5% of the spatial step size.The computational domain was equal to one wavelength and a uniform mesh system was used.The number of nodes for the one-dimensional grid was chosen to be 20 which rendered a grid-independent convergent numerical solution.

Experimental Hardware and Diagnostics
The Water-Nitrogen Injection Spray Test Facility (WNIST) is located at Marshall Space Flight Center in Huntsville, AL.The purpose of WNIST is to provide a high-pressure environment in which injector sprays can be studied and characterized using propellant simulants of liquid water and gaseous nitrogen.The chamber is modular in design, allowing for a variety of injector types, and scales, to be installed and studied.De-ionized water, fed from a pressurized tank, simulates liquid oxygen (LO x ) at flow rates up to 1 lb m /sec.Gaseous nitrogen is delivered from a 1440-psi regulation system located upstream of the chamber, giving total flow rates up to 3 lb m /sec.The nitrogen's temperature is controllable via a hot water heater system upstream of the chamber.Thus, the gas's density can be set by appropriate selection of temperature and chamber pressure.Using remote-operated flow control valves, flow rates for both the water and nitrogen can be controlled over a wide range.
In order to set chamber pressure, nitrogen gas is flowed through additional feed lines into the working volume of the chamber.A pair of exhaust valves limits chamber venting to provide the desired working backpressure within the chamber.Using this control system, chamber pressures of up to 1400 psia, can be reached within the facility.Judicious selection of water/nitrogen flow parameters can give a range of liquid oxidizer/gaseous fuel propellants simulated at realistic levels of thrust chamber operating pressures and flow rates.These propellant combinations are found in rocket systems such as the Space Shuttle Main Engine (SSME), which uses liquid oxygen and gaseous hydrogen.Figure 1 is a picture of the windowed chamber configuration.
A pair of high-strength fused-silica windows, oriented 180° from each other, gives optical access to the chamber.By using a low-frequency strobe light and a digital camera, shadow graph images can be taken for various injector sprays within the chamber.Images presented in this paper are those from a ¼, geometrically-scaled, SSME injector spray.Figure 2 gives the dimensions of the injector tip.We , and chamber pressure.
Pictures of the spray flow field were taken at various combinations of the parameters listed in Table 1.The camera used is a Kodak DCS ProSLR/n model with digital imaging capability.This camera was used in conjunction with a high-quality strobe light set a 1 Hz flash frequency.The camera and strobe were oriented across the chamber from each other, each set in front of the two windows.Once the injector spray was established at the desired flow parameters, images were taken with the camera.Picture resolution was set at 4500 pixels long by 3000 pixels wide, giving a spatial length scale of 5.5*10 -4 inch/pixel.Four sample images of the spray are given in Figs. 3 and 4. The images are for a fixed Re l = 25,000, g We = 500, and varying chamber pressure, P c .It is clear from Figs. 3 and 4 that increasing the chamber pressure promotes the liquid jet atomization due to greater gas-liquid aerodynamic interaction.The greater the chamber pressure is the higher the gas density and hence the shearing action of gas on the liquid jet surface leading to breaking the liquid bulk into smaller drops.
Maximum and minimum droplet sizes were measured within the focal region of the pictures taken.The largest droplets seen along the centerline of the spray were taken as the maximum droplets.Small satellite droplets, taken either along or around the centerline, were taken as the minimum droplets.Diameters of these droplets were found using the conversion given earlier to go from number of pixels to an absolute value of inches.Details about this measurement process are given in the Results section.

Results
In the present computations, the liquid and gas properties are taken as those of water and nitrogen gas, respectively, to enable simulations of the experimental conditions.The liquid is assumed incompressible and the gas to obey the perfect gas equation.The initial jet radius is equal to half the liquid inner diameter, so that a = 0.046/2 in.= 0.023 in.as given in Fig. 2. Dimensionless parameters, e.g., Weber number, Reynolds number, and gas to liquid density ratio are therefore determined based on these properties, the initial jet radius, and the gas-liquid relative velocities.The fact that the experimental data were obtained as function of aerodynamic (gas) Weber number, g We , based on jet diameter but the model equations involved the liquid Weber number, We , based on jet radius was taken into account by using the conversion g We We / ρ 2 = , as alluded to earlier.No reference of Reynolds number is needed since it is directly related to Weber number for known fluid properties.
The time-history of the disturbance growth rate may be indicative of the convergence of the numerical solutions.Therefore, the temporal variation of dimensionless growth rate is scrutinized in Fig. 5 for dimensionless wave number, K = 0.7 which corresponds to maximum growth rate for Rayleigh jet.The gas-liquid relative velocity, U 0 , is set at a low value of 0.1 m/s so that it will be applicable to Rayleigh jet.The gas to liquid density ratio is taken as ρ = 0.001.The dimensionless growth rate plotted in Fig. 5 is calculated at each time step from the numerical method suggested by Mansour & Lundgren (1990).
( ) The value of η 0 used in Eq. ( 19) is taken as, η 0 = 0.01 as will be discussed later.It is seen in Fig. 5 that the growth rate history exhibits an initial interval of rapid increase followed by a fairly flat region, which indicates a stable solution.
The initial disturbance amplitude used in the computations may have an effect on the growth rate of disturbances.It is desirable to keep the initial disturbance amplitude small to be able to observe the initial stages of the evolution of disturbances on the liquid jet surface and minimize numerical instability.However, employing smaller initial disturbance amplitudes results in longer computational times.To decide upon an acceptable level of the initial amplitude of disturbance its influence on disturbance growth rate is examined in Fig. 6 at K = 0.7, ρ = 0.001, and U 0 = 0.1 m/s.It is seen from Fig. 6 that the dimensionless initial amplitude doesn't have much effect on the dimensionless growth rate when η 0 ≤ 0.01.Therefore, the amplitude of the initial disturbance is set at 1% of the unperturbed jet radius throughout this study.
The evolution of interfacial disturbances on the liquid jet surface is displayed in Fig. 7 at K = 0.7, Re = 1,000 and ρ = 0.001.Note that only Reynolds number is specified since Weber number is directly related to Reynolds number through fluid properties.Figure 7 illustrates the spatial growth of the infinitesimal initial disturbances to a finite size.It should be indicated that the interfacial lines in Fig. 7 are separated by 50 time steps.Recall that the dimensionless time step is set at ΔT = 0.05 ΔZ = 0.05 (λ/an) = 0.05 (2π/Kn), where n is the number of grid points taken as 20.Since the total dimensional domain investigated is one wave length, λ = 2π/K, the dimensionless axial distance plotted in Fig. 7 is normalized by K/π so that it varies from 0 to 2.0 as shown.The dimensionless time of interfacial evolution may be calculated by multiplying the number of lines between the initial and final state of the interface by 50 ΔT.
It is observed in Fig. 7 that the jet surface distortion leads to formation of large size main drops with smaller satellite drops interspaced between them, as has been reported in past studies (Sellens, 1992, Yi & Reitz, 2002).
When the necking portion of the jet touches the centerline, the main and satellite drops will detach from the jet.Using this breakup criterion, the drop radii can be estimated assuming the resultant drops would be spherical.In order to investigate the effect of liquid viscosity on the development of jet disturbance, all parameters are kept the same as those in relation to Fig. 7 except that Reynolds number is decreased to 100 signifying a 10 fold rise in liquid viscosity.The corresponding results shown in Fig. 8 illustrate that increasing the liquid viscosity yield reduction in the magnitude of disturbances for the same evolution time as compared with Fig. 7.This result is expected due to the dissipative nature of liquid viscosity, which leads to suppression of the growth of surface waves.Therefore, it is concluded that the viscous term is properly represented in the present model.The effects of gas-liquid aerodynamic interactions maybe explored by holding all parameters fixed at their values relative to Fig. 7 but increasing the gas to liquid density ratio by one order of magnitude to ρ = 0.01. Figure 9 reveals that the enhanced aerodynamic interaction at higher gas density (or pressure) promotes the destabilization of the liquid jet and yields formation of smaller size main and satellite drops.So, even though the gas pressure term in the model is linearized, it appears to encapsulate the correct aerodynamic behavior.Using dimensions of the injector for scaling, droplet diameters of the spray were determined graphically from the experimental images.These values were evaluated at a reference distance of 1.6" axially downstream of the injector face.The possible largest and smallest drops are identified from the experimental images at this axial position and their sizes estimated by counting the pixels on the high-resolution images.Subsequently, the experimental mean drop size is calculated by a simple arithmetic averaging.It is therefore, logical to compare measured maximum, minimum, and mean drop sizes to computed main, satellite, and mean ones.
Figures 10 and 11 portray such comparisons of measured and computed drop sizes at chamber pressures, P c , of 300 and 150 psia, respectively, as functions of aerodynamic Weber number.Both the experimental and theoretical drop size tend to decrease with aerodynamic Weber number as depicted in Figs. 10 and 11.By comparing the results of Figs. 10 and 11 it is deduced that the drop sizes are reduced as the gas pressure drops from 300 to 150 psia.This result may seem at odds with the fact that smaller drop sizes are produced at higher gas pressures (or gas density) due to the enhancement of gas-liquid aerodynamic interaction as was discussed in relation with Fig. 9.However, it should be kept in mind that in order to maintain the aerodynamic Weber number range the same in both Figs. 10 and 11, the gas-liquid relative velocity must be increased at higher gas pressures.Therefore, for the same aerodynamic Weber number, the aerodynamic interactions are greater at smaller gas pressures.
It is evident from Figs. 10 and 11 that the present model yield predictions that are in relatively favorable agreement with the experimental data particularly at gas pressure 300 psia and for the lower range of Weber number.The agreement between measurements and simulations is best for satellite drop sizes.The experimental trend that the drop sizes decrease with the aerodynamic Weber number is also mimicked by the numerical simulations.However, the model accuracy suffers at higher aerodynamic Weber number and at a gas pressure of 150 psia, as can be seen in Fig. 11, due to the adoption of a linearized aerodynamic term in Eq. ( 8).The linearized aerodynamic term is incapable of rendering accurate representation of strong aerodynamic interactions present at higher gas density which are highly nonlinear.Therefore, implementation of a nonlinear aerodynamic term is planned to improve the model's performance.

Conclusions
A well coordinated experimental and computational study of the atomization of shear coaxial liquid jets has been presented.The experimental technique involves innovative digital imaging processes.The experimental data cover a wide range of parameters relevant to realistic SSME operating conditions.A rigorous mathematical formulation and numerical solution scheme of the computational model equations are demonstrated.
The computational model appropriately captures the effects of capillary, viscous and aerodynamic forces.
Increasing viscous forces dampens the instability of the jet.A larger gas to liquid density ratio boosts jet atomization due to a more pronounced aerodynamic interaction.Reasonably good predictions of the measured drop sizes are obtained especially at lower gas-liquid relative velocity.
The present model represents a first cut-effort at developing a computationally effective and robust model to simulate the atomization processes in shear coaxial liquid jets.The model could benefit from incorporating a nonlinear aerodynamic term to replace the current linearized version.This maybe achieved via applying the source panel method where the jet surface is approximated by panels with associated source strengths (Fletcher, 1991, Yong & Reitz, 2004).In addition, a boundary condition model for a pinched-off end must be added to permit following of the secondary breakup of the ligaments produced at low wave number and to examine the behavior of longer ligaments subject to more complex disturbance functions.Further validation of the model against experimental data and empirical relations is also intended.

FigureFigure 5 .Figure 7 .
Figure 1.Schematic of WNIST facility Three parameters were chosen to determine the injector's flow range to be imaged: Reynold's number of the Note that the experimental Reynolds number, Re l , is based on the jet diameter and liquid velocity, U l , not jet radius and relative velocity, U 0 , as Reynolds number, Re , that appears in Eq. (11) of the present model.Since the gas-liquid relative velocity not the liquid velocity appears in the model equations, it was only possible to match Weber numbers when comparing numerical and experimental results.Test ranges consisted of holding one of these three values constant, while varying the other two numbers at determined intervals.Table 1 below shows the test ranges for the values of Re l , =.