Numerical Simulation of Jet Injection and Species Mixing under High-Pressure Conditions

Direct Numerical Simulation (DNS) and Large Eddy Simulation (LES) are performed of a round fluid jet entering a high-pressure chamber. The chemical compositions and temperatures of the jet and that of the fluid in the chamber are initially prescribed. The governing equations consist of the conservation equations for mass, momentum, species and energy, and are complemented by a real-gas equation of state. The fluxes of species and heat are written in the framework of fluctuation-dissipation theory and include Soret and Dufour effects. For more than two species, the full mass diffusion and thermal diffusion matrices are computed using high-pressure mixing rules which utilize as building blocks the corresponding binary diffusion coefficients. The mixture viscosity and thermal conductivity are computed using standard mixing rules and corresponding states theory. To evaluate the physical model and numerical method, LES is employed first to simulate a supercritical N2 jet injected into N2. Time averaged results show reasonable agreement with the experimental data. Then, DNS is conducted to study the spatial evolution of a supercritical N2 jet injected into CO2. Time averaged results are used to compute the length of the potential core and the species diffusion characteristics. Spectral analysis is then applied on a time series data obtained at several axial locations and a dominant frequency is observed inside the potential core.


Introduction
Numerous propulsion devices -diesel, rocket, gas turbine and HCCI engines -rely on mixing and combustion of chemical species under high-pressure conditions. What distinguishes highpressure (high-p) ‡uid behavior from its low-pressure counterpart is the nature of the molecular interactions; under high-p the molecules are much tightly packed and may experience van der Waals interactions and caging e¤ects. Macroscopically, this means that the perfect gas equation of state (EOS) is no longer appropriate to describe this ‡uid and should be replaced by a real-gas EOS. Moreover, under high-p conditions the ‡uxes of mass and heat must assume their full form given by the ‡uctuation-dissipation theory [1], that includes Soret and Dufour e¤ects. For multispecies situations the mass di¤usion and thermal di¤usion matrices must be fully populated to account for interactions among species rather than interactions between species or interactions of one species with its global complement in the mixture. Clearly, simulating such ‡ows requires complex models. Evaluating such models and simulations with experimental data is crucial in order to gain con…dence that the model captures the essential aspects of the situation and that the numerical method is appropriate. Despite the de…nite interest in understanding multispecies situations, experimental data for only single and binary species exist so far. Selected representative experimental studies are those of Oschwald and Schik [2] and Meyer et al. [3] who obtained data for N 2 -into-N 2 injection; Oschwald et al. [4] who performed experiments with N 2 injection into H 2 ; Chehroudi et al. [5] who injected N 2 into N 2 , O 2 , He or a 50%/50% mixture of CO and N 2 ; Segal and Polikhov [6], Roy and Segal [7] and Roy et al. [8] all of whom obtained data for a ‡uoroketone jet injected into N 2 ; and Falgout et al. [9,10] who observed the injection of n-dodecane in air containing 21% O 2 , and butanol, n-dodecane and n-hexadecane in N 2 , respectively. No experiments with multi-species (i.e. more than two species) mixing exist, except those of Manin et al. [11] for n-heptane, n-dodecane and n-hexadecane injected into a hot ‡uid mixture resulting from the burning of C 2 H 2 and H 2 in air, but the percent composition of this mixture is not documented although it can generally be assumed to be of H 2 O, CO 2 and N 2 in unknown proportions. Most experimental observations are of visual type trying to unravel the species-speci…c conditions separating the single-phase and two-phase regimes. Quantitative data for validation are rare and restricted to few studies. We ask here whether the available data is su¢ cient to validate models, i.e. distinguish between the predictions of models that are di¤erent and indicate the model which best represents the physics. Failing to do would mean that experimental data of di¤erent type than one acquired so far must be obtained, and would open up the discussion regarding which data would be appropriate for model evaluation and whether the diagnostics to acquire this data in an incontrovertible way are available. Among the above-cited investigations, we are particularly interested in the data of Mayer et al. [3] for N 2 -into-N 2 injection and of forthcoming data for N 2 injected into CO 2 . This last case is partially relevant to engines using exhaust gas re-circulation and of paramount importance for understanding conditions encountered by a lander on Venus. This paper is organized as follows: First, the previously derived governing equations in the most general form [12] are described, and we refer to other publications for the detailed computation of transport properties. Then, we conduct both LES to evaluate the model by comparison with experimental data and DNS to understand the details of ‡ow evolution in a spatial con…guration rather than the temporal con…guration of the past [12]. Model evaluation must be necessarily conducted in a spatial rather than temporal realm. Since all our previous studies (e.g. [12]) were performed in the context of a temporal mixing layer, the con…guration, numerical method and boundary conditions including in ‡ow conditions are described in detail. The description and analysis of the results follows. A summary and conclusions are …nally o¤ered.

