- Email: [email protected]
Numerical study of hydrogen–air mixing in turbulent compressible coaxial jets
i n t e r n a t i o n a l j o u r n a l o f h y d r o g e n e n e r g y 4 0 ( 2 0 1 5 ) 9 5 3 9 e9 5 5 4
Available online at www.sciencedirect.com
ScienceDirect journal homepage: www.elsevier.com/locate/he
Numerical study of hydrogeneair mixing in turbulent compressible coaxial jets R. Ouzani*, M. Si-Ameur LESEI Laboratory, Mechanical Engineering Department, University of Batna, Algeria
article info
abstract
Article history:
Numerical simulations are carried out to study the space development of hydrogeneAir
Received 4 March 2015
mixing jets discharging in to a confined environment. The full Navier-Stokes equations are
Received in revised form
solved with a high order Godunov's scheme piecewise parabolic method (PPM) with the
11 May 2015
approximate Riemann solver of Roe. The case considered in this study is based on the
Accepted 15 May 2015
experimental work of Eggers (1971). A great attention has been paid to the computation of
Available online 10 June 2015
the dynamic and mixture field, in order to analyse deeply by means of flow visualizations and statistics the characteristics phenomena of the turbulent mixing. Hydrogen species is
Keywords:
seeded in the inner jet and the air in the supersonic co-flow. The aim is to comprehend the
Hydrogen
turbulent structures effect on the mixing process of hydrogen and air. The approach is
Turbulent mixing
mainly based on Monotone Integrated Large Eddy Simulations (MILES). The study shows
Coaxial jet
the MILES ability to track mixing process in such configuration as the grid resolution is
PPM scheme
important.
MILES Structure function
Furthermore, it is found that the turbulent mixing activity is subjected to an intermittent specificity of the coherent vortices: ring structures with intermittent mushroom-type structures characterizing the ejection phenomena due to the longitudinal counter-rotating vortices. The H2 is quasi-completely dissipated as it is entrained outward the rings and transported by the irrotational stream (air). To overcome this situation, it is recommended to lock the inner jet by an outer jet. The studies with and without confinement display the same results even in 2D as the convective Mach number is subsonic Mc ¼ 0:44. Copyright © 2015, Hydrogen Energy Publications, LLC. Published by Elsevier Ltd. All rights reserved.
Introduction The increased interest in hydrogen utilization as principal or addition fuel has known a great importance from the point view of research studies as well as industrial applications and space vehicle propulsion system. The present investigation is focused on H2eair turbulent mixing process in
compressible coaxial jets. Gas jets are widely used in flow field of a several chemical engineering process. The comprehension of the turbulent mixing included hydrogen helps for the new design and performance improvement of devices such combustion chamber, injection systems. A better hydrogen and air turbulent mixing is of crucial importance to lead to higher combustion efficiency and performance improvement.
* Corresponding author. E-mail addresses: [email protected] (R. Ouzani), [email protected] (M. Si-Ameur). http://dx.doi.org/10.1016/j.ijhydene.2015.05.095 0360-3199/Copyright © 2015, Hydrogen Energy Publications, LLC. Published by Elsevier Ltd. All rights reserved.
9540
i n t e r n a t i o n a l j o u r n a l o f h y d r o g e n e n e r g y 4 0 ( 2 0 1 5 ) 9 5 3 9 e9 5 5 4
Important investigations dealing with coaxial jets have been exerted by several researchers, but in incompressible flows situation. Indeed, coaxial jets in transition have been extensively studied in the past [1e7], in order to comprehend the dynamical mechanism and mixing including such shear layers interaction, the vortex shedding frequencies and recirculation. However, limited studies on multiple species mixing in high speed coaxial jets are available in the open literature. Furthermore, hydrogen is the appropriate fuel in supersonic combustion applications for air breathing high-speed propulsion. Experimental laboratory research studies are delicate and require many resources, the open door is then numerical simulations. In this context, the present study tracks the development of the mixing mechanism in H2eair compressible jets by means of flow visualization and statistics. The aim is to better understand the turbulent mixing by using numerical simulations in three-dimensional configuration as a tool of investigation. This study can be regarded as complement to available investigations [8e10] dealing with supersonic combustion of H2eAir species. Undeniably, turbulent combustion results from the interactions between: (i) thermodynamics and heat release associated to chemical reactions, (ii) complex chemistries between several species through many reaction steps, (iii) advective transport and diffusion. The velocity field is not only transports the various species but also depends on heat released by reactions, which affects the density and then the momentum balance. The mean time step depends also on several chemical reaction rates which change with concentration and temperature fluctuations. A large number of elementary chemical reactions consume H2 and Air which modify greatly the interface turbulent diffusion process (turbulent mixing). Reactive flows involve less intricate interactions, especially when heat release is insignificant [11], as in highly diluted liquid flows or for turbulent diffusion of reactive pollutant in the atmosphere and ocean. Accordingly, it is important to carry out studies only based on turbulent mixing in parallel to combustion works. Better
comprehension of turbulence-mixing interactions lead to improve the control of the flame ignition and extinction processes in high speed flows. Thus, chemical reactions are avoided in this present study in order to concentrate on the effect of the turbulent structures on the H2eAir mixing process. A great attention has been paid to detailed evolution of the spatial and temporal evolution of the mixing fraction, while benefiting from the advantage of the spatially developing approach. Numerical description is strengthened by means of snapshots visualization of the flow and statistics based on temporal and spatial information of the compressible flows. The numerical simulations are complex because of supersonic/subsonic flows transition, mixing mechanism between hydrogen and air, sound waves, convection of species. We are interested in the far-off mixing jets (Fig. 1) without in quest of to represent the mechanisms of generation of shearing associated to the nozzle exit. The experiment carried out by Eggers [12] gives us an opportunity for comparing numerical simulations data and laboratory high Reynolds flows. In order to achieve significance numerical computations with well reproduced features especially in spatial development, the practical and adequate technique is the Large Eddy Simulation (LES) approach. The present 3-D computations at high Reynolds number are performed using the Monotone Integrated Large Eddy Simulation (MILES) technique. The nonlinear flux limiting PPM algorithm, initially conceived to solve gasodynamics shocked problems, have been used with increasing assurance to compute several high-Reynolds number flows in astrophysics [13]. The effects of subgrid scales on the resolved flow are implicitly included into the convective fluxes; the PPM intrinsically acts as a high frequency filter and work adequately near the grid scale cut-off. In this paper, an academic test (temporal mixing layer) is carry out to analyse and to justify the adopted approach. It is shown that the large scales are insensible to an explicit subgrid model which leads to worthlessness contribution and reduced performance in the turbulence spectrum. Contrary to classical schemes which produce numerical viscosity, the PPM advection scheme by approximating the flow variable evolution along each cell with a monotonous
Fig. 1 e (a) schematic configuration of the injectors (Eggers [12]), (b) Sketch of spatially developing jets and computational domain.
i n t e r n a t i o n a l j o u r n a l o f h y d r o g e n e n e r g y 4 0 ( 2 0 1 5 ) 9 5 3 9 e9 5 5 4
parabola, behaves as a subgrid scale filter for scales smaller than several grid sizes. This subgrid model adequately couples the large resolved scales with the unresolved subgrid motion through conservation laws. This method is highly localized in space and time, induces insignificant dissipation for the well resolved scales. The motivated tests of a compressible turbulent flow [14] using PPM scheme, without an explicit subgrid turbulence model, have outlined the fact that the solenoı¨dal modes (incompressible modes) build-up a Kolmogorov k5/3 spectrum. Also, it is found that PPM displays an effective dissipation that reliable approximately the fourth power of the wave number, which means that the residual dissipation handles mainly the small scales and do not affect the correct energy cascade. Karaca [15] has also shown that Implicit large eddy (ILES) simulation with 5th order WENO scheme lead to good result compared to classic LES as a grid convergence is reached. The present paper is organised as follows. In section Governing equations and the numerical method an overview of the mathematical model of the physical problem is outlined, with emphasis on the numercial method used for the below computations. In section High speed H2eAir jets, the numerical results are presented through the dynamics and scalar mixing fields, comparaison with laboratory experiments of Eggers [12] are given. Finally, the conclusion is drawn in section Conclusion.
Governing equations and the numerical method Numerical simulations are performed here, by solving the complete compressible Navier-Stokes equations. The fourth order accurate piecewise parabolic method (PPM) together with the approximate Riemann Solver of Roe solve the inviscid equations. The overall scheme is second order accurate in time. The diffusion terms are regrouped in a source term of Navier-Stokes equations. Computations of high-speed jets are performed using a code specifically developed to numerically simulate gas flows whose composition may evolve with location ! x and time t. It is based on the usual set of gas flow equations, with the vector of the conserved variables is: Uð! x ; tÞ ¼ fr; rU1 ; rU2 ; rU3 ; re; rYi gT
(1)
r is the density, ðU1 ; U2 ; U3 Þ the three velocity components, e the total energy per unit mass, Yi is the mass fraction of the species i. The total energy of the flow is defined by: e¼
1 2 U þ U22 þ U23 þ Cv T 2 1
(2) Pns
where the average heat capacity per unit mass Cv ¼ i¼1 CviYi (ns is the number of involved species). Thermo-dynamics of the gas mixture is modelled with an ideal gas equation. The pressure P and the temperature T are related by: P RT ¼ r M
9541
Each component of Uð! x ; tÞ corresponds to a conserved quantity like mass, momentum or energy. The total energy e is the sum of the mechanical and internal energies of gas mixture. The equations for conservation of mass, momentum, energy and detailed mass of species i are written in conservative form. 3 vFj ðUð! x ; tÞÞ vUð! x ; tÞ X þ ¼ SðUð! x ; tÞÞ vt vxj j¼1
(4)
8 9 rUj > > > > > > > rUj U1 þ dj1 P > > > > > > > > < rUj U2 þ dj2 P > = ! Fj ðUð x ; tÞÞ ¼ rUj U3 þ dj3 P > > > > ðre þ pÞUj > > > > > > > > > > rYi Uj > > : ; … where: 8 > > > > > > > > > > > > > > > > > > > > > > > > > > > <
0 9 r > > > j¼3 > P v Re t1j > > > > > vx > j j¼1 > > > r > > > > j¼3 > P v Re t2j > > > > > vx > j j¼1 > > > = r t v SðUÞ ¼ j¼3 3j P Re > > > > > > > > vxj > > j¼1 > > > > > > > > > > > > > Cpr r ! > > > > > VT þ V$ U :t V$ > > > > Pr:Re Re > > > > > > > > > > r > > > > VY V$ > > i > > > > Sc:Re > > : ; … The right-hand side term of equation (4) contains terms of molecular flux (momentum, heat and mass). P3 vUk vUj 2 i tij ¼ vU þ d . k¼1 vxk vx vx 3 ij j
i
€ necker symbol, Re,Pr and Sc are the Reynolds dij is the Kro number, Prandtl number and Schmidt number of species i, respectively. They are defined in section High speed H2eAir jets. m, l, Di are the dynamical viscosity coefficient, heat conduction coefficient, species diffusion coefficient. Cp is the heat capacity at constant pressure of the mixture per unit mass. The values of Cp, Cv and M evolve and have to be calculated at each ð! x ; tÞ. The Cp value changes with the local composition of the H2eAir mixture and is dogged by a mass weighted over the heat capacities of the H2 and Air. The full compressible Navier-Stokes equation (4) are split up in three one-dimensional inviscid operators Lj ðj ¼ 1; 2; 3Þ: Lj :
! vUð! x ; tÞ vFj ðUð x ; tÞÞ ¼0 þ vt vxj
(5)
(3)
where R is the universal gas constant and M is the average molar mass calculated from the molar mass Mi of each species Ps 1 i: M ¼ ni¼1 Yi =Mi .
and three-dimensional viscous-diffusive operator: j:
vUð! x ; tÞ ¼ SðUð! x ; tÞÞ vt
(6)
9542
i n t e r n a t i o n a l j o u r n a l o f h y d r o g e n e n e r g y 4 0 ( 2 0 1 5 ) 9 5 3 9 e9 5 5 4
The fourth order accurate piecewise parabolic method (PPM) [16] together with the approximate Riemann solver of Roe [17] solve the equation (5). The spatial discretization of equation (6) is done by finite volume. The integration in time is accomplished by second order accurate low storage Runge-Kutta time stepping. The averaged values Ui over the ith grid cell are updated by computed the fluxes through the Riemann Solver between two discontinuities at each cell interface over one timestep. The space higher-order precision is achieved by taking into account the spatial variation of the discretized function in order to built two states UL (left) and UR (right) as an average of the information which is crossing to the interface over one timestep. The PPM algorithm generates a parabolic interpolation function for the primitive dynamic and scalar fields over each cell. The parabola is defined by its average over the cell and the border values which are computed with a fourth-order differencing formula. They are corrected where they generate new extrememums in order to realize monotonicity. The values of the primitive variables from both sides to each interface during one timestep, are acquired by averaging the piecewise parabolic profile over the appropriate domain of dependence. For a wave travelling with the speed c, the average is defined by: L ðvðxÞ; cÞ rightward waves : fiþ1=2 8 xiþ1=2 Z > > > 1 > < vðxÞdx if c > 0 ¼ cDt x cDt > iþ1=2 > > > : viþ1=2 if c 0
8 > > > > < 1 R cDt leftward waves : fiþ1=2 ðvðxÞ; cÞ ¼ > > > > :
(7)
Z
vðxÞdx viþ1=2
R
RLk ðxÞeLk ; VðxÞ Viþ1=2 ¼
k¼1
if c < 0
k¼3 X
RRk ðxÞeRk
(9)
(10)
k¼1
where Rk ðxÞ is the amplitude and ek ðxÞ the eigenvector of the linear wave associated with the kth characteristic. The Roe fluxes are determined after a Roe average by:
Fi±1=2 ¼
F ULi±1=2 þ F URi±1=2 2
(14)
The spatial discretization of these terms is done in finitevolume form. The fluxes and the stress tensor ! P3 vUk vUj 2 i t tij ¼ vU þ d are computed on the cell ink¼1 vxk vx vx 3 ij j
i
terfaces and are differenced conservatively. The combined solution of the three hyperbolic Lj ðj ¼ 1; 2; 3Þ and parabolic j problems are computed as follows: 1 1 1 1 Unþ1 ¼ j2 L1 L2 L3 j2 Un , Unþ1 ¼ j2 L3 L2 L1 j2 Un alternatively. Variables at time-step n þ 1 are computed from variables at time n. The timestep control is based on known criteria especially the CFL number which permits to obtain an individual time is step for convection operator Lj ðj ¼ 1; 2; 3Þ. Dtx < min clDx þcr the CFL number (cl , cr : left and right sound velocity). An optimal global value of the time step Dt is selected at each time inside the limit of stability to guarantee efficient computations. Dtmax ¼ minðDtx ; Dty ; Dtz ; Dtj Þ. The key param-
2Dmax
; Dmax ¼ max Di ; mr;
1
1 þ 1 þ 1 Dx2 Dy2 Dz2
l Cpr
.
if c 0
R
k¼3 X
(13)
xiþ1=2
L R Viþ1=2 ¼ fiþ1=2 ðVðxÞ; ui þ ci Þ ; Vi1=2 ¼ fi1=2 ðVðxÞ; ui ci Þ
L
(12)
/ 1 vr! u ¼ V $ $t vt Re vre ! Cpr ! ! 1 ! ¼ V$ VT þ V$ U :t vt Pr:Re Re
heat). Dtj <
xiþ1=2 cDt
The primitive vector of the dynamical and scalar fields on the left and on the right of the cells border, are decomposed into upstream and downstream travelling pressure wave characteristics. The left and right states of the primitive vector are linearised about:
VðxÞ Viþ1=2 ¼
vrYi ! r ! V Yi ¼ V$ Sc:Re vt
eters taken into account are the computational grid size, convection, and diffusion process (momentum, mass and !
(8)
L
The equations for species, momentum diffusion and energy conservation (sources terms of the equation (6)) are respectively:
k¼3 X a* R* e* k k k
(11)
k¼1
The scheme is of fourth order in time and in space where the solution is smooth and second order in space and third order in time near discontinuities.
Initial and boundary conditions In order to close the equation (4), boundary and initial conditions must be specified. In this free flow, boundary conditions behave as a model of the flow outside the computational domain. These conditions are given under a piecewise formulation. The boundaries are tailored in order to let mass fluxes; waves and flow structures freely escape the domain when mass fluxes only are permitted to enter the computation domain. Since the PPM-scheme integrates a differential equation along each ðx; tÞ-characteristics with an upwind scheme, the boundary conditions are written in a natural way. The inflow is sensitive to the upstream boundary condition with the mean shear originating the instabilities which develop downstream. A white noise is superimposed to a " !#! basic profile: uðrÞ ¼ Ua þ
Uj Ua 2
1 tanh
Rj 4q
Rj r
Rrj
in order
to favour unstable modes and model the residual turbulence of the upstream flow. Rj is the injector jet Raduis. A CroccoeBusemann's relation is also implemented to model the heating due to molecular dissipation as an effect of
i n t e r n a t i o n a l j o u r n a l o f h y d r o g e n e n e r g y 4 0 ( 2 0 1 5 ) 9 5 3 9 e9 5 5 4
9543
Fig. 2 e Test on outflow boundary condition.
ra compressibility rðrÞ ¼1þ
1rra j
2
" tanh
Rj 4q
!# Rj r
Rrj
! 1 .
The initial conditions have little physical significance although they are needed to start a computation. They can be somewhat arbitrary provided that the analysis of results starts at a time which is large enough to let phenomena develop with little dependence on initial conditions. The flow field is initialized over the whole domain with its upstream values. In cross stream and spanwise directions, free slip walls (or outflow without wall) and periodic boundary conditions are imposed respectively. At the downstream boundary, all the characteristics ðx; y; z; tÞ are transported outwards, to avoid the dependence between the flow and the outflow boundary. For the exit characteristic a Newmann condition is imposed: Unx þ1 ¼ Unx . For the incoming characteristic a Dirichlet condition is fixed to zero to avoid the perturbation of the upstream flow. The H2 central jet is subsonic, the outflow boundary should be efficient and consistently to avoid parasite reflections which can move back to upstream and then perturb the all the flow
Fig. 3 e Energy spectrum: computation with and without explicit subgrid model.
jets. To be cautiously and adequately controlled, a numerical test is performed as [18] which is based on innermost surface. The primitive dynamical variables are initialized by: 8 > > > > > > <
U1 ð! xÞ ¼ 0 ! U ðxÞ ¼ 0 2
U3 ð! xÞ ¼ 0 > > >
> > > : P ¼ P0 þ ε exp ln 2 r2 b2 As it is shown on Fig. 2, low reflection is generated; the sphere wave is submitted to negligible distortion, which gives confidence to our outflow boundary.
Academic test: temporal mixing layer A test configuration (temporal mixing layer) is analysed to justify and validate the approach (MILES) of the present heart article. Indeed, an explicit subgrid model based on filtered structure function [19] coupled to the flux limiters of the PPM
Fig. 4 e Isosurface of vorticity modulus: Grid resolution 64 64 64.
9544
i n t e r n a t i o n a l j o u r n a l o f h y d r o g e n e n e r g y 4 0 ( 2 0 1 5 ) 9 5 3 9 e9 5 5 4
Fig. 5 e H2 field in temporal mixing layer at different resolution: a) 48 48 48 b) 64 64 64 c) 96 96 96 d).128 128 128.
Fig. 6 e Energy spectrum at different grid resolutions.
Fig. 7 e Radial velocity profile: grid convergence test.
scheme is considered. The numerical implementation is accomplished on the basis of conventional LES for compressible flows approach. The Favre averaging is explicitly
~ ¼ rf. The filtered conserved vector is: taken into account by: f r ~ ~ ~ ~ i Þ. ~ U ¼ ðr; rU1 ; rU2 ; rU3 ; re; rY The filtered Navier-Stokes equations are written under this new conservative form.
i n t e r n a t i o n a l j o u r n a l o f h y d r o g e n e n e r g y 4 0 ( 2 0 1 5 ) 9 5 3 9 e9 5 5 4
Table 1 e CPU for different grid resolution. Grid resolution
CPU/Step/Node/Iteration
875 300 300 875 150 150 436 75 75
2.8271143E-06 (s) 2.2263755E-06 (s) 1.4447446E-06 (s)
Table 2 e Flow parameters, data from Eggers [12]. Parameter
Hydrogen jet
Air jet
1074 260=300 0:886 46:68 102 1:32 105 0:44 0.9 0.7
394 222=300 1:32
Uðm=sÞ Tstat =Ttot ðKÞ Mach rj U2j =ra U2a r U Djet Rejet ¼ jet mjet jet Ujet Uair Mc ¼ ajet þaair m ScH2 ¼ rH DH Cp m 2 2 Pr ¼ l H H 2 2
3 vFj ðUÞ vU X þ ¼ SðUÞ vt vxj j¼1
(15)
The filtered convective flux is: 8 9 > > > > > > ej > > rU > > > > > > > ~ > rUj U1 þ rRTdj1 > > > > > > < = ~ rUj U2 þ rRTdj2 Fj ðUÞ ¼ ~ j3 > > > rUj U3 þ rRTd > > > > > > > > > rðe þ RTÞU j > > > > > > > > rYi Uj > > : ; i h ~iU ~iU ~ j þ rUi Uj rU ~j rUi Uj ¼ rU |fflfflfflfflfflfflfflfflfflfflfflfflffl{zfflfflfflfflfflfflfflfflfflfflfflfflffl} Tij
h i ~ U ~ j þ rðe þ RTÞUj r e~ þ RT ~ U ~j rðe þ RTÞUj ¼ r e~ þ RT |fflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflffl{zfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflffl} Qij
i h ~ j þ rYi Uj rYi U ~j rYi Uj ¼ rYi U |fflfflfflfflfflfflfflfflfflfflfflfflffl{zfflfflfflfflfflfflfflfflfflfflfflfflffl} Gij
Tij , Qij , Gij are respectively the subgrid flux terms (momentum, heat, mass). They are modelled on the basis of eddy
9545
viscosity nt and eddy diffusivity Dt formulation. This model is built on proportionality between the subgrid stress and a local strain calculated from the large scale velocity field. In this case, eddy viscosity and diffusivity concepts are supported by physical arguments about damping of large scales by subgrid energy and k5/3 equilibrium spectra, nevertheless it misrepresent the so-called backscatter from the short waves to large ones and the interactions between short waves. The boundary conditions and the physical parameters are as follows: periodic in streamwise and spanwise. In transversal direction, an outflow condition is Imposed. The Reynolds number is about 10þ6, convective Mach number Mc ¼ 0:3. Four cases are analysed hereafter regarding the grid resolutions. Fig. 3 shows the turbulence spectrum with and without an additional physical subgrid scale model. The spectrum decay of the turbulent kinetic energy is satisfactory reproduced in the two cases. However, added an explicit subgrid model lead to worthlessness contribution and reduced performance. Furthermore, removing the limiters will affect the stability of the scheme as non-physical oscillations (Gibbs instabilities in the presence of strong gradients) grow, especially at low resolutions. The structure of the flow field is illustrated in the Fig. 4, the formation of coherent structures even at high Reynolds number. The flow topology resembles to that of Comte et al. [20]. In addition, the presence of vortex structures at small scales is obviously detected. Fig. 5 shows that PPM detect adequately, as the grid resolution is sufficiently improved, all the main compressible turbulent shear layer features including turbulent mixing. Increasing the resolution let to show more details and highly topology flow quality. The analysis of the spectra shows clearly and quantitatively its good tendency with increasing resolution (Fig. 6). In inertial range scales, a well characterized 5/3 power enveloppe over about one decade is visible which fill out the spectrum progressively by smaller scales as the grid resolutions is increased. These observations indicate that the large scales are well represented with a minimum dependence against the short one as the computational grid increases. This result corroborate the fact that this sophistical powerful scheme is eligible for the computation of high Reynolds compressible turbulent flows shocked or not, as nonlinear filter built-in the algorithm work adequately especially near the grid scale cut-off.
Fig. 8 e Cross section of H2 mass fraction on the jet axis.
9546
i n t e r n a t i o n a l j o u r n a l o f h y d r o g e n e n e r g y 4 0 ( 2 0 1 5 ) 9 5 3 9 e9 5 5 4
Fig. 9 e Cross section of vorticity field on the jet axis.
High speed H2eAir jets The numerical simulations were performed on 8 processor DELL Station with a total run time of 336 h. Several grid resolutions are tested and the appropriate one consists of 875 300 300 nodes, respectively in the streamwise, transverse and lateral directions. This choice is based on compromise between the available computation materials, the accuracy of steep dynamics and scalar gradients, robustness and significance of the results. The CPU are (see Table 1): A test convergence is carried out showing that the computations converge satisfactory to the experimental measurements as the resolution increases. Fig. 7 shows the progressive decay of the velocity as far as the grid resolution improves reducing thus the excess contribution of the unresolved scales. The target computation grid is then used in the analysis henceforth with confidence. The experiment carried out by Eggers [12] gives us an opportunity for comparing numerical simulations results and laboratory real flow scales. Computation parameters were designed to fit the real jets configuration. The jet at the injector exit is made of a Hydrogen at Mach number Mj ¼ 0.886, the adjacent jet is made of air at Mach number Ma ¼ 1.32. All gases are assumed to be ideal. Velocity, density and pressure are made non-dimensional by using the values
at the injector exit. Lengths are scaled on the diameter Dj of the hydrogen injector exit. Table 2 summarizes the flow parameters (the data from Eggers used in this study).
Results To qualitatively point up mixing with a set of pictures once the flow becomes statistically quasi-stationary and the effect of the initial condition has decayed, the instantaneous planar flow visualization is shown in Fig. 8. The natural quasi-linear region until x/Dj ¼ 5, where the instabilities grow freely, as the inlet perturbation is based on white noise, is obviously visible. The introduction of white noise perturbation avoids difficulties in specifying random perturbation with a given spectrum; this region is not representative for confrontation to experimental investigations as the causes of disturbances is linked to several sources like vibrations, pressure waves… However, it has been shown [3,4,21] that the consideration of hyperbolic profile let to develop easily the main instabilities and the flow features. For x/Dj > 7, the interaction of the two distinct mixing layers in the transition region, followed by a merging process into an inflate enclosure with a spate big spots of hydrogen. This enclosure deflates forthwith in the next nearer spatial positions which an indicative of constriction and relaxation
Fig. 10 e Visualisation of the potential core and rings vortex.
i n t e r n a t i o n a l j o u r n a l o f h y d r o g e n e n e r g y 4 0 ( 2 0 1 5 ) 9 5 3 9 e9 5 5 4
9547
Fig. 11 e Isosurfaces of Q criterion coloured by: (a) H2 masse fraction (b) streamwise velocity.
process. Downstream, one can see mixing over the entire thickness of the overall jet development. The hydrogen is sometimes dispersed top down, following the dynamics of the pulse enclosure. Spatially, the intermittent extent of the mixing is evident. The induction of the enclosure is also visible on the cut x-y vorticity field (Fig. 9) for 10 < x/Dj < 15. As for incompressible jet, the shear layer mode is evident; we can observe the development of structures on the two shear layers of the jet because of Kelvin-Helmoltz instabilities.
Beyond, eddies remain well organized and keep a high level of turbulence energy. The contrast between the vorticity field and mixing field is visible; the animations indicate different mechanisms. The potential cone is well visible in Fig. 10, the longitudinal and radial velocities remain constant near the centerline. The first structures (rings) appear at x/Dj ¼ 5 inducing a radial momentum diffusion otherwise velocity gradient. The vortex merging process linked to incompressible jet as observed by several investigations is not visible in our compressible
Fig. 12 e (a) Radial profile of the mean longitudinal velocity in the fully turbulent region (b) Radial profile of the rms of the longitudinal velocity component in the fully turbulent region.
9548
i n t e r n a t i o n a l j o u r n a l o f h y d r o g e n e n e r g y 4 0 ( 2 0 1 5 ) 9 5 3 9 e9 5 5 4
Fig. 13 e Cut through the (z,y) plane of the H2 mass fraction and vorticity fields.
situation. Furthermore, the varicose mode seems higher than the sinuous one as no helical structures are observed. Indeed, Fig. 11 shows an instantaneous view of threedimensional positive Q (is the second invariant of the velocity-gradient tensor, it is well known as a being a good indicator of turbulent structures) isosurfaces coloured by: H2 masse fraction (Fig. 11a) and streamwise velocity (Fig. 11b). The rings primary vortices appear in the early transition zone between the jet and the co-flow. Pairs of alternate streamwise vortices appear between two consecutive rings, in agreement with the classical phenomena of three-dimensionalization in
free-shear layers. The streamwise vortices play a dominant role in the transition processes towards a fully-turbulent state, and their location is very dependent of the inlet conditions [5]. The apparition of streamwise vortices is finally followed by an abrupt increase in the level of small-scale turbulence. Then the flow quickly reaches a state of fullydeveloped turbulence. A quantitative analysis is carried out, the mean and rms profiles are obtained through a temporal average of the dynamic and scalar fields. The Favre averaging is explicitly taken into account. The time averaging period is set to two times the
i n t e r n a t i o n a l j o u r n a l o f h y d r o g e n e n e r g y 4 0 ( 2 0 1 5 ) 9 5 3 9 e9 5 5 4
Fig. 14 e Downstream evolution: (a) azimuthal turbulent intensity (b) radial turbulent energy (continuous line) and azimuthal turbulent energy (dashed line).
Fig. 15 e H2 mass fraction isosurface (a) 0.07 (b) 0.2 (c) 0.5.
9549
9550
i n t e r n a t i o n a l j o u r n a l o f h y d r o g e n e n e r g y 4 0 ( 2 0 1 5 ) 9 5 3 9 e9 5 5 4
Fig. 16 e (a) Mean axial velocity profiles at several downstream sections. (b) rms axial velocity profiles at several downstream sections.(c) Mean radial velocity profile: numerical and experimental comparison.
flow time, where one flow-through time is defined as the time of the flow to traverse the computational domain. In order to compute the turbulence intensities, the relation: sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi r42 þ is used. The temporal average: 40 ¼ 42 24r4 r r Z fðx; y; zÞ ¼ T1
t0 þT t0
fðx; y; z; tÞdt. Where: fðx; y; zÞ ¼ Ui¼1;2;3 or YH2
T ¼ 4Lx =Uj is the time averaging period, t0 ¼ 4Lx =Uj is the time required to erase the effect of the initial condition. The mean and turbulent velocity profiles indicate a self similarity sate (Fig. 12). Moreover, there is a strong momentum transfer from jet to coflow. This transfer implies a spreading of the jet into co-flow. One can see high level of longitudinal turbulent intensity through the width of the jet and downstream position which is a consequence of the longitudinal vortices. However, the azimuthal instabilities are strong as the Reynolds number is high in our numerical simulations and substantially different to incompressible case. Indeed, Fig. 13 shows clearly the growing process of the azimuthal instabilities as several pairs of mushroom structures, which are responsible of diffusion, and dispersion of H2 on the jet periphery. These mushroom structures are located in the inner and outer rings. The inner ones shattered and eject H2 inside
the drained jet surface. However, the outer ones entrained air leading to the growth of mixture envelope. The inner can be due to the supersonic co-flow stream pushing the braids region within. Furthermore, Fig. 14 shows clearly an abrupt increase of azimuthal fluctuations in transition region because of a high mushrooms eddies activity. In the fully turbulent state, a decrease of these fluctuations is attributed to locking process of rings in enclosure followed by an alteration of the azimuthal rotation and the radial expansion of the rings. The production of intense ejections, in which mixture is expelled outwards from the rings are caused by the streamwise vortices, which also pump irrotational fluid (air), by entrainment. The estimated Strouhal number value is about 0.04 which obtained after a period of acquisition of the pressure field followed by a Fourier transform. The emergence of the streamwise vortices is followed by the creation of small turbulent structures where indicate the beginning of the transition to turbulence. The interaction and induction of inner and outer streamwise vortices lead to transition to turbulence which pump high level of air inside shear layer, The H2 is quasi completely dissipated as it is entrained outward the rings and transported by the irrotational stream (air). Indeed, a low values of H2 species as well as vorticity are clearly observed for x/Dj > 20. To overcome this situation, it is then recommended to lock the inner jet by an
i n t e r n a t i o n a l j o u r n a l o f h y d r o g e n e n e r g y 4 0 ( 2 0 1 5 ) 9 5 3 9 e9 5 5 4
9551
Fig. 17 e (a) Mean radial H2 mass fraction profiles at several downstream sections. (b) rms radial H2 mass fraction profiles at several downstream sections.(c) rms axial H2 mass fraction profile (d) Mean radial H2 mass fraction profile: numerical and experimental comparison.
outer jet (replacing the co-flow jet) to avoid the dissipation of H2 in turbulent zone. This recommendation is inspired from the work of Balarac and Si-Ameur [5] who found that the inner jet remain confined by the outer one in a configuration of incompressible coaxial jets. Fig. 15 shows iso-mixture fractions which are valuable for fractal dimensions and mixing efficiency. Depiction of the isosurface of H2 is important for modelling turbulent mixing. Thus, the turbulent mixing starts with an engulfment of two species through the shear layers implied by the KelvinHelmholtz vortices. The isosurface is highly convoluted after the amplification instabilities, with evidence of a mixing envelope in these jets for low value of H2 Fig. 15a, unlike the high value
where spatial intermittence is palpable Fig. 15b and c. This result show clearly that H2 is highly diffused and dispersed far from the central jet until his complete quasi dissipation. This disadvantage is due to significant momentum and mass transfer by the longitudinal eddies, towards the coflow jet. This observation constitutes an important stage in the evolution of H2-air mixing. This finding is useful for combustion process especially in the phase of flame ignition as molecular mixing is primordial to ensure a normal ignition process. The Fig. 16a shows the radial variation of the mean streamwise velocity profile at several downstream positions. At the beginning of transition, the jet maintains a profile analogous to the inlet boundary profile. Continuous decreases
Fig. 18 e Transverse jet section coloured by probability density of H2.
9552
i n t e r n a t i o n a l j o u r n a l o f h y d r o g e n e n e r g y 4 0 ( 2 0 1 5 ) 9 5 3 9 e9 5 5 4
Fig. 19 e Mean radial (velocity and H2 mass fraction) profiles: numerical and experimental comparison a) x/D ¼ 9.58, b) x/ D ¼ 15.44, c) x/D ¼ 25.2.
of the maximum mean velocity along the jet axis, followed by a spreading of jet into the co-flow as a consequence of momentum transfer, are evident and similar to classical incompressible jets. However, the turbulent activity is located at the interface at the beginning of transition (Fig. 16b), then it spread on the total width. Two peaks are observed which is consistent with the visualizations of two mushrooms structures (sublayers inner and outer). A continuous increase is located in the enclosure followed by abrupt decrease of streamwise turbulent intensity. Our observations are
corroborating by the good agreement with the measurements of Eggers, as it is evident on Fig. 16c. The Fig. 17a illustrates the average radial profile of the H2 at with various sections of the jet. Two stages can be distinguished; upstream the H2 profiles evolves gradually; the increases in mixture zone is due to radial and longitudinal evolutions noticeably the pulsatile movement of the ring structures permitting the penetration of air in H2 jet; the rms of the mixing fraction is very low in this region (Fig. 17b). Two peaks are also evident which are the signature of the two
i n t e r n a t i o n a l j o u r n a l o f h y d r o g e n e n e r g y 4 0 ( 2 0 1 5 ) 9 5 3 9 e9 5 5 4
9553
Fig. 20 e 2D computations: a) H2 field, b) pressure field.
sublayers mushrooms structures. Fig. 17c highlights a continuous increase followed by abrupt decrease of the H2 rms, indicating blockage of radial and azimuthal rings development (pulsating movement is altered in the enclosure). It follows longitudinal ejections which only transfer the H2 in the co-flow and inducing its quasi dissipation. The comparison between the experimental measurements of Eggers (1971) and the present numerical computations (Figs. 16c and 17d) shows a good agreement qualitatively and quantitatively. The turbulent mixing activity is confined in narrow thickness jr/Djj<1 in transition region, however it extends downstream over twice, as it is visible on Fig. 16a. The figure shows that a significant quantity of air invaded the jet center. There is thus a strong transfer of H2 to ambient jet towards the central jet implied by the transfer of momentum that we already evoked. The study of the fluctuations of H2 (Fig. 17b) reveals good tendency compared to experimental data. Moreover, even the small values of H2 mainly within the periphery of the jet downstream, a considerable scalar fluctuations, on the axis of the jet, are detected. That is explained by a strong intermittency of the turbulent mixing activity, as it is shown on Fig. 18, which presents a transverse jet section coloured by probability density of H2. The probability to find a high mass fraction of H2 is small on all the jet width, which shows clearly his quasi-dissipation mentioned above. Furthermore, the strong intermittency in the area of the central jet, undoubtedly due to the injections implied by the longitudinal vortices. A further confrontation of the present computations against experimental at different sections have revealed good agreement. Fig. 19 shows clearly this the good performance of the PPM to track dynamics and mixing in high-speed shear flows comparatively to WENO 5 scheme [15].
Confinement effect All the results of the subsection Results are obtained in confinement environment. However, another study has been undertaken by removing the walls and adopting outflow boundary conditions in transverse direction. The exact same results are obtained for unconfined situation which means that the walls do not play any role in the dynamic development of the jets. This can be explained by the value of the
convective Mach number about 0.44. The wall effects can be important when Mc > 1 as it is demonstrated by Si-Ameur & Chollet [22]. An examination of the 2D computations has revealed the same finding. Indeed the Kelvin-Helmholtz instabilities, the interaction of the large turbulent structures with the mixing process are unchanged in the two situation. The interaction of the mixing zone with eventually Mach wave is insignificant in the case of the confinement (Fig. 20).
Conclusion A three-dimensional spatially developing jets, of two supersonic-subsonic streams of Hydrogen and air gas has been performed using the Monotone Integrated Large eddy simulation approach. The computations explicitly include the near field and the far field of the jets development with an emphasis on mixing. A locking process of rings in enclosure followed by an alteration of their azimuthal rotation and radial expansion are observed, through flow visualizations and statistics. Furthermore, H2 is highly diffused and dispersed far from the central jet until his complete quasi-dissipation. After the transition, the flow structure is directed preferentially in the streamwise direction due to the longitudinal turbulent scales, which are ejected from the enclosure. These structures particularly entrain air stream continuously inside jet leading to quasi-complete alteration of H2. The examination of the transverse jet coloured by PDFs has revealed the existence of strong intermittency in the area of the central jet, undoubtedly due to the injections implied by the longitudinal vortices. A continuous increase followed by abrupt decrease of the H2 rms, indicating blockage of radial and azimuthal rings development (pulsating movement is altered in the enclosure). It follows longitudinal ejections which only transfer the H2 in the co-flow and inducing its near dissipation. The confined environment do not affect the hydrogen-air mixing mechanism as the convective Mach number is subsonic. A comparison with laboratory experiment has shown a good agreement as well as for the mean and turbulent
9554
i n t e r n a t i o n a l j o u r n a l o f h y d r o g e n e n e r g y 4 0 ( 2 0 1 5 ) 9 5 3 9 e9 5 5 4
dynamic and scalar fields. On se basis of this study, the mixing aireH2 jets of this configuration is not suitable for combustion application, as the fuel can leave completely the central jet and ejected outflow by the supersonic co-flow. It is recommended to lock the inner jet by an outer jet to avoid this probable situation.
references
[1] Villermaux E, Rehab H. Mixing in coaxial jets. J Fluid Mech 2000;425:161e85. [2] Ritchie BD, Mujumdar DR, Seitzman JM. Mixing in coaxial jets using synthetic jet actuators. 2000. AIAA Paper. tais O, Lesieur M. Direct numerical [3] Balarac G, Si-Ameur M, Me simulations of high velocity ratio coaxial jets: mixing properties and upstream conditions influence. J Turbul 2007;8(22). tais O, Lesieur M. Large Eddy [4] Balarac G, Si-Ameur M, Me simulation of coaxial jets: coherent structures and mixing properties. Direct Large Eddy Simul VI 2006;10:277e84. [5] Balarac G, Si-Ameur M. Mixing and coherent vortices in turbulent coaxial jets. C R Mec 2005;333:622e7. [6] Reynier P, Minh HH. Numerical prediction of unsteady compressible turbulent coaxial jets. Comput Fluids 1998;27(2):239e54. [7] Klioutchnikov I, Olivier H, Odenthal J. Numerical investigation of coaxial jets entering into a hot environment. Comput Fluids 2013;86:490e9. [8] Jin Tai, Luo Kun, Lu Shuqiang, Fan Jianren. Direct numerical simulation of a supersonic lifted hydrogen jet flame: a priori study on combustion models. Acta Astronaut 2015;109:52e64. [9] Luo Kun, Jin Tai, Lu Shuqiang, Fan Jianren. DNS analysis of a three-dimensional supersonic turbulent lifted jet flame. Fuel 2013;108:691e8.
[10] Moule Y, Sabelnikov V, Mura A. Highly resolved numerical simulation of combustion in supersonic hydrogeneair coflowing jets. Combust Flame 2014;161:2647e68. [11] Si Ameur M. Isothermal reactive mixing layer: numerical study. Comput Therm Sci 2011;3:227e35. [12] Eggers JM. Turbulent mixing of coaxial compressible hydrogeneair jets. Tech. Rep. NASA-TN D-6487. 1971. [13] Porter DH, Pouquet. Kolmogorov-like spectra in decaying three-dimensional supersonic flows. Phys Fluids 1994;6:2133e42. [14] Porter DH, Pouquet A, Woodward PR. A numerical study of supersonic turbulence. Theor Comput Fluid Dyn 1992;4:13e49. [15] Karaca M, Lardjane N, Fedioun I. Implicit large eddy simulation of high speed non reacting and reacting air/H2 jets with a 5th order WENO scheme. Comput Fluids 2012;62:25e44. [16] Colella P, Woodward PR. The piecewise parabolic method (PPM) for gas-dynamical simulations. J Comput Phys 1984;54:174e201. [17] Roe PL. Approximate Riemann solvers, parameter vectors and difference schemes. J Comput Phys 1981;43:357e72. [18] Lodato G, Domingo P, Vervisch L. Three-dimensional boundary conditions for direct and large-eddy simulation of compressible viscous flows. J Comput Phys 2008;227:5105e43. [19] Ducros F, Comte P, Lesieur M. Large-eddy simulation of transition to turbulence in a boundary layer developing spatially over a flat plate. J Fluid Mech 1996;326:1e36. [20] Comte P, Lesieur M, Lamballais E. Large and small-scale stirring of vorticity and a passive scalar in a 3D temporal mixing layer. Phys Fluids A 1992;4:2761e78. rique des effects de tempe rature dans [21] Daviller G. Etude nume les jets simples et coaxiaux. PhD thesis, ISAE-ENSMA Poitiers. 2010. [22] Si Ameur M, Chollet JP. Transition to turbulence of supersonic mixing layer. C R Acad Sci Paris Iib 1998;326:315e22.
Available online at www.sciencedirect.com
ScienceDirect journal homepage: www.elsevier.com/locate/he
Numerical study of hydrogeneair mixing in turbulent compressible coaxial jets R. Ouzani*, M. Si-Ameur LESEI Laboratory, Mechanical Engineering Department, University of Batna, Algeria
article info
abstract
Article history:
Numerical simulations are carried out to study the space development of hydrogeneAir
Received 4 March 2015
mixing jets discharging in to a confined environment. The full Navier-Stokes equations are
Received in revised form
solved with a high order Godunov's scheme piecewise parabolic method (PPM) with the
11 May 2015
approximate Riemann solver of Roe. The case considered in this study is based on the
Accepted 15 May 2015
experimental work of Eggers (1971). A great attention has been paid to the computation of
Available online 10 June 2015
the dynamic and mixture field, in order to analyse deeply by means of flow visualizations and statistics the characteristics phenomena of the turbulent mixing. Hydrogen species is
Keywords:
seeded in the inner jet and the air in the supersonic co-flow. The aim is to comprehend the
Hydrogen
turbulent structures effect on the mixing process of hydrogen and air. The approach is
Turbulent mixing
mainly based on Monotone Integrated Large Eddy Simulations (MILES). The study shows
Coaxial jet
the MILES ability to track mixing process in such configuration as the grid resolution is
PPM scheme
important.
MILES Structure function
Furthermore, it is found that the turbulent mixing activity is subjected to an intermittent specificity of the coherent vortices: ring structures with intermittent mushroom-type structures characterizing the ejection phenomena due to the longitudinal counter-rotating vortices. The H2 is quasi-completely dissipated as it is entrained outward the rings and transported by the irrotational stream (air). To overcome this situation, it is recommended to lock the inner jet by an outer jet. The studies with and without confinement display the same results even in 2D as the convective Mach number is subsonic Mc ¼ 0:44. Copyright © 2015, Hydrogen Energy Publications, LLC. Published by Elsevier Ltd. All rights reserved.
Introduction The increased interest in hydrogen utilization as principal or addition fuel has known a great importance from the point view of research studies as well as industrial applications and space vehicle propulsion system. The present investigation is focused on H2eair turbulent mixing process in
compressible coaxial jets. Gas jets are widely used in flow field of a several chemical engineering process. The comprehension of the turbulent mixing included hydrogen helps for the new design and performance improvement of devices such combustion chamber, injection systems. A better hydrogen and air turbulent mixing is of crucial importance to lead to higher combustion efficiency and performance improvement.
* Corresponding author. E-mail addresses: [email protected] (R. Ouzani), [email protected] (M. Si-Ameur). http://dx.doi.org/10.1016/j.ijhydene.2015.05.095 0360-3199/Copyright © 2015, Hydrogen Energy Publications, LLC. Published by Elsevier Ltd. All rights reserved.
9540
i n t e r n a t i o n a l j o u r n a l o f h y d r o g e n e n e r g y 4 0 ( 2 0 1 5 ) 9 5 3 9 e9 5 5 4
Important investigations dealing with coaxial jets have been exerted by several researchers, but in incompressible flows situation. Indeed, coaxial jets in transition have been extensively studied in the past [1e7], in order to comprehend the dynamical mechanism and mixing including such shear layers interaction, the vortex shedding frequencies and recirculation. However, limited studies on multiple species mixing in high speed coaxial jets are available in the open literature. Furthermore, hydrogen is the appropriate fuel in supersonic combustion applications for air breathing high-speed propulsion. Experimental laboratory research studies are delicate and require many resources, the open door is then numerical simulations. In this context, the present study tracks the development of the mixing mechanism in H2eair compressible jets by means of flow visualization and statistics. The aim is to better understand the turbulent mixing by using numerical simulations in three-dimensional configuration as a tool of investigation. This study can be regarded as complement to available investigations [8e10] dealing with supersonic combustion of H2eAir species. Undeniably, turbulent combustion results from the interactions between: (i) thermodynamics and heat release associated to chemical reactions, (ii) complex chemistries between several species through many reaction steps, (iii) advective transport and diffusion. The velocity field is not only transports the various species but also depends on heat released by reactions, which affects the density and then the momentum balance. The mean time step depends also on several chemical reaction rates which change with concentration and temperature fluctuations. A large number of elementary chemical reactions consume H2 and Air which modify greatly the interface turbulent diffusion process (turbulent mixing). Reactive flows involve less intricate interactions, especially when heat release is insignificant [11], as in highly diluted liquid flows or for turbulent diffusion of reactive pollutant in the atmosphere and ocean. Accordingly, it is important to carry out studies only based on turbulent mixing in parallel to combustion works. Better
comprehension of turbulence-mixing interactions lead to improve the control of the flame ignition and extinction processes in high speed flows. Thus, chemical reactions are avoided in this present study in order to concentrate on the effect of the turbulent structures on the H2eAir mixing process. A great attention has been paid to detailed evolution of the spatial and temporal evolution of the mixing fraction, while benefiting from the advantage of the spatially developing approach. Numerical description is strengthened by means of snapshots visualization of the flow and statistics based on temporal and spatial information of the compressible flows. The numerical simulations are complex because of supersonic/subsonic flows transition, mixing mechanism between hydrogen and air, sound waves, convection of species. We are interested in the far-off mixing jets (Fig. 1) without in quest of to represent the mechanisms of generation of shearing associated to the nozzle exit. The experiment carried out by Eggers [12] gives us an opportunity for comparing numerical simulations data and laboratory high Reynolds flows. In order to achieve significance numerical computations with well reproduced features especially in spatial development, the practical and adequate technique is the Large Eddy Simulation (LES) approach. The present 3-D computations at high Reynolds number are performed using the Monotone Integrated Large Eddy Simulation (MILES) technique. The nonlinear flux limiting PPM algorithm, initially conceived to solve gasodynamics shocked problems, have been used with increasing assurance to compute several high-Reynolds number flows in astrophysics [13]. The effects of subgrid scales on the resolved flow are implicitly included into the convective fluxes; the PPM intrinsically acts as a high frequency filter and work adequately near the grid scale cut-off. In this paper, an academic test (temporal mixing layer) is carry out to analyse and to justify the adopted approach. It is shown that the large scales are insensible to an explicit subgrid model which leads to worthlessness contribution and reduced performance in the turbulence spectrum. Contrary to classical schemes which produce numerical viscosity, the PPM advection scheme by approximating the flow variable evolution along each cell with a monotonous
Fig. 1 e (a) schematic configuration of the injectors (Eggers [12]), (b) Sketch of spatially developing jets and computational domain.
i n t e r n a t i o n a l j o u r n a l o f h y d r o g e n e n e r g y 4 0 ( 2 0 1 5 ) 9 5 3 9 e9 5 5 4
parabola, behaves as a subgrid scale filter for scales smaller than several grid sizes. This subgrid model adequately couples the large resolved scales with the unresolved subgrid motion through conservation laws. This method is highly localized in space and time, induces insignificant dissipation for the well resolved scales. The motivated tests of a compressible turbulent flow [14] using PPM scheme, without an explicit subgrid turbulence model, have outlined the fact that the solenoı¨dal modes (incompressible modes) build-up a Kolmogorov k5/3 spectrum. Also, it is found that PPM displays an effective dissipation that reliable approximately the fourth power of the wave number, which means that the residual dissipation handles mainly the small scales and do not affect the correct energy cascade. Karaca [15] has also shown that Implicit large eddy (ILES) simulation with 5th order WENO scheme lead to good result compared to classic LES as a grid convergence is reached. The present paper is organised as follows. In section Governing equations and the numerical method an overview of the mathematical model of the physical problem is outlined, with emphasis on the numercial method used for the below computations. In section High speed H2eAir jets, the numerical results are presented through the dynamics and scalar mixing fields, comparaison with laboratory experiments of Eggers [12] are given. Finally, the conclusion is drawn in section Conclusion.
Governing equations and the numerical method Numerical simulations are performed here, by solving the complete compressible Navier-Stokes equations. The fourth order accurate piecewise parabolic method (PPM) together with the approximate Riemann Solver of Roe solve the inviscid equations. The overall scheme is second order accurate in time. The diffusion terms are regrouped in a source term of Navier-Stokes equations. Computations of high-speed jets are performed using a code specifically developed to numerically simulate gas flows whose composition may evolve with location ! x and time t. It is based on the usual set of gas flow equations, with the vector of the conserved variables is: Uð! x ; tÞ ¼ fr; rU1 ; rU2 ; rU3 ; re; rYi gT
(1)
r is the density, ðU1 ; U2 ; U3 Þ the three velocity components, e the total energy per unit mass, Yi is the mass fraction of the species i. The total energy of the flow is defined by: e¼
1 2 U þ U22 þ U23 þ Cv T 2 1
(2) Pns
where the average heat capacity per unit mass Cv ¼ i¼1 CviYi (ns is the number of involved species). Thermo-dynamics of the gas mixture is modelled with an ideal gas equation. The pressure P and the temperature T are related by: P RT ¼ r M
9541
Each component of Uð! x ; tÞ corresponds to a conserved quantity like mass, momentum or energy. The total energy e is the sum of the mechanical and internal energies of gas mixture. The equations for conservation of mass, momentum, energy and detailed mass of species i are written in conservative form. 3 vFj ðUð! x ; tÞÞ vUð! x ; tÞ X þ ¼ SðUð! x ; tÞÞ vt vxj j¼1
(4)
8 9 rUj > > > > > > > rUj U1 þ dj1 P > > > > > > > > < rUj U2 þ dj2 P > = ! Fj ðUð x ; tÞÞ ¼ rUj U3 þ dj3 P > > > > ðre þ pÞUj > > > > > > > > > > rYi Uj > > : ; … where: 8 > > > > > > > > > > > > > > > > > > > > > > > > > > > <
0 9 r > > > j¼3 > P v Re t1j > > > > > vx > j j¼1 > > > r > > > > j¼3 > P v Re t2j > > > > > vx > j j¼1 > > > = r t v SðUÞ ¼ j¼3 3j P Re > > > > > > > > vxj > > j¼1 > > > > > > > > > > > > > Cpr r ! > > > > > VT þ V$ U :t V$ > > > > Pr:Re Re > > > > > > > > > > r > > > > VY V$ > > i > > > > Sc:Re > > : ; … The right-hand side term of equation (4) contains terms of molecular flux (momentum, heat and mass). P3 vUk vUj 2 i tij ¼ vU þ d . k¼1 vxk vx vx 3 ij j
i
€ necker symbol, Re,Pr and Sc are the Reynolds dij is the Kro number, Prandtl number and Schmidt number of species i, respectively. They are defined in section High speed H2eAir jets. m, l, Di are the dynamical viscosity coefficient, heat conduction coefficient, species diffusion coefficient. Cp is the heat capacity at constant pressure of the mixture per unit mass. The values of Cp, Cv and M evolve and have to be calculated at each ð! x ; tÞ. The Cp value changes with the local composition of the H2eAir mixture and is dogged by a mass weighted over the heat capacities of the H2 and Air. The full compressible Navier-Stokes equation (4) are split up in three one-dimensional inviscid operators Lj ðj ¼ 1; 2; 3Þ: Lj :
! vUð! x ; tÞ vFj ðUð x ; tÞÞ ¼0 þ vt vxj
(5)
(3)
where R is the universal gas constant and M is the average molar mass calculated from the molar mass Mi of each species Ps 1 i: M ¼ ni¼1 Yi =Mi .
and three-dimensional viscous-diffusive operator: j:
vUð! x ; tÞ ¼ SðUð! x ; tÞÞ vt
(6)
9542
i n t e r n a t i o n a l j o u r n a l o f h y d r o g e n e n e r g y 4 0 ( 2 0 1 5 ) 9 5 3 9 e9 5 5 4
The fourth order accurate piecewise parabolic method (PPM) [16] together with the approximate Riemann solver of Roe [17] solve the equation (5). The spatial discretization of equation (6) is done by finite volume. The integration in time is accomplished by second order accurate low storage Runge-Kutta time stepping. The averaged values Ui over the ith grid cell are updated by computed the fluxes through the Riemann Solver between two discontinuities at each cell interface over one timestep. The space higher-order precision is achieved by taking into account the spatial variation of the discretized function in order to built two states UL (left) and UR (right) as an average of the information which is crossing to the interface over one timestep. The PPM algorithm generates a parabolic interpolation function for the primitive dynamic and scalar fields over each cell. The parabola is defined by its average over the cell and the border values which are computed with a fourth-order differencing formula. They are corrected where they generate new extrememums in order to realize monotonicity. The values of the primitive variables from both sides to each interface during one timestep, are acquired by averaging the piecewise parabolic profile over the appropriate domain of dependence. For a wave travelling with the speed c, the average is defined by: L ðvðxÞ; cÞ rightward waves : fiþ1=2 8 xiþ1=2 Z > > > 1 > < vðxÞdx if c > 0 ¼ cDt x cDt > iþ1=2 > > > : viþ1=2 if c 0
8 > > > > < 1 R cDt leftward waves : fiþ1=2 ðvðxÞ; cÞ ¼ > > > > :
(7)
Z
vðxÞdx viþ1=2
R
RLk ðxÞeLk ; VðxÞ Viþ1=2 ¼
k¼1
if c < 0
k¼3 X
RRk ðxÞeRk
(9)
(10)
k¼1
where Rk ðxÞ is the amplitude and ek ðxÞ the eigenvector of the linear wave associated with the kth characteristic. The Roe fluxes are determined after a Roe average by:
Fi±1=2 ¼
F ULi±1=2 þ F URi±1=2 2
(14)
The spatial discretization of these terms is done in finitevolume form. The fluxes and the stress tensor ! P3 vUk vUj 2 i t tij ¼ vU þ d are computed on the cell ink¼1 vxk vx vx 3 ij j
i
terfaces and are differenced conservatively. The combined solution of the three hyperbolic Lj ðj ¼ 1; 2; 3Þ and parabolic j problems are computed as follows: 1 1 1 1 Unþ1 ¼ j2 L1 L2 L3 j2 Un , Unþ1 ¼ j2 L3 L2 L1 j2 Un alternatively. Variables at time-step n þ 1 are computed from variables at time n. The timestep control is based on known criteria especially the CFL number which permits to obtain an individual time is step for convection operator Lj ðj ¼ 1; 2; 3Þ. Dtx < min clDx þcr the CFL number (cl , cr : left and right sound velocity). An optimal global value of the time step Dt is selected at each time inside the limit of stability to guarantee efficient computations. Dtmax ¼ minðDtx ; Dty ; Dtz ; Dtj Þ. The key param-
2Dmax
; Dmax ¼ max Di ; mr;
1
1 þ 1 þ 1 Dx2 Dy2 Dz2
l Cpr
.
if c 0
R
k¼3 X
(13)
xiþ1=2
L R Viþ1=2 ¼ fiþ1=2 ðVðxÞ; ui þ ci Þ ; Vi1=2 ¼ fi1=2 ðVðxÞ; ui ci Þ
L
(12)
/ 1 vr! u ¼ V $ $t vt Re vre ! Cpr ! ! 1 ! ¼ V$ VT þ V$ U :t vt Pr:Re Re
heat). Dtj <
xiþ1=2 cDt
The primitive vector of the dynamical and scalar fields on the left and on the right of the cells border, are decomposed into upstream and downstream travelling pressure wave characteristics. The left and right states of the primitive vector are linearised about:
VðxÞ Viþ1=2 ¼
vrYi ! r ! V Yi ¼ V$ Sc:Re vt
eters taken into account are the computational grid size, convection, and diffusion process (momentum, mass and !
(8)
L
The equations for species, momentum diffusion and energy conservation (sources terms of the equation (6)) are respectively:
k¼3 X a* R* e* k k k
(11)
k¼1
The scheme is of fourth order in time and in space where the solution is smooth and second order in space and third order in time near discontinuities.
Initial and boundary conditions In order to close the equation (4), boundary and initial conditions must be specified. In this free flow, boundary conditions behave as a model of the flow outside the computational domain. These conditions are given under a piecewise formulation. The boundaries are tailored in order to let mass fluxes; waves and flow structures freely escape the domain when mass fluxes only are permitted to enter the computation domain. Since the PPM-scheme integrates a differential equation along each ðx; tÞ-characteristics with an upwind scheme, the boundary conditions are written in a natural way. The inflow is sensitive to the upstream boundary condition with the mean shear originating the instabilities which develop downstream. A white noise is superimposed to a " !#! basic profile: uðrÞ ¼ Ua þ
Uj Ua 2
1 tanh
Rj 4q
Rj r
Rrj
in order
to favour unstable modes and model the residual turbulence of the upstream flow. Rj is the injector jet Raduis. A CroccoeBusemann's relation is also implemented to model the heating due to molecular dissipation as an effect of
i n t e r n a t i o n a l j o u r n a l o f h y d r o g e n e n e r g y 4 0 ( 2 0 1 5 ) 9 5 3 9 e9 5 5 4
9543
Fig. 2 e Test on outflow boundary condition.
ra compressibility rðrÞ ¼1þ
1rra j
2
" tanh
Rj 4q
!# Rj r
Rrj
! 1 .
The initial conditions have little physical significance although they are needed to start a computation. They can be somewhat arbitrary provided that the analysis of results starts at a time which is large enough to let phenomena develop with little dependence on initial conditions. The flow field is initialized over the whole domain with its upstream values. In cross stream and spanwise directions, free slip walls (or outflow without wall) and periodic boundary conditions are imposed respectively. At the downstream boundary, all the characteristics ðx; y; z; tÞ are transported outwards, to avoid the dependence between the flow and the outflow boundary. For the exit characteristic a Newmann condition is imposed: Unx þ1 ¼ Unx . For the incoming characteristic a Dirichlet condition is fixed to zero to avoid the perturbation of the upstream flow. The H2 central jet is subsonic, the outflow boundary should be efficient and consistently to avoid parasite reflections which can move back to upstream and then perturb the all the flow
Fig. 3 e Energy spectrum: computation with and without explicit subgrid model.
jets. To be cautiously and adequately controlled, a numerical test is performed as [18] which is based on innermost surface. The primitive dynamical variables are initialized by: 8 > > > > > > <
U1 ð! xÞ ¼ 0 ! U ðxÞ ¼ 0 2
U3 ð! xÞ ¼ 0 > > >
> > > : P ¼ P0 þ ε exp ln 2 r2 b2 As it is shown on Fig. 2, low reflection is generated; the sphere wave is submitted to negligible distortion, which gives confidence to our outflow boundary.
Academic test: temporal mixing layer A test configuration (temporal mixing layer) is analysed to justify and validate the approach (MILES) of the present heart article. Indeed, an explicit subgrid model based on filtered structure function [19] coupled to the flux limiters of the PPM
Fig. 4 e Isosurface of vorticity modulus: Grid resolution 64 64 64.
9544
i n t e r n a t i o n a l j o u r n a l o f h y d r o g e n e n e r g y 4 0 ( 2 0 1 5 ) 9 5 3 9 e9 5 5 4
Fig. 5 e H2 field in temporal mixing layer at different resolution: a) 48 48 48 b) 64 64 64 c) 96 96 96 d).128 128 128.
Fig. 6 e Energy spectrum at different grid resolutions.
Fig. 7 e Radial velocity profile: grid convergence test.
scheme is considered. The numerical implementation is accomplished on the basis of conventional LES for compressible flows approach. The Favre averaging is explicitly
~ ¼ rf. The filtered conserved vector is: taken into account by: f r ~ ~ ~ ~ i Þ. ~ U ¼ ðr; rU1 ; rU2 ; rU3 ; re; rY The filtered Navier-Stokes equations are written under this new conservative form.
i n t e r n a t i o n a l j o u r n a l o f h y d r o g e n e n e r g y 4 0 ( 2 0 1 5 ) 9 5 3 9 e9 5 5 4
Table 1 e CPU for different grid resolution. Grid resolution
CPU/Step/Node/Iteration
875 300 300 875 150 150 436 75 75
2.8271143E-06 (s) 2.2263755E-06 (s) 1.4447446E-06 (s)
Table 2 e Flow parameters, data from Eggers [12]. Parameter
Hydrogen jet
Air jet
1074 260=300 0:886 46:68 102 1:32 105 0:44 0.9 0.7
394 222=300 1:32
Uðm=sÞ Tstat =Ttot ðKÞ Mach rj U2j =ra U2a r U Djet Rejet ¼ jet mjet jet Ujet Uair Mc ¼ ajet þaair m ScH2 ¼ rH DH Cp m 2 2 Pr ¼ l H H 2 2
3 vFj ðUÞ vU X þ ¼ SðUÞ vt vxj j¼1
(15)
The filtered convective flux is: 8 9 > > > > > > ej > > rU > > > > > > > ~ > rUj U1 þ rRTdj1 > > > > > > < = ~ rUj U2 þ rRTdj2 Fj ðUÞ ¼ ~ j3 > > > rUj U3 þ rRTd > > > > > > > > > rðe þ RTÞU j > > > > > > > > rYi Uj > > : ; i h ~iU ~iU ~ j þ rUi Uj rU ~j rUi Uj ¼ rU |fflfflfflfflfflfflfflfflfflfflfflfflffl{zfflfflfflfflfflfflfflfflfflfflfflfflffl} Tij
h i ~ U ~ j þ rðe þ RTÞUj r e~ þ RT ~ U ~j rðe þ RTÞUj ¼ r e~ þ RT |fflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflffl{zfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflffl} Qij
i h ~ j þ rYi Uj rYi U ~j rYi Uj ¼ rYi U |fflfflfflfflfflfflfflfflfflfflfflfflffl{zfflfflfflfflfflfflfflfflfflfflfflfflffl} Gij
Tij , Qij , Gij are respectively the subgrid flux terms (momentum, heat, mass). They are modelled on the basis of eddy
9545
viscosity nt and eddy diffusivity Dt formulation. This model is built on proportionality between the subgrid stress and a local strain calculated from the large scale velocity field. In this case, eddy viscosity and diffusivity concepts are supported by physical arguments about damping of large scales by subgrid energy and k5/3 equilibrium spectra, nevertheless it misrepresent the so-called backscatter from the short waves to large ones and the interactions between short waves. The boundary conditions and the physical parameters are as follows: periodic in streamwise and spanwise. In transversal direction, an outflow condition is Imposed. The Reynolds number is about 10þ6, convective Mach number Mc ¼ 0:3. Four cases are analysed hereafter regarding the grid resolutions. Fig. 3 shows the turbulence spectrum with and without an additional physical subgrid scale model. The spectrum decay of the turbulent kinetic energy is satisfactory reproduced in the two cases. However, added an explicit subgrid model lead to worthlessness contribution and reduced performance. Furthermore, removing the limiters will affect the stability of the scheme as non-physical oscillations (Gibbs instabilities in the presence of strong gradients) grow, especially at low resolutions. The structure of the flow field is illustrated in the Fig. 4, the formation of coherent structures even at high Reynolds number. The flow topology resembles to that of Comte et al. [20]. In addition, the presence of vortex structures at small scales is obviously detected. Fig. 5 shows that PPM detect adequately, as the grid resolution is sufficiently improved, all the main compressible turbulent shear layer features including turbulent mixing. Increasing the resolution let to show more details and highly topology flow quality. The analysis of the spectra shows clearly and quantitatively its good tendency with increasing resolution (Fig. 6). In inertial range scales, a well characterized 5/3 power enveloppe over about one decade is visible which fill out the spectrum progressively by smaller scales as the grid resolutions is increased. These observations indicate that the large scales are well represented with a minimum dependence against the short one as the computational grid increases. This result corroborate the fact that this sophistical powerful scheme is eligible for the computation of high Reynolds compressible turbulent flows shocked or not, as nonlinear filter built-in the algorithm work adequately especially near the grid scale cut-off.
Fig. 8 e Cross section of H2 mass fraction on the jet axis.
9546
i n t e r n a t i o n a l j o u r n a l o f h y d r o g e n e n e r g y 4 0 ( 2 0 1 5 ) 9 5 3 9 e9 5 5 4
Fig. 9 e Cross section of vorticity field on the jet axis.
High speed H2eAir jets The numerical simulations were performed on 8 processor DELL Station with a total run time of 336 h. Several grid resolutions are tested and the appropriate one consists of 875 300 300 nodes, respectively in the streamwise, transverse and lateral directions. This choice is based on compromise between the available computation materials, the accuracy of steep dynamics and scalar gradients, robustness and significance of the results. The CPU are (see Table 1): A test convergence is carried out showing that the computations converge satisfactory to the experimental measurements as the resolution increases. Fig. 7 shows the progressive decay of the velocity as far as the grid resolution improves reducing thus the excess contribution of the unresolved scales. The target computation grid is then used in the analysis henceforth with confidence. The experiment carried out by Eggers [12] gives us an opportunity for comparing numerical simulations results and laboratory real flow scales. Computation parameters were designed to fit the real jets configuration. The jet at the injector exit is made of a Hydrogen at Mach number Mj ¼ 0.886, the adjacent jet is made of air at Mach number Ma ¼ 1.32. All gases are assumed to be ideal. Velocity, density and pressure are made non-dimensional by using the values
at the injector exit. Lengths are scaled on the diameter Dj of the hydrogen injector exit. Table 2 summarizes the flow parameters (the data from Eggers used in this study).
Results To qualitatively point up mixing with a set of pictures once the flow becomes statistically quasi-stationary and the effect of the initial condition has decayed, the instantaneous planar flow visualization is shown in Fig. 8. The natural quasi-linear region until x/Dj ¼ 5, where the instabilities grow freely, as the inlet perturbation is based on white noise, is obviously visible. The introduction of white noise perturbation avoids difficulties in specifying random perturbation with a given spectrum; this region is not representative for confrontation to experimental investigations as the causes of disturbances is linked to several sources like vibrations, pressure waves… However, it has been shown [3,4,21] that the consideration of hyperbolic profile let to develop easily the main instabilities and the flow features. For x/Dj > 7, the interaction of the two distinct mixing layers in the transition region, followed by a merging process into an inflate enclosure with a spate big spots of hydrogen. This enclosure deflates forthwith in the next nearer spatial positions which an indicative of constriction and relaxation
Fig. 10 e Visualisation of the potential core and rings vortex.
i n t e r n a t i o n a l j o u r n a l o f h y d r o g e n e n e r g y 4 0 ( 2 0 1 5 ) 9 5 3 9 e9 5 5 4
9547
Fig. 11 e Isosurfaces of Q criterion coloured by: (a) H2 masse fraction (b) streamwise velocity.
process. Downstream, one can see mixing over the entire thickness of the overall jet development. The hydrogen is sometimes dispersed top down, following the dynamics of the pulse enclosure. Spatially, the intermittent extent of the mixing is evident. The induction of the enclosure is also visible on the cut x-y vorticity field (Fig. 9) for 10 < x/Dj < 15. As for incompressible jet, the shear layer mode is evident; we can observe the development of structures on the two shear layers of the jet because of Kelvin-Helmoltz instabilities.
Beyond, eddies remain well organized and keep a high level of turbulence energy. The contrast between the vorticity field and mixing field is visible; the animations indicate different mechanisms. The potential cone is well visible in Fig. 10, the longitudinal and radial velocities remain constant near the centerline. The first structures (rings) appear at x/Dj ¼ 5 inducing a radial momentum diffusion otherwise velocity gradient. The vortex merging process linked to incompressible jet as observed by several investigations is not visible in our compressible
Fig. 12 e (a) Radial profile of the mean longitudinal velocity in the fully turbulent region (b) Radial profile of the rms of the longitudinal velocity component in the fully turbulent region.
9548
i n t e r n a t i o n a l j o u r n a l o f h y d r o g e n e n e r g y 4 0 ( 2 0 1 5 ) 9 5 3 9 e9 5 5 4
Fig. 13 e Cut through the (z,y) plane of the H2 mass fraction and vorticity fields.
situation. Furthermore, the varicose mode seems higher than the sinuous one as no helical structures are observed. Indeed, Fig. 11 shows an instantaneous view of threedimensional positive Q (is the second invariant of the velocity-gradient tensor, it is well known as a being a good indicator of turbulent structures) isosurfaces coloured by: H2 masse fraction (Fig. 11a) and streamwise velocity (Fig. 11b). The rings primary vortices appear in the early transition zone between the jet and the co-flow. Pairs of alternate streamwise vortices appear between two consecutive rings, in agreement with the classical phenomena of three-dimensionalization in
free-shear layers. The streamwise vortices play a dominant role in the transition processes towards a fully-turbulent state, and their location is very dependent of the inlet conditions [5]. The apparition of streamwise vortices is finally followed by an abrupt increase in the level of small-scale turbulence. Then the flow quickly reaches a state of fullydeveloped turbulence. A quantitative analysis is carried out, the mean and rms profiles are obtained through a temporal average of the dynamic and scalar fields. The Favre averaging is explicitly taken into account. The time averaging period is set to two times the
i n t e r n a t i o n a l j o u r n a l o f h y d r o g e n e n e r g y 4 0 ( 2 0 1 5 ) 9 5 3 9 e9 5 5 4
Fig. 14 e Downstream evolution: (a) azimuthal turbulent intensity (b) radial turbulent energy (continuous line) and azimuthal turbulent energy (dashed line).
Fig. 15 e H2 mass fraction isosurface (a) 0.07 (b) 0.2 (c) 0.5.
9549
9550
i n t e r n a t i o n a l j o u r n a l o f h y d r o g e n e n e r g y 4 0 ( 2 0 1 5 ) 9 5 3 9 e9 5 5 4
Fig. 16 e (a) Mean axial velocity profiles at several downstream sections. (b) rms axial velocity profiles at several downstream sections.(c) Mean radial velocity profile: numerical and experimental comparison.
flow time, where one flow-through time is defined as the time of the flow to traverse the computational domain. In order to compute the turbulence intensities, the relation: sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi r42 þ is used. The temporal average: 40 ¼ 42 24r4 r r Z fðx; y; zÞ ¼ T1
t0 þT t0
fðx; y; z; tÞdt. Where: fðx; y; zÞ ¼ Ui¼1;2;3 or YH2
T ¼ 4Lx =Uj is the time averaging period, t0 ¼ 4Lx =Uj is the time required to erase the effect of the initial condition. The mean and turbulent velocity profiles indicate a self similarity sate (Fig. 12). Moreover, there is a strong momentum transfer from jet to coflow. This transfer implies a spreading of the jet into co-flow. One can see high level of longitudinal turbulent intensity through the width of the jet and downstream position which is a consequence of the longitudinal vortices. However, the azimuthal instabilities are strong as the Reynolds number is high in our numerical simulations and substantially different to incompressible case. Indeed, Fig. 13 shows clearly the growing process of the azimuthal instabilities as several pairs of mushroom structures, which are responsible of diffusion, and dispersion of H2 on the jet periphery. These mushroom structures are located in the inner and outer rings. The inner ones shattered and eject H2 inside
the drained jet surface. However, the outer ones entrained air leading to the growth of mixture envelope. The inner can be due to the supersonic co-flow stream pushing the braids region within. Furthermore, Fig. 14 shows clearly an abrupt increase of azimuthal fluctuations in transition region because of a high mushrooms eddies activity. In the fully turbulent state, a decrease of these fluctuations is attributed to locking process of rings in enclosure followed by an alteration of the azimuthal rotation and the radial expansion of the rings. The production of intense ejections, in which mixture is expelled outwards from the rings are caused by the streamwise vortices, which also pump irrotational fluid (air), by entrainment. The estimated Strouhal number value is about 0.04 which obtained after a period of acquisition of the pressure field followed by a Fourier transform. The emergence of the streamwise vortices is followed by the creation of small turbulent structures where indicate the beginning of the transition to turbulence. The interaction and induction of inner and outer streamwise vortices lead to transition to turbulence which pump high level of air inside shear layer, The H2 is quasi completely dissipated as it is entrained outward the rings and transported by the irrotational stream (air). Indeed, a low values of H2 species as well as vorticity are clearly observed for x/Dj > 20. To overcome this situation, it is then recommended to lock the inner jet by an
i n t e r n a t i o n a l j o u r n a l o f h y d r o g e n e n e r g y 4 0 ( 2 0 1 5 ) 9 5 3 9 e9 5 5 4
9551
Fig. 17 e (a) Mean radial H2 mass fraction profiles at several downstream sections. (b) rms radial H2 mass fraction profiles at several downstream sections.(c) rms axial H2 mass fraction profile (d) Mean radial H2 mass fraction profile: numerical and experimental comparison.
outer jet (replacing the co-flow jet) to avoid the dissipation of H2 in turbulent zone. This recommendation is inspired from the work of Balarac and Si-Ameur [5] who found that the inner jet remain confined by the outer one in a configuration of incompressible coaxial jets. Fig. 15 shows iso-mixture fractions which are valuable for fractal dimensions and mixing efficiency. Depiction of the isosurface of H2 is important for modelling turbulent mixing. Thus, the turbulent mixing starts with an engulfment of two species through the shear layers implied by the KelvinHelmholtz vortices. The isosurface is highly convoluted after the amplification instabilities, with evidence of a mixing envelope in these jets for low value of H2 Fig. 15a, unlike the high value
where spatial intermittence is palpable Fig. 15b and c. This result show clearly that H2 is highly diffused and dispersed far from the central jet until his complete quasi dissipation. This disadvantage is due to significant momentum and mass transfer by the longitudinal eddies, towards the coflow jet. This observation constitutes an important stage in the evolution of H2-air mixing. This finding is useful for combustion process especially in the phase of flame ignition as molecular mixing is primordial to ensure a normal ignition process. The Fig. 16a shows the radial variation of the mean streamwise velocity profile at several downstream positions. At the beginning of transition, the jet maintains a profile analogous to the inlet boundary profile. Continuous decreases
Fig. 18 e Transverse jet section coloured by probability density of H2.
9552
i n t e r n a t i o n a l j o u r n a l o f h y d r o g e n e n e r g y 4 0 ( 2 0 1 5 ) 9 5 3 9 e9 5 5 4
Fig. 19 e Mean radial (velocity and H2 mass fraction) profiles: numerical and experimental comparison a) x/D ¼ 9.58, b) x/ D ¼ 15.44, c) x/D ¼ 25.2.
of the maximum mean velocity along the jet axis, followed by a spreading of jet into the co-flow as a consequence of momentum transfer, are evident and similar to classical incompressible jets. However, the turbulent activity is located at the interface at the beginning of transition (Fig. 16b), then it spread on the total width. Two peaks are observed which is consistent with the visualizations of two mushrooms structures (sublayers inner and outer). A continuous increase is located in the enclosure followed by abrupt decrease of streamwise turbulent intensity. Our observations are
corroborating by the good agreement with the measurements of Eggers, as it is evident on Fig. 16c. The Fig. 17a illustrates the average radial profile of the H2 at with various sections of the jet. Two stages can be distinguished; upstream the H2 profiles evolves gradually; the increases in mixture zone is due to radial and longitudinal evolutions noticeably the pulsatile movement of the ring structures permitting the penetration of air in H2 jet; the rms of the mixing fraction is very low in this region (Fig. 17b). Two peaks are also evident which are the signature of the two
i n t e r n a t i o n a l j o u r n a l o f h y d r o g e n e n e r g y 4 0 ( 2 0 1 5 ) 9 5 3 9 e9 5 5 4
9553
Fig. 20 e 2D computations: a) H2 field, b) pressure field.
sublayers mushrooms structures. Fig. 17c highlights a continuous increase followed by abrupt decrease of the H2 rms, indicating blockage of radial and azimuthal rings development (pulsating movement is altered in the enclosure). It follows longitudinal ejections which only transfer the H2 in the co-flow and inducing its quasi dissipation. The comparison between the experimental measurements of Eggers (1971) and the present numerical computations (Figs. 16c and 17d) shows a good agreement qualitatively and quantitatively. The turbulent mixing activity is confined in narrow thickness jr/Djj<1 in transition region, however it extends downstream over twice, as it is visible on Fig. 16a. The figure shows that a significant quantity of air invaded the jet center. There is thus a strong transfer of H2 to ambient jet towards the central jet implied by the transfer of momentum that we already evoked. The study of the fluctuations of H2 (Fig. 17b) reveals good tendency compared to experimental data. Moreover, even the small values of H2 mainly within the periphery of the jet downstream, a considerable scalar fluctuations, on the axis of the jet, are detected. That is explained by a strong intermittency of the turbulent mixing activity, as it is shown on Fig. 18, which presents a transverse jet section coloured by probability density of H2. The probability to find a high mass fraction of H2 is small on all the jet width, which shows clearly his quasi-dissipation mentioned above. Furthermore, the strong intermittency in the area of the central jet, undoubtedly due to the injections implied by the longitudinal vortices. A further confrontation of the present computations against experimental at different sections have revealed good agreement. Fig. 19 shows clearly this the good performance of the PPM to track dynamics and mixing in high-speed shear flows comparatively to WENO 5 scheme [15].
Confinement effect All the results of the subsection Results are obtained in confinement environment. However, another study has been undertaken by removing the walls and adopting outflow boundary conditions in transverse direction. The exact same results are obtained for unconfined situation which means that the walls do not play any role in the dynamic development of the jets. This can be explained by the value of the
convective Mach number about 0.44. The wall effects can be important when Mc > 1 as it is demonstrated by Si-Ameur & Chollet [22]. An examination of the 2D computations has revealed the same finding. Indeed the Kelvin-Helmholtz instabilities, the interaction of the large turbulent structures with the mixing process are unchanged in the two situation. The interaction of the mixing zone with eventually Mach wave is insignificant in the case of the confinement (Fig. 20).
Conclusion A three-dimensional spatially developing jets, of two supersonic-subsonic streams of Hydrogen and air gas has been performed using the Monotone Integrated Large eddy simulation approach. The computations explicitly include the near field and the far field of the jets development with an emphasis on mixing. A locking process of rings in enclosure followed by an alteration of their azimuthal rotation and radial expansion are observed, through flow visualizations and statistics. Furthermore, H2 is highly diffused and dispersed far from the central jet until his complete quasi-dissipation. After the transition, the flow structure is directed preferentially in the streamwise direction due to the longitudinal turbulent scales, which are ejected from the enclosure. These structures particularly entrain air stream continuously inside jet leading to quasi-complete alteration of H2. The examination of the transverse jet coloured by PDFs has revealed the existence of strong intermittency in the area of the central jet, undoubtedly due to the injections implied by the longitudinal vortices. A continuous increase followed by abrupt decrease of the H2 rms, indicating blockage of radial and azimuthal rings development (pulsating movement is altered in the enclosure). It follows longitudinal ejections which only transfer the H2 in the co-flow and inducing its near dissipation. The confined environment do not affect the hydrogen-air mixing mechanism as the convective Mach number is subsonic. A comparison with laboratory experiment has shown a good agreement as well as for the mean and turbulent
9554
i n t e r n a t i o n a l j o u r n a l o f h y d r o g e n e n e r g y 4 0 ( 2 0 1 5 ) 9 5 3 9 e9 5 5 4
dynamic and scalar fields. On se basis of this study, the mixing aireH2 jets of this configuration is not suitable for combustion application, as the fuel can leave completely the central jet and ejected outflow by the supersonic co-flow. It is recommended to lock the inner jet by an outer jet to avoid this probable situation.
references
[1] Villermaux E, Rehab H. Mixing in coaxial jets. J Fluid Mech 2000;425:161e85. [2] Ritchie BD, Mujumdar DR, Seitzman JM. Mixing in coaxial jets using synthetic jet actuators. 2000. AIAA Paper. tais O, Lesieur M. Direct numerical [3] Balarac G, Si-Ameur M, Me simulations of high velocity ratio coaxial jets: mixing properties and upstream conditions influence. J Turbul 2007;8(22). tais O, Lesieur M. Large Eddy [4] Balarac G, Si-Ameur M, Me simulation of coaxial jets: coherent structures and mixing properties. Direct Large Eddy Simul VI 2006;10:277e84. [5] Balarac G, Si-Ameur M. Mixing and coherent vortices in turbulent coaxial jets. C R Mec 2005;333:622e7. [6] Reynier P, Minh HH. Numerical prediction of unsteady compressible turbulent coaxial jets. Comput Fluids 1998;27(2):239e54. [7] Klioutchnikov I, Olivier H, Odenthal J. Numerical investigation of coaxial jets entering into a hot environment. Comput Fluids 2013;86:490e9. [8] Jin Tai, Luo Kun, Lu Shuqiang, Fan Jianren. Direct numerical simulation of a supersonic lifted hydrogen jet flame: a priori study on combustion models. Acta Astronaut 2015;109:52e64. [9] Luo Kun, Jin Tai, Lu Shuqiang, Fan Jianren. DNS analysis of a three-dimensional supersonic turbulent lifted jet flame. Fuel 2013;108:691e8.
[10] Moule Y, Sabelnikov V, Mura A. Highly resolved numerical simulation of combustion in supersonic hydrogeneair coflowing jets. Combust Flame 2014;161:2647e68. [11] Si Ameur M. Isothermal reactive mixing layer: numerical study. Comput Therm Sci 2011;3:227e35. [12] Eggers JM. Turbulent mixing of coaxial compressible hydrogeneair jets. Tech. Rep. NASA-TN D-6487. 1971. [13] Porter DH, Pouquet. Kolmogorov-like spectra in decaying three-dimensional supersonic flows. Phys Fluids 1994;6:2133e42. [14] Porter DH, Pouquet A, Woodward PR. A numerical study of supersonic turbulence. Theor Comput Fluid Dyn 1992;4:13e49. [15] Karaca M, Lardjane N, Fedioun I. Implicit large eddy simulation of high speed non reacting and reacting air/H2 jets with a 5th order WENO scheme. Comput Fluids 2012;62:25e44. [16] Colella P, Woodward PR. The piecewise parabolic method (PPM) for gas-dynamical simulations. J Comput Phys 1984;54:174e201. [17] Roe PL. Approximate Riemann solvers, parameter vectors and difference schemes. J Comput Phys 1981;43:357e72. [18] Lodato G, Domingo P, Vervisch L. Three-dimensional boundary conditions for direct and large-eddy simulation of compressible viscous flows. J Comput Phys 2008;227:5105e43. [19] Ducros F, Comte P, Lesieur M. Large-eddy simulation of transition to turbulence in a boundary layer developing spatially over a flat plate. J Fluid Mech 1996;326:1e36. [20] Comte P, Lesieur M, Lamballais E. Large and small-scale stirring of vorticity and a passive scalar in a 3D temporal mixing layer. Phys Fluids A 1992;4:2761e78. rique des effects de tempe rature dans [21] Daviller G. Etude nume les jets simples et coaxiaux. PhD thesis, ISAE-ENSMA Poitiers. 2010. [22] Si Ameur M, Chollet JP. Transition to turbulence of supersonic mixing layer. C R Acad Sci Paris Iib 1998;326:315e22.