Simulations of highly under-expanded jets

Simulations of highly under-expanded jets R. Buttay, P.J. Mart´ınez Ferrer, G. Lehnasch and A. Mura Institut Pprime UPR 3346 CNRS, ISAE - ENSMA and University of Poitiers, FRANCE 1. Introduction Highly resolved numerical simulations of underexpanded jets are conducted. Such high speed jets may result from the accidental release of high pressure flammable mixtures into the quiescent atmosphere, which poses important concerns related to explosion hazards. The nature of the hazard will depend on the stability of any jet fire resulting from the under-expanded release of fuel. More precisely, it will depend on whether or not combustion can be sustained in the vicinity of the release or there is a delay during which an explosive cloud may form. The description of such under-expanded jets covers also a broad range of applications related to spacecraft propulsion including hypervelocity Scramjets or rocket engines. For instance, an underexpanded torch jet is used to initiate combustion in expander cycle engines. The description of scalar mixing downstream of the Mach bottle thus appears as an essential issue, which is central to the present paper. It constitutes a preliminary step before a more detailed analysis of the effects of heat release, chemical kinetics and self-ignition on such compressible jet structures. 2. Numerical methods The present study is carried out with a numerical solver able to describe compressible multicomponent reactive mixtures. We therefore consider the compressible Navier-Stokes equations written for a reactive multicomponent mixture. The treatment of the inviscid component of the transport equation for the conservative vector relies on the seventh-order accurate Weighted Essentially Non-Oscillatory (WENO7) reconstruction of the characteristic fluxes. In practice, the numerical solver uses a seventh-order accurate centered finite difference scheme, and the application of the WENO7 scheme is conditioned to a smoothness criterion which involves the local values of the normalized spatial variations of both pressure and density. The viscous and molecular diffusion flux functions are determined thanks to an eighthorder centered difference scheme. The temporal integration is performed by using an explicit thirdorder TVD Runge-Kutta algorithm. Further details about the numerical methods as well as an exhaustive verification of the solver are provided by Martinez Ferrer et al (2014). 3. Numerical setup The studied simulation consists in air released from a high pressure vessel into the quiescent atmosphere. Air is considered as a two-species mixture (O2and N2) described with variable heat capacities and transport properties thus avoiding the resort to simplifying hypotheses such as constant heat capacity ratio value. The corresponding release velocity is 630 m/s at the exit (Ma = 1). The flow field is initialized with the mixture characteristic of air at a pressure of 1 atm and a temperature of 300 K. The diameter of the injector exit is set to D= 0.001 m, which corresponds to a Reynolds number of 77500. The inflow parameters retained to perform the present simulation are listed in table 1 and correspond to a sonic under-expanded jet with a nozzle pressure ratio (NPR) based on static pressure of fifteen. Table 1. Under-expanded jet flow parameters. Injector Free-stream P(atm) 15.0 1.0 T(K) 1000.0 300.0 M a 1.0 0.05 u(m/s) 630.0 20.0 YO20.233 0.233 YN20.767 0.767 The computational domain dimension, nondimensionalized by the diameter of injection D, are L∗ x1 = 14 and L∗ x2 =L∗ x3 = 6. This domain is discretized with Nx1×Nx2×Nx3= 880×449×449 nodes Cartesian grid, which corresponds to approximately 180 millions nodes. Sponge regions combining both grid coarsening and explicit filtering are used in order to avoid spurious numerical wave reflections and to make easier the processing of open boundary conditions. The resolution in the highly resolved region is Δx= Δy= Δz= D/60. Perfectly non-reflecting boundaries conditions are applied at the top, bottom, backside and frontside boundaries. Partially non-reflecting boundary condition is imposed at the outflow. The value of the CFL number is set to 0.75 and the Fourier number Fo is set to 0.1. y/D Figure 1. Numerical schlieren image of the near-field of the jet. 21st Intl. Shock Interact. Symp. 189 3 - 8 Aug. 2014, Riga, Latvia 4. Description of the velocity field The structure of the axisymmetric free jet expanding through the small orifice into quiescent atmosphere is displayed in Fig. 1. The whole compressible structure is clearly identified. As the flow leaves the nozzle the high pressure mismatch causes it to expand and accelerate (Fig. 2). Expansion waves originate near the expansion point, propagate and meet the outer boundary of the jet, where they are reflected as compression waves. Coalescence of these pressure waves results in a curved barrel shock surrounding the immediate supersonic region. The reflection of the incident shock is not regular and a Mach disk pattern appears in the near-field of the jet. The flow is subsonic just behind the Mach disk, while it remains supersonic downstream of the barrel shock. The triple point connects various discontinuities and becomes the origin of a new slip line, which gives rise to a supersonic shear layer. 0 0.5 1 1.5 2 2.5 3 3.5 4 4.5 5 5.5 6 0 2 4 6 8 10 12 14 P* Ma x/D Figure 2. Axial profiles of Mach number and normalized pressure P∗=P/Pa. 0 1 2 3 0 2 4 6 8 10 12 14 Present study Lehnasch et al. (2005) Velikorodny et al. (2012) Yuceil et al. (2003) Chauveau et al. (2006) u∗ x/D Figure 3. Normalized axial velocity profiles u∗=u/ue along the centerline of the jet. Figure 2 shows centerline values of the normalised pressure P∗=P/Paand Mach number. Due to the expansion of the gas, the pressure and the temperature decrease significantly while the velocity and the Mach number increase (Fig. 2 and Fig. 3). The Mach number drop indeed allows to delineate the Mach disk location at approximately x/D = 3.55. The position of the Mach disk is checked against the empirical correlation of Ashkenas et al (1966) xDM /D = 0.67�P0/Pa, which provides a similar estimate: xDM /D = 3.60. Downstream of the Mach disk, the pressure stabilizes around the atmospheric pressure while the flow reaccelerates progressively and becomes supersonic at a distance x/D ≈13 from the injector. Figure 3 reports comparisons of streamwise velocity centerline values with previous experimental and numerical data. Considering the difficulties associated with high velocity measurements above x/D = 2, the present velocity profile displays a satisfactory level of agreement with experimental data. 5. Scalar field: turbulent mixing We now investigate scalar mixing in such highly compressible flow. To this purpose, we solve an additional transport equation for a passive scalar ξwich is advected with an unity Lewis number. The molecular diffusivity of ξis taken equal to the thermal diffusivity and Lewis number effects associated to differential diffusion are thus dropped off from the present analysis. This tracer ξis defined to be unity in the jet and zero elsewhere. Figure 4 presents the mean Mach number field superimposed with four iso-contours of ξ(0.1,0.4,0.7,0.95). 0 5.3 x/D Figure 4. Mean Mach number field superimposed with ξiso-contours. 0 1 x/D Figure 5. Mean ξfield superimposed with ξisocontours. 21st Intl. Shock Interact. Symp. 190 3 - 8 Aug. 2014, Riga, Latvia Figure 5 reports the mean field of the tracer ξ. The boundary of the mean jet behaves like an impermeable membrane. In this configuration, the turbulence develops at this boundary and grows in the supersonic shear layers. Figure 6 and 7 display the radial profiles of the mean value and variance of the passive scalar at different crosssections of the jet. The quantity � ξ��2characterizes the dispersion of the scalar from its average value. 0 0.2 0.4 0.6 0.8 1 0 0.5 1 1.5 2 2.5 3 x/D = 6 x/D = 8 x/D = 10 x/D = 12 x/D = 14 � ξ r/D Figure 6. Radial profiles of � ξ. 0 0.005 0.01 0.015 0.02 0.025 0.03 0.035 0.04 0 0.5 1 1.5 2 2.5 3 x/D = 6 x/D = 8 x/D = 10 x/D = 12 x/D = 14 � ξ��2 r/D Figure 7. Radial profiles of � ξ��2. Production of variance reflects the inhomogeneity of the local mixture. Conversely, its destruction characterizes the action of molecular processes trough the mean value of the SDR (scalar dissipation rate) χξ=D(∂ξ/∂xk)(∂ξ/∂xk) of the passive tracer. The scalar mixing takes place in the supersonic shears layers as shown in Fig. 6. One can notice that at x/D = 12, the value of � ξon the symmetry axis (r/D = 0) is lower than unity. The shear layers indeed intersect the centerline of the jet at a distance of x/D ≈12 which corresponds approximately to the length of the subsonic throat. Figure 7 illustrates the mixture homogenization and the destruction of variance taking place between the plane x/D = 6 and the plane x/D = 14. The variance � ξ��2is also plotted against the mean value � ξin Fig. 8. Thus, the resulting profiles can be compared to the maximum realizable value of the scalar variance, as given by � ξ(1 −� ξ). Figure 9 and 10 represent respectively ξ�2and χξalong ξiso-contours defined previously. For x/D ≤3, profiles are less representative due of the lack of point in this region to capture the dynamics which may explain the persistence of residual oscillations on the profiles. The peak of both ξ�2and χξat x/D ≈4 perceptible on the iso-contours ξ= 0.4, 0.7 and 0.95 are explained by the impact of the reflected shock wave while the iso-contour ξ=0.1 does not cross the reflected shock wave and does not exhibit the presence of such a peak. This is consistent with the previous investigation of Huh et al (1996) who pointed out the shock wave enhancement of mixing by deflecting streamlines in mixing region and generating shock-generated vorticity. 0 0.005 0.01 0.015 0.02 0.025 0.03 0.035 0.04 0 0.2 0.4 0.6 0.8 1 x/D = 6 x/D = 8 x/D = 10 x/D = 12 x/D = 14 � ξ��2 � ξ Figure 8. Scalar variance plotted versus � ξ. 0 0.005 0.01 0.015 0.02 0.025 0.03 0.035 0.04 0.045 0 2 4 6 8 10 12 14 ξ= 0.1 ξ= 0.4 ξ= 0.7 ξ= 0.95 - - - - ξ�2 x/D Figure 9. Profiles of ξ�2along differents ξiso-contours. Probability density functions (PDFs) of ξare 21st Intl. Shock Interact. Symp. 191 3 - 8 Aug. 2014, Riga, Latvia reported on figure 11. These PDFs are evaluated on iso-contours of ξfor different values of x/D. The shape of the PDF, according to their position in the jet, are consistent with their theoretical counterparts (Bilger (1980)). On the boundary of the jet and internal side of the jet, i.e. ξ= 0.1 and ξ= 0.95, the influence of molecular mixing effects on ξis less appreciable. In contrast, inside the shear layer, i.e. ξ= 0.4 and ξ= 0.7, the PDF shapes tighten around the mean value. Scalar dissipation rate is known to play a crucial role for non-premixed conditions since mixing is a prerequisite before combustion occurs. In the field of turbulent combustion modeling approach, it remains a common practice to close the average SDR that appears in the RHS of the scalar variance transport equation (1) (see appendix) by invoking a similarity hypothesis between scalar and velocity turbulence spectra which results in the classical approximation τξ�Cξτtwith Cξa modeling constant. This leads to the Linear relaxation model (LRM) which consists in �εξ=� ξ��2/τξ� � ξ��2/Cξτt. 0 1000 2000 3000 4000 5000 6000 7000 0 2 4 6 8 10 12 14 ξ= 0.1 ξ= 0.4 ξ= 0.7 ξ= 0.95 - - - - χξ x/D Figure 10. Profiles of χξalong differents ξisocontours. The validity of the simplified LRM closure is analysed by investigating the proportionality constant between τξand τt, i.e. the scalar to turbulence time scale ratio Cξ=τξ/τtin the present configuration. Figure 12 reports the scalar mixing time scale defined as τξ=� ξ��2/�εξalong the isocontours of the mean value ξ. Figure 13 reports the turbulence time scale defined as τt=k/ε along ξiso-contours. The quantity kis the turbulent kinetic energy while εdenotes its dissipation rate. Finally, Fig. 14 reports the profile of the scalar to turbulence time scale ratio. In this figure, it is noteworthy that the shock-wave impact on the mixing layer does not significantly impact the scalar to turbulence time scale ratio. Moreover, it is also remarkable that, in highly compressible situations such as those considered therein, the hypothesis of a constant value, as implied by the LRM closure, is rather well verified and not significantly altered by compressibility effects. 0 1 2 3 4 5 6 7 0 0.2 0.4 0.6 0.8 1 x/D = 6 x/D = 8 x/D = 10 x/D = 12 x/D = 14 P(ξ) ξ (a) ξ= 0.1 0 0.5 1 1.5 2 2.5 3 3.5 0 0.2 0.4 0.6 0.8 1 x/D = 6 x/D = 8 x/D = 10 x/D = 12 x/D = 14 P(ξ) ξ (b) ξ= 0.4 0 0.5 1 1.5 2 2.5 3 3.5 4 4.5 0 0.2 0.4 0.6 0.8 1 x/D = 6 x/D = 8 x/D = 10 x/D = 12 x/D = 14 P(ξ) ξ (c) ξ= 0.7 0 10 20 30 40 50 60 70 80 0 0.2 0.4 0.6 0.8 1 x/D = 6 x/D = 8 x/D = 10 x/D = 12 x/D = 14 P(ξ) ξ (d) ξ= 0.95 Figure 11. Pdf of ξon ξ= 0.1 (a), ξ= 0.4 (b), ξ= 0.7 (c) and ξ= 0.95 (d) iso-contours. 21st Intl. Shock Interact. Symp. 192 3 - 8 Aug. 2014, Riga, Latvia 0 5e-06 1e-05 1.5e-05 2e-05 2.5e-05 3e-05 3.5e-05 0 2 4 6 8 10 12 14 ξ= 0.1 ξ= 0.4 ξ= 0.7 ξ= 0.95 - - - - τξ x/D Figure 12. Profiles of τξalong differents ξisocontours. 0 5e-06 1e-05 1.5e-05 2e-05 2.5e-05 3e-05 3.5e-05 0 2 4 6 8 10 12 14 ξ= 0.1 ξ= 0.4 ξ= 0.7 ξ= 0.95 - - - - τt x/D Figure 13. Profiles of τtalong differents ξiso-contours. 0 0.2 0.4 0.6 0.8 1 1.2 0 2 4 6 8 10 12 14 ξ= 0.4 ξ= 0.7 - - Cξ x/D Figure 14. Profiles of Cξ=τξ/τtalong differents ξ iso-contours. 6. Conclusions Highly resolved numerical simulations of highly under-expanded turbulent gas jets have been conducted. The comparisons performed between the present results and experimental data sets or empirical correlations give rise to a satisfactory level of agreement. Special emphasis has been placed on the description of turbulent mixing downstream of the Mach disk structure. The study has been focused on the applicability of the linear relaxation model (LRM) as a possible closure of the mean scalar dissipation rate and especially on the mapping of the scalar to turbulence time scale ratio Cξ. The obtained results show that the hypothesis of a constant value, as implied by the LRM closure, is rather well verified and not significantly altered by compressibility effects. Future works will consist in computing highly under-expanded hydrogen/air jet so as to determine the flammability index (FI) map in such conditions. This will provide very valuable insights into security issues relevant to hydrogen explosion hazards Acknowledgements The present work is part of the PhD thesis of Romain Buttay, financially supported by CNRS and Region Poitou-Charentes. This work was granted access to the HPC resources of IDRIS under the allocations x20142a0912 and x20142b7251 made by GENCI (Grand Equipement National de Calcul Intensif). Appendix The transport equation for the scalar variance writes : ∂ ∂t(ρξ��2) + ∂F ξ��2 k ∂xk =−2ρD ∂ξ�� ∂xk ∂ξ�� ∂xk −2ρu�� kξ�� ∂� ξ ∂xk (1) with the scalar flux Fξ��2 k= (ρukξ��2−ρD ∂ξ��2 ∂xk ). In this transport equation, the first term in the left hand side is the accumulation term and the second is the (conservative) flux term (convection and diffusion). On the right hand side of eq. (1), the first term corresponds to mean SDR while the second is the production associated to mean concentration gradients. References Ashkenas H., Sherman F.S. (1966), Structure and utilization of supersonic free jets in low density wind tunnels, Rarefied Gas Dynamics, 2 84– 105. Beguier C., Deskeyser L., Launder B.E. (1986), Ratio of scalar and velocity dissipation time scales in shear flow turbulence, Physics of Fluids A 21, 307310. Bilger R. (1980), ”Turbulent flows with nonpremixed reactants”, in Turbulent Reacting Flows, Topics in Applied physics. 21st Intl. Shock Interact. Symp. 193 3 - 8 Aug. 2014, Riga, Latvia Chauveau C., Davidenko D.M., Sarh B., Gokalp I., Avrashkov V., Fabre C. (2006), PIV measurements in an under-expanded hot free jet, 10th International Symposium on Application of Laser Techniques to Fluid Mechanics, Lisbon, Portugal, 26–29 june 2006. Ewan B.C.R., Moodie K. (1986), Structure and velocity measurements in under-expanded jets, Combustion Science and Technology, 45 275. Gomet L., Robin V., Mura A. (2012), Influence of residence and scalar mixing time scales in non-premixed combustion in supersonic turbulent flows, Combustion Science and Technology, 184:10-11, 1471–1501. Huh H., Driscoll J.F. (1996), Shock wave enhancement of the mixing and the stability limits of supersonic hydrogen-air jet flames, Proceedings of the 26th International Symposium on Combustion, 1996, 2933–2939. Lehnasch G. (2005), Contribution a l’etude numerique des jets supersoniques sous-detendus, PhD Thesis, University of Poitiers, 2005. Martinez Ferrer P.J., Buttay R., Lehnasch G., Mura A. (2014), A detailed verification procedure for compressible reactive multicomponent Navier-Stokes solver, Computers & Fluids, 89 88–110. Velikorodny A., Kudriakov S. (2012), Numerical study of the near-field of highly underexpanded turbulent gas jet, International Journal of Hydrogen Energy, 37 17390–17399. Yuceil K.B., Otugen M.V., Arik E. (2003), Interferometric Rayleigh scattering and PIV measurments in the near-field of under-expanded sonic jets, 41st Aerospace Science Meeting and Exhibit, Reno, Nevada, 6–9 January 2003. 21st Intl. Shock Interact. Symp. 194 3 - 8 Aug. 2014, Riga, Latvia