Description of the governing equations
The governing equations are the conservation equations -continuity, momentum, species and total energy -complemented by the Peng-Robinson (PR) real-gas equation of state (EOS). In the conservation equations, the ‡uxes of species and heat are written according to ‡uctuationdissipation theory [1]. The species ‡uxes consist of the multicomponent species mass di¤usion with a fully populated di¤usion matrix and the Soret e¤ect; the heat ‡ux contains the Fourier term, Dufour e¤ects and the enthalpy transported by the species ‡uxes. All transport properties -viscosity, mass di¤usion coe¢ cients, thermal di¤usion coe¢ cients and thermal conductivityare computed using high-pressure-valid mixing rules [13], and thus transport properties are functions of all thermodynamic variables, (p; T; X ) where T is the temperature and X is the mole fraction of species : The PR EOS is well known to be inaccurate for species other than hydrocarbons and therefore the volume correction of [14] computed from a reference Gibbs free energy is used to enhance its accuracy to the same level as the Lee-Kesler EOS. Details on the entire formulation can be found in [12].
According to the 'alternate approach' suggested by Jordan [15], to obtain the present LES Table 1. Initial conditions and other characteristics of the simulations. The abbreviation "comp." denotes the chemical composition.
conservation equations, the governing di¤erential equations are …rst transformed into curvilinear coordinates and then …ltered with a lowpass …lter rather than opposite (i.e. the …ltering occurs before the transformation); the two approaches lead to the same mathematical equations. This …ltering results in the formation of unclosed terms in the conservation equations that contain the subgrid activity; these are the subgrid-scale (SGS) terms. The LES equations have been presented in general form in [16], and contain the conventional subgrid terms as well as unconventional subgrid terms arising from the steep thermodynamic-variable gradients which are inherent in the ‡ow ( [12]). Since LES is only here used in the context of single ‡uid composition, the unconventional SGS terms identi…ed for a mixture of species initially at di¤erent temperatures and the unconventional SGS terms arising from the species ‡uxes are neglected. Similarly, since the available data is only global, that is, the spatial details of the ‡ow cannot be validated, we also neglect, in a …rst approximation, the SGS terms arising from the pressure gradient in the momentum equation and the heat ‡ux in the energy equation that were shown to be important to recover in LES the spatial ‡ow details [17,18]. The constant-coe¢ cient Smagorinsky model is used for the SGS conventional ‡uxes since it has been shown that if the interest is only in the mean ‡ow, the nature of the SGS model is not as important as the thermodynamics [19].

Con…guration, numerical method and boundary conditions
The con…guration of interest is that of a three-dimensional round jet injected into a chamber at high-p conditions. The chamber dimensions for DNS in the (x 1 ; x 2 ; x 3 ) directions are 25D, 10D and 10D respectively where D is the diameter of the jet, here set to be 0.005m. The domain extent for LES are 50D, 25D and 25D in the (x 1 ; x 2 ; x 3 ) directions respectively. The domain is large enough to assume that the injected ‡uid is e¤ectively a free jet. The boundary conditions in the x 2 and x 3 directions are of non-re ‡ecting type for real-gases [20]. In the x 1 axial direction, the boundary conditions are in ‡ow and out ‡ow type. Acoustic re ‡ections due to non-orthogonal wave impingement on the boundaries is prevented by applying sponge boundary conditions [21] and by arti…cially increasing viscosity near the out ‡ow boundary; this region is not used when presenting statistical averaging results. The jet in ‡ow velocity pro…le is given by where U j is the speci…ed jet centerline velocity, r = p x 2 2 + x 2 3 , r j is the jet radius and is the initial shear layer thickness chosen to be 0.05r j . Similar hyperbolic tangent pro…les are used at the in ‡ow for the temperature and the species mass fraction. For the DNS, the domain is initialized using a method previously adopted by Nichols et al. [22] for simulating variable density jets at atmospheric pressure. In this method, the same pro…le as that of the in ‡ow (i.e. x 1 = 0) for velocity, temperature and mass fraction is initially speci…ed at all x 1 locations throughout the domain. Random perturbations of amplitude 5% of U j , in the form of white noise are added to the initial velocity …eld in the shear layer. After this initial forcing, the ‡ow evolves with no further perturbations being added. For the LES, the in ‡ow perturbations are provided using the method adopted by Freund [23].
The governing equations are temporally discretized using a fourth-order explicit Runge-Kutta time integration and spatially discretized utilizing a sixth-order compact scheme. In both DNS and LES, stability was achieved by …ltering the conservative variables using a tenth-order …lter [24]. The …ltering was performed every 4 time steps in DNS and every time step in LES; the rationale for these …ltering frequencies was explained in detail elsewhere [17].

Results
4.1. LES of N 2 =N 2 jet LES of N 2 supercritical jets injected into a chamber containing N 2 were performed in the past. Zong et al. [25] performed two-dimensional jet simulations and compared the radial variation of the density pro…le at two axial locations to the experimental data of Chehroudi et al. [26]; the agreement with the experimental data was reasonable despite the lack of portraying stretching/tilting e¤ects which are typically considered important in vorticity production. Schmitt et al. [27] performed three-dimensional LES of supercritical jets using a Wall Adapted Local Eddy viscosity(WALE) model and the centerline density pro…le was compared with the experimental data of Mayer et al. [3] with good agreement. More recently Müller et al. [19] performed LES of Mayer's experiments using the Smagorinsky and Vreman models and found good agreement for the …rst order statistics between the models and a good agreement with the experiment for the centerline density pro…le. Considering that di¤erent models led to good agreement with data indicates that the data so far available is insu¢ cient to distinguish between, or among, models -a model based on appropriate and relevant physics versus a model containing inappropriate simpli…cations -and therefore none of the favorable comparisons with available data can be considered a validation of that model. Thus, we propose that good agreement with existing data is only a necessary, but not su¢ cient condition, for model validation. In fact, we propose that one of the roles of simulations should be to highlight quantities which would constitute discriminating tests able to select the best model to describe high-pressure mixing and combustion.
With the goal of showing the necessary evaluation with the data, LES of a supercritical N 2 jet into a N 2 -…lled chamber is performed and compared to the experiments of Mayer et al. [3].   Table 1. Figure 1(a) shows the x 1 variation of the time averaged density along the jet centerline obtained from experiments and simulation. The time averaged quantities (denoted by <> hereafter) is de…ned over a …xed time interval as where t 0 = t 2 t 1 is the integration time. The jet centerline is de…ned as the line along which x 2 = 0 and x 3 = 0, and the values along this line are computed by interpolating from the nodes closest to the line at all x 1 locations. The simulation results agree reasonably with the experimental results. We also compare our results with those of Schmitt et al. [27] and a good agreement is obtained between the two simulations that use very di¤erent numerical schemes and turbulence models. This corroborates the observation of Müller et al. [19] that the …rst order statistics are generally insensitive to the numerical schemes and turbulence models and further reiterates the need for experimental data pertaining to higher order statistics. Figure  1(b) depicts the instantaneous (t = 0:0416 s) temperature contours in the x 3 = 0 plane. The presence of small scales downstream of the potential core is evident. Stretched entities emanating from the potential core can be observed that are characteristic of a single-phase mixing layer observed at the edge of supercritical jets.
4.2. DNS of N 2 /CO 2 jet DNS of a supercritical N 2 jet injected into a chamber containing CO 2 is performed at Re j =1000.
To our knowledge, this is the …rst DNS of a jet having a speci…ed composition that is injected into a chamber having a di¤erent composition. The jet exit and chamber conditions are listed in Table 1. Figure 2 shows the isocontours of the Q-criterion which is a commonly used quantity to visualize vortical structures, where ! i is the ith component of the vorticity vector and S ij is the strain-rate tensor. The isocontours show the vortex ring near the in ‡ow that further breaks up to produce a range of small scale vortices. Figure 2 shows the jr j contours in the x 1 x 2 plane at x 3 = 0. These high-value jr j regions are a signi…cant feature of high-p ‡ows and representing them accurately is an important criterion measuring the success of subgrid scale models in the context of high-p ‡ows. As a manifestation of mixing due to molecular di¤usion enhanced by turbulence, shown on the x 1 x 3 plane (at x 2 = 0) are the Y N 2 contours, hY N 2 i ; where Y is the mass fraction. Since the injected N 2 ‡uid is lighter than the surrounding, the encountering of the heavier ‡uid leads to a rapid three-dimensional breakdown of the jet and therefore to enhanced mixing. The edges of mass fraction contours also show evidence of a single-phase mixing layer which is characteristic of supercritical jets, as opposed to spherical droplets formed due to surface tension at subcritical pressures. Quantifying this aspect requires further investigation.
To examine if a potential core exists, statistics are obtained by performing time averaging for approximately 100tU j =D. Figure 3(a) shows the variation of h i = j along the jet center line. A time averaged potential core length x 1 =D = 5 is obtained after which rapid mixing of the jet is observed. Figure 3(b) shows hY N 2 i which also exhibits a potential core length that extends 5 diameters downstream of the injector; the decreased hY N 2 i past x 1 =D = 5 is due to mixing with the heavier CO 2 that is responsible for the increase of h i = j past that axial location. To further obtain an understanding of the ‡ow, we investigate time-averaged statistics along the radial direction at four di¤erent streamwise locations x 1 =D = 2, 10, 15 and 20. Figure 4(a) illustrates the variation of the mean axial velocity, non-dimensionalized by the jet centerline velocity, along the radial direction at di¤erent streamwise locations. The x 1 =D = 2 location is inside the potential core and hence the pro…le retains a shape similar to the inlet velocity pro…le. However, at further downstream locations, rapid mixing causes the pro…le to ‡atten out as the jet ‡uid mixes with the ambient ‡uid and the jet spreads out. Figure 4(b) illustrates the variation of hY N 2 i along the radial direction. By comparing the mean velocity and mean mass fraction pro…les it can be observed that the decay rate of the jet centerline value is larger for the species than that of axial velocity; however the spreading rate in the radial direction  remains approximately the same. Whether this axial evolution and radial spreading comparison are universal, and if so, the reasons for these e¤ects should be investigated in future work.
At atmospheric p, it is generally considered that mixing occurs primarily due to molecular mass di¤usion enhanced by turbulence. Under high-pressure conditions, the Soret e¤ect may also play an important role in species di¤usion. The species mass ‡ux is given by where D T; ; D p; and D 0 are di¤usion coe¢ cients and m is the molar mass (see [12] for details). The contribution of the Soret term (the …rst term in Eqn. 4) along with the other two terms are plotted in Figure 5. Since a statistical representation of these terms would involve computing budgets, a calculation which is beyond the scope of this investigation, we defer that calculation to future study and present here only the magnitude of the instantaneous values on x 3 = 0 plane. The largest contribution to J comes from the mass fraction gradient term ( Figure 5(c)) followed by the Soret term ( Figure 5(a)) the magnitude of which is approximately two orders of magnitude smaller than the largest term. The contribution of rp ( Figure 5(b)) is negligible for the current case, as generally accepted. By comparing the spatial distribution of these terms with that of ! 3 contours in Figure 5(d), we can conjecture that these di¤usion e¤ects play an important role in the formation of small scales through mixing induced formation of large magnitude jr j regions which act as a solid boundary at which vorticity forms. A detailed evaluation is necessary to further shed light on this aspect. To understand the turbulent characteristics of the ‡ow, a spectral analysis is performed by applying a Fourier transform to , Y N 2 and u 1 at four locations (x 1 =D = 4; 8; 12 and 16) along the jet centerline. A dominant frequency at a Strouhal number (St D = f D=U j , where f is the frequency (f = 1=t)) of 0.33 is obtained. This is the frequency at which a continuous series of vortex rings are generated and is particularly visible as peaks in the spectra at x 1 =D = 4 in Figure 6(a). As x 1 =D increases, (e.g. x 1 =D = 8; 12 and 16) in Figures 6(b), (c) and (d), the spectra becomes smoother, a fact which is manifestation of a ‡ow with turbulent characteristics. The spectrum shows a monotonic fall in the energy at higher frequencies indicating that all relevant frequencies are captured in the simulation.

Summary and Conclusions
Direct Numerical and Large Eddy Simulations are used to study high-pressure jet injection and mixing. LES is …rst used to evaluate the physical model and numerical method for a cryogenic N 2 jet injected into a N 2 -…lled chamber. The time averaged density variation along the streamwise direction is compared to the experimental data and a reasonable agreement is obtained therefore ful…lling a necessary condition for suggesting the validity of the physical model and numerical method. The …rst DNS of a jet of a speci…ed chemical composition injected into a chamber …lled with a ‡uid of di¤erent composition is conducted: the DNS is of a supercritical N 2 jet injected into a CO 2 -…lled chamber. Instantaneous results show that the employed grid and model are able to capture the high density gradient magnitude regions that are characteristic of a high-p ‡ow. The radial variation of mean axial velocity and mean N 2 mass fraction indicate that the centerline decay rate of species concentration is larger than the axial velocity decay rate; whether this fact is universal and its causes will be the subject of future research. The contribution of the Soret term is approximately two orders of magnitude smaller than the mass fraction gradient term; the spatial distribution of all di¤usion terms indicate their important role in the formation of smaller turbulent scales. Statistical results show that the mean core length is about 5D. Spectral analysis shows a dominant frequency of St D = 0:33, the frequency at which vortex rings are produced. Spectra from locations downstream of the potential core show a broadband nature, characteristic of a turbulent ‡ow.