On the role of bulk viscosity in compressible reactive shear layer developments

Tenth International Conference on Computational Fluid Dynamics (ICCFD10), Barcelona, Spain, July 9-13, 2018 ICCFD10-2018-0086 On the role of bulk viscosity in compressible reactive shear layer developments Radouan Boukharfane∗, Pedro J. Martínez Ferrer ∗, Arnaud Mura∗ ∗PPRIME UPR 3346 CNRS, ENSMA, 86961 Futuroscope Chasseneuil Cedex, France Vincent Giovangigli∗∗ ∗∗ CMAP UMR 7641 CNRS, Ecole Polytechnique, 91128 Palaiseau Cedex, France Corresponding author: arnaud.m[email protected] Abstract Despite 150 years of research after the reference work of Stokes, it should be acknowledged that some confusion still remains in the literature regarding the importance of bulk viscosity effects in flows of both academic and practical interests. On the one hand, it can be readily shown that the neglection of bulk viscosity (i.e., κ= 0) is strictly exact for mono-atomic gases. The corresponding bulk viscosity effects are also unlikely to alter the flowfield dynamics provided that the ratio of the shear viscosity µto the bulk viscosity κremains sufficiently small. On the other hand, for polyatomic gases, the scattered available experimental and numerical data show that it is certainly not zero and actually often far from negligible [13]. Therefore, since the ratio κ/µ can display significant variations and may reach very large valuesa, it remains unclear to what extent the neglection of κ holds [3]. The purpose of the present study is thus to analyze the mechanisms through which bulk viscosity and associated processes may alter a canonical turbulent flow. In this context, we perform direct numerical simulations (DNS) of spatially-developing compressible non-reactive and reactive hydrogen-air shear layers interacting with an oblique shock wave. The corresponding flowfield is of special interest for various reactive high-speed flow applications, e.g., Scramjets. The corresponding computations either neglect the influence of bulk viscosity (κ= 0) or take it into consideration by evaluating its value using the EGlib library [15]. The qualitative inspection of the results obtained for two-dimensional cases in either the presence or the absence of bulk viscosity effects shows that the local and instantaneous structure of the mixing layer may be significantly altered when taking bulk viscosity into account. This contrasts with some mean statistical quantities, e.g., the vorticity thickness growth rate, which do not exhibit any significant sensitivity to the bulk viscosity. Enstrophy, Reynolds stress components, and turbulent kinetic energy (TKE) budgets are then evaluated from three-dimensional reactive simulations. Slight modifications are put into evidence on the energy transfer and dissipation contributions. From the obtained results, one may expect that refined large-eddy simulations (LES) may be rather sensitive to the consideration of bulk viscosity, while Reynolds-averaged Navier-Stokes (RANS) simulations, which are based on statistical averages, are not. The filtering of the present dataset may provide further insights so as to assess (or not) such a conclusion. Keywords: Bulk Viscosity, Shear Layer, Direct Numerical Simulation, Molecular Transport aIt can exceed thirty for dihydrogen. 1 1 Introduction The bulk (or volume) viscosity κ, which can related to the second (or dilatational) viscosity coefficient λ, is related to the vibrational and rotational energy of the molecules. From the macroscopic viewpoint, it characterizes the resistance to dilatation of an infinitesimal bulk element at constant shape [2]. It is strictly zero only for dilute monoatomic gases and this theoretical result is often used to discard it, regardless of the nature or internal structure of the fluid as well as the flowfield conditions. However, acoustic absorption measurements performed at room temperature have shown that the ratio of the volume to the shear viscosity κ/µ may be up to thirty for dihydrogen [10], and recent analyses of reactive multicomponent high-speed flows have confirmed that it is not justified to neglect it, except for the sake of simplicity [3]. The dilatational viscosity is important in describing sound attenuation in gaseous media, and the absorption of sound energy into the fluid depends itself on the sound frequency, i.e., the rate of fluid expansion and compression. For polyatomic gases, the available measurements of κ, which remains quite seldom due to the complexity of its determination, show that it is certainly not zero and actually far from negligible. It is also noteworthy that theoretical analyses do show that κ/η is at least of the order of unity. Therefore, since the ratio κ/µ can display significant variations and may reach very large values, it is unclear to what extent the Stockes hypothesis (i.e., λ=−2µ/3or κ= 0) may hold for compressible and turbulent flows of gases featuring a ratio κ/µ greater than unity. In either an expansion or a contraction of the gas mixture, the work done by the pressure modifies immediately the translational energy of the molecules, while a certain time-lag is needed for the translational and internal energy to re-equilibrate through inelastic collisions [8]. This can be described through a system of two coupled partial differential equations written for the internal and translational temperatures, with a pressure-dilatation term that acts as a source term in the translational temperature budget. The volume (or bulk) viscosity is associated to this relaxation phenomenon and it is evaluated from this internal energy relaxation time-lag. The evaluation of this property for a mixture of polyatomic gases is far from being an easy task since the kinetic theory of gases does not yield an explicit expression for this transport coefficient, but instead linear systems that must be solved [14]. The corresponding systems are derived from polynomial expansions of the species’ perturbed distribution functions. The bulk viscosity is obtained here using the library EGlib developed by Ern and Giovangigli [13,15]. It is evaluated as a linear combination of the pure species volume viscosities, which require the evaluation of various collision integrals [14]. The impact of bulk viscosity effects has been previously analysed in several situations including shock-hydrogen bubble interactions [3], turbulent flames [17], compressible boundary layers [11], shock-boundary layer interaction [1], and planar shock-wave [9]. All these studies confirm that the bulk viscosity effects may be significant. The purpose of the present work is to assess its influence in regard to both the instantaneous and statistical features of canonical compressible turbulent multicomponent flows. Using direct numerical simulation (DNS), we investigate the impact of the bulk viscosity coefficient κon the spatial development of reactive and non-reactive compressible mixing layers interacting with an oblique shock wave. Such a canonical flowfield is typical of the shock-mixing layer interactions that take place in compressible flows of practical interest. For instance, supersonic jets at high nozzle-pressure ratio (NPR) give rise to a complex cellular structures, where shocks and expansions waves interact with the turbulent outer shear layer [7]. It is also encountered in Scramjet intakes and combustors, where shock waves interact with the shear layers issued from the injection systems. On 2 the one hand, it is clear that the occurrence of shock waves in supersonic combustors induces pressure losses that cannot be avoided but, on the other hand, the resulting shock interactions with mixing layers contribute to scalar dissipation (i.e., mixing) rates enhancement [4], and may favor combustion stabilization in high-speed flows. The present manuscript is organized as follows: the mathematical model is presented in the next section (i.e., §2), which also includes a short description of the numerical methods. The details of the computational setup are subsequently provided in section §3. Section §4gathers all the results issued from (i) two-dimensional numerical simulations of both inert (§4.1) and reactive (§4.2) cases, and (ii) the three-dimensional case, which is analysed in §4.3. Finally, some concluding remarks and perspectives for future works are presented in section §5. 2 Mathematical description and computational model In this work, the in-house massively parallel DNS solver CREAMS is used. It solves the unsteady, three-dimensional set of compressible Navier-Stokes equations for multicomponent reactive mixtures [23]: ∂t(ρ) + ∇·(ρu)=0,(1a) ∂t(ρu) + ∇·(ρu⊗u) = ∇·σ,(1b) ∂t(ρEt) + ∇·(ρuEt) = ∇·(σ·u−J),(1c) ∂t(ρYα) + ∇·(ρuYα) = −∇·(ρVαYα) + ρ˙ωα,(1d) where tdenotes the time, ∇is the spatial derivative operator, uis the flow velocity, ρis the density, Et=e+u·u/2is the total specific energy (obtained as the sum of the internal specific energy, e, and kinetic energy), Yαis the mass fraction of chemical species α(with α∈S={1,...,Nsp}), Vαis the diffusion velocity of species α,Jthe heat is flux vector and ˙ωαrepresents the chemical production rate of species α. The integer Nsp denotes the number of chemical species. The above set of conservation equations (1) requires to be completed by constitutive laws. In this respect, the ideal gas mixture equation of state (EoS), P=ρRT/Wwith Rthe universal gas constant, is used to relate the pressure Pto the temperature T. In this expression, the quantity Wdenotes the molar weight of the multicomponent mixture, which is obtained as the sum of the molecular mass of each individual species W−1=PNsp α=1 Yα/Wα. Within the framework of the kinetic theory of dilute polyatomic gas mixtures, the molecular diffusion velocity vector Vα,α∈S, heat flux vector J, and second-order stress tensor σare expressed as follows: ρVαYα=−X β∈S ρYαDα,β (dβ+χβXβ∇( log T) ) ,(2a) J=X α∈S ρVαYαhα+RTχα Wα−λT∇T, (2b) σ=−PI+τ=−PI+µ(∇u+∇u|) + λ(∇·u)I,(2c) where Dα,β,(α, β )∈S2, are the multicomponent diffusion coefficients, dα,α∈S, the species diffusion driving forces, χα,α∈S, the rescaled thermal diffusion ratios, Xα,α∈S, the species mole fractions, hαthe enthalpy per unit mass of the α-th species, and λTthe thermal conductivity. The diffusion driving force dαof the α-th species is given by dα= 3 ∇Xα+ ( Xα−Yα)∇( log P). The quantity µdenotes the shear viscosity and λdenotes the second (or dilatation) viscosity coefficient. The bulk viscosity coefficient κappears explicitly in the expression of the viscous stress tensor τ. A relationship between the bulk viscosity κand viscosity coefficients µand λcan be deduced from the expression of the total pressure, which can be evaluated as the component of the spherical tensor based on the trace of the total stress tensor σ: −tr(σ) 3=− i=3 X i=1 σii 3=P−λ+2 3µ∇·u=P−κ∇·u(3) The second term in the right-hand-side of the above expression is the dilational contribution, which defines the bulk viscosity as κ=λ+2µ/3. As mentioned above, the Stokes’ hypothesis, stating that λ=−2µ/3(and hence κ= 0), is often retained as a simplifying assumption. Many efforts have been devoted to the derivation of relationships between the bulk viscosity and fundamental fluid properties [21,29]. If we consider a single polyatomic gas with a unique internal energy mode, the internal energy relaxation time τint can be related to the bulk viscosity [8,6]: κ=PR/c2 v·cint τint,(4) where cint denotes the internal heat capacity and cvthe specific heat at constant volume. When there are several internal energy modes and/or several species present in the mixture, the above simple expression is replaced by the solution to a linear system [12]. Within the Monchick and Mason approximation [26], neglecting complex collisions characterized by more than one quantum jump, the reduced system is diagonal and yields κ[3]: κ=PR/c2 v·X k∈P Xkcint kτint k,(5) where P= 1,· · · , npis the polyatomic species indexing set. The average relaxation time for internal energy of the k-th species τint kis then expressed as: cint k/τint k=X m∈N cm k/τm k,(6) where τm kdenotes the average relaxation time of internal energy mode mfor the k-th species, and Nis the internal energy mode indexing set. The CREAMS solver is coupled with the EGlib library to estimate transport coefficients from the kinetic theory of gases [16]. In this library, the optimized subroutines EGSKmare used to evaluate the bulk viscosity. The integer m∈J2,6Kassociated to the subroutine name refers to retained level of approximation. The higher the value of m, the more expensive the algorithm but also the more accurate the bulk viscosity expression. Following the work of Billet et al. [3], the value m= 3 is retained for the purpose of the present study. The shear viscosity and difusion velocities are evaluated with the routines EGFE3 and EGFYV, respectively. EGFLCT3 is used to determine the thermal conductivity λTand rescaled thermal diffusion ratios χα. The above system (1) is discretized on a Cartesian grid. A seventh-order accurate WENO scheme is used to approximate inviscid fluxes, while an eighth-order accurate centered difference scheme is retained to approximate viscous and diffusive contributions. Time integration is performed with a third-order accurate TVD Runge-Kutta scheme. The stiffness associated 4 to the wide range of time scales involved in the description of the chemical system is addressed using the Sundials CVODE solver [20]. A standard splitting operator technique, similar to the one previously retained in reference [32], is used. β U1 U2 U1 U2 x2 x1 @R Mixing layer  Transmitted shock Shock generator @I Computational domain @I Incident shock @R Reflected shock Figure 1: Sketch of the twodimensional shock–mixing layer interaction geometry. The computational domain dimensions are Lx1×Lx2= 275.0×120.0 in inlet vorticity thickness units (i.e., δω,0). It is uniformly discretized using Nx1×Nx2= 1640 ×720 grid points. 3 Problem statement and computational setup We study the interaction of an oblique shock with a spatially-developing shear layer. The upper stream corresponds to the fuel inlet, i.e., a mixture containing hydrogen, and the bottom inlet stream to vitiated air. Both twoand three-dimensional computations are performed. Figure 1 provides a typical sketch of the corresponding computational geometry and Table 1 gathers the values of the main parameters relevant to the present numerical simulation. The flow initialization is similar to the one retained in reference [23]. Assuming equal free-stream specific heat capacity ratios, the convective Mach number may be evaluated from Mc= (U1− U2)/(a1+a2), where a1and a2denotes the sonic speeds of streams 1 (oxidizer inlet stream) and 2 (fuel inlet stream) respectively. For the present set of computations, it is equal to Mc= 0.48. Fuel Oxidizer T(K) 545.0 1475.0 u1(m/s) U2U1 u2(m/s) 0.0 0.0 u3(m/s) 0.0 0.0 ρ(kg/m3)0.354 0.203 YH2(−)0.05 0.0 YO2(−)0.0 0.278 YH2O(−)0.0 0.17 YH(−)0.0 5.60 ·10−7 YO(−)0.0 1.55 ·10−4 YOH (−)0.0 1.83 ·10−3 YHO2(−)0.0 2.50 ·10−7 YN2(−)0.95 0.55 Table 1: Parameters of the shock–mixing layer interaction case. The mixing layer flow is impinged by an oblique shock wave that is issued from the oxidizer inlet stream (1) at the bottom boundary. The oblique shock wave angle is β= 33◦, see Figure 1. The geometrical parameters relevant to the present set of numerical simulations are provided in Table 2. The quantities Lx1,Lx2, and Lx3denote the computational domain lengths in each direction normalized by the initial vorticity thickness δω,0, while Nx1,Nx2, and Nx3are 5 the corresponding numbers of grid points. In the two-dimensional computations, only the x1and x2-directions are considered. Table 2: Computational mesh description. Lx1Lx2Lx3Nx1Nx2Nx3δω,0(m) 280 130 15 1640 750 180 1.44e−4 The flow is initialized with a hyperbolic tangent profile for the streamwise velocity component, while the other velocity components are set at zero. Species mass fractions and density are also set according to the following general expression: ϕ(x1, x2, x3) = ϕ1+ϕ2 2+ϕ1−ϕ2 2tanh 2x2 δω,0,(7) where ϕdenotes any of the flow variables mentioned above (i.e., species mass fraction or streamwise velocity component). The value of the Reynolds number Reω, based on the initial vorticity thickness and inlet velocity difference ∆U = U1−U2is Reδω= 640. Dirichlet boundary conditions are applied at the two supersonic inlets, perfectly non-reflecting boundary conditions are set at the outflow, and periodic boundary conditions are settled along the x3-direction. A slip boundary condition is imposed at the top, while the bottom boundary condition is set by using Rankine-Hugoniot relations, generalized for a multicomponent mixture [25]. In order to trigger flow transition, a slight white noise fluctuation is superimposed to the transverse velocity component along the line (x1, x2) = (4δω,0,0). The value of the CFL number is set to 0.75. Reactive flow simulations are conducted with the detailed mechanism of O’Conaire et al. [27]. It consists of nine chemical species (H2, O2, H2O, H, O, OH, HO2, H2O2, and N2) and 21 elementary reaction steps. The concentrations of these species at the inlet have been determined from equilibrium conditions so as to reach favorable self-ignition conditions within the extension of the computational domain. u1(m/s) u2(m/s) T(K) P(Pa) ρ(kg/m3) Fuel 1634.0 0.0 1475.0 94232.25 0.354 Oxidizer 1526.0 156.7 1582.6 129951.6 0.421 Bottom 973.0 0.0 545.0 94232.25 0.203 Table 3: Flow parameters of the shock–mixing layer interaction. Throughout this manuscript, the Reynolds and Favre averages of any quantity ϕare denoted by ϕand eϕ, while the corresponding statistical fluctuations are denoted by single and double primes, i.e., ϕ0and ϕ00, respectively. Averaging is performed over both transverse directions of statistical homogeneity (x1and x2) and time t. 4 Analysis of computational results 4.1 Inert two-dimensional mixing layer Figure 2 displays the instantaneous field of the ratio κ/µ computed for the present flow conditions. From this figure, it is noteworthy that (i) this ratio reaches values significantly larger than unity, (ii) it exhibits important spatial variations, the most significant of which are related to mixture composition. This contrasts with its sensitivity to pressure variations (i.e., 6 shock waves), which seems to remain rather moderate. Considering the values of κ/µ, as well as the amplitude of its variations, one may expect some remarkable effects of the bulk viscosity on this inert flowfield. It is the objective of this preliminary section to study to what extent the bulk viscosity may influence the instantaneous and statistical characteristics of the two-dimensional mixing layer development. Figure 2: Instantaneous field of the ratio κ/µ in the case κ6= 0 at t∆U/δω,0= 75.0. Instantaneous flow visualizations are very revealing of some local features of the shear layer, which are filtered out once fields or cross-stream profiles of averaged quantities are considered instead. For instance, quantitative comparisons of the onset of the streamwise vortices formation can be obtained from the instantaneous fields of the dimensionless magnitude of the density gradient, i.e., “numerical Schlieren”, reported in Figure 3. (a) κ6= 0 (b) κ= 0 Figure 3: Instantaneous field of the numerical density-based Schlieren at t∆U/δω,0= 75.0. In this figure, it is remarkable that, in the absence of bulk viscosity, the normalized abscissa at which the vortex roll-up processes take place is approximately x1/δω,0= 103.0, while in the situation featuring non-zero value of the bulk viscosity, it can be estimated to be x1/δω,0= 116.0. This can be explained by the same argument as the one invoked by Billet et al. [3] in their study of shock/hydrogen bubble interaction. The shear layer is also a diffusion layer associated to a density gradient, the absence of bulk viscosity makes the baroclinic term ∇P×∇ρ/ρ2greater at the shock / mixing layer interaction location, thus favoring the birth of velocity fluctuations. When volume viscosity is taken into account, the shock is much smoother in agreement with the physical theory of shock wave internal structure. As a consequence, the baroclinic production term is lower than in absence of volume viscosity. The interaction between the mixing layer and the shock wave can be further assessed by considering the Richtmyer-Meshkov instability development, the key point of which is the baroclinic 7 effect [5]. The Richtmyer-Meshkov instability indeed takes place when two fluids having different densities are impulsively accelerated, similarly to what occurs when the shock wave impinges the mixing layer in the present study. This process may be analysed by considering the transport equation for the enstrophy Ω = |ω|2/2. Such a transport equation is readily obtained by (i) taking the curl of the momentum transport equation and subsequently (ii) multiplying each term of the resulting equation by the vorticity vector itself: DΩ Dt=∂tΩ + u·∇Ω =1 2ω·(∇u+∇u|)·ω | {z } E −Ω∇·u | {z } D +ω·∇P×∇ρ ρ2 | {z } B +ω·∇×∇·τ ρ | {z } V (8) The production / destruction terms on the right hand side (RHS) of (8) are associated to vortex-stretching, dilatation, baroclinic torque, and viscous dissipation. The baroclinic contribution ∇P×∇ρplays a significant role in the enstrophy and vorticity production and it appears as one of the main sources of vorticity in supersonic flows [31]. (a) κ6= 0 (b) κ= 0 Figure 4: Instantaneous field of the mixing zone colored by k∇P×∇ρk. Figure 4 displays the instantaneous field of the magnitude of ∇P×∇ρwithin a region restricted to mixture fraction values such that ξ∈[0.05,0.95]. This passive scalar ξis defined on the basis of conserved elemental (i.e., atomic) mass fractions [28]. The mass fraction of chemical element γ, denoted aγ, is readily deduced from the chemical species mass fractions: aγ= Nsp X α=1 YαNα,γAγ Wα ,(9) where Aγis the atomic weight associated to element γand Nα,γ denotes the number of γ atoms present in each molecule of chemical species α. For a two-feeding inlet system such as 8 the one considered here1, the mixture fraction is then obtained by summing over all elemental mass fractions and normalizing the result ξ=Pγ|aγ−aγ,O| Pγ|aγ,F−aγ,O|,(10) where aγ,Oand aγ,Fdenote the mass fractions of atom γin the oxidizer and fuel inlet streams, respectively. Figure 4 shows that the term ∇P×∇ρdisplays large values at the periphery of the mixing layer where the pressure gradient and the density gradient are significantly misaligned, thus promoting the development of the mixing layer. The misalignment of the pressure gradient – that is imposed by the shock wave – and local density gradient – that is associated to the mixing layer – serves as a basis to vorticity generation through the baroclinic term. It may contribute significantly to mixing enhancement. 0.2 0.3 0.4 0.5 0.6 0.7 0.8 20 40 60 80 100 || P× || 0.2 0.3 0.4 0.5 0.6 0.7 0.8 20 40 60 80 100 0.0 20.8 41.7 62.5 83.4 104.2 125.0 145.9 166.7 187.6 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.2 0.3 0.4 0.5 0.6 0.7 0.8 20 40 60 80 100 20 40 60 80 100 20 60 100 140 180 ξ ξ k∇P×∇ρk Figure 5: JPDF of ξand k∇P×∇ρkfor κ6=(left) and κ= 0 (right). To quantify the vorticity production that is induced by the baroclinic term in the presence or in the absence bulk viscosity, Figure 5 reports the joint probability density function (JPDF) of the passive scalar ξand baroclinic term k∇P×∇ρkobtained for both cases. This figure reveals that the production of vorticity by baroclinic effects is concentrated around the stoichiometric mixture fraction value ξst = 0.43. It also confirms that the neglection of bulk viscosity tends to enhance large values of this production term. One of the fundamental statistical quantities that characterizes the mixing layer development is its normalized growth rate [30]. Although the definition of this growth rate is not unique, its most standard expression relies on the vorticity thickness definition: δω(x1) = U1−U2 ∂eu1/∂x2|max (11) Figure 6 displays the spatial evolution of normalized vorticity thickness for both cases. It reveals that the interaction of the reflected shock wave with the mixing zone changes significantly the mixing layer growth rate (i.e., the slope) and that the evolution of the mixing layer 1Other situations featuring more than two inlets have been recently addressed by Gomet et al. [18]. 9 Table 4: Values of (U1/U2)dδω/dx1obtained from the three-dimensional mixing layer computations. 1st region 2nd region 3rd region κ= 0 0.023 0.223 0.142 κ6= 0 0.022 0.192 0.137 P/Pin, at the same locations. It is noteworthy that the minimum levels achieved by the average pressure are larger for the case κ6= 0. This finding is in line with the previous results of Gonzalez and Emanuel [19] concerning the high sensitivity of the pressure field to the Stokes hypothesis and the strong modification of the pressure distribution obtained for large values of the ratio κ/µ. 0 5 10 15 20 25 −0.4 −0.2 0.0 0.2 0.4 0.6 0.8 x2/δω,0 (eu−U2)/∆U x1/δω,0= 140 x1/δω,0= 170 x1/δω,0= 200 (a) 0 5 10 15 20 25 0.0 0.2 0.4 0.6 0.8 1.0 x2/δω,0 e Z x1/δω,0= 140 x1/δω,0= 170 x1/δω,0= 200 (b) Figure 17: Comparison between the normalized velocity and passive scalar profiles obtained for both case κ6= 0 and κ= 0 plotted at three abscissae. The solid lines correspond to κ6= 0 and the lines with symbol (++) to κ= 0. 6 8 10 12 14 1.34 1.35 1.36 1.37 x2/δω,0 e P/Pin x1/δω,0= 140 (a) 8 10 12 14 16 1.84 1.85 1.86 1.87 1.88 x2/δω,0 e P/Pin x1/δω,0= 170 (b) 6 8 10 12 14 16 18 1.86 1.87 1.88 1.89 x2/δω,0 e P/Pin x1/δω,0= 200 (c) Figure 18: Longitudinal profiles of the mean pressure normalized by its value at the inlet Pin, same symbols as those used in Figure 17. Since it has been shown that the first-order statistical moments do not display a significant sensitivity to the bulk viscosity, a closer look is now taken at some second-order moments that characterize the velocity and passive scalar fluctuations. Figure 19 reports the longitudinal evolution of the maximum values of the normalized TKE. It is noteworthy that the threedimensional character of the present set simulation slightly modifies the conclusion that were 16 previously drawn from the two-dimensional computations: the influence of bulk viscosity is noticeable. Up to the abscissa x1/δω,0≈150.0, neglecting the bulk viscosity only leads to a very slight overestimate compared to the case where the effects of κare considered. The region that extends from x1/δω,0= 150.0until the interaction with the reflected shock is characterized by a significant change of behavior and the values obtained with κ6= 0 are larger than those obtained with κ= 0. This region is characterized by strong pressure wave reflection from the upper limit of the computational domain. After the interaction with the reflected shock wave (and up to the abscissa x1/δω,0= 275.0), the maxima of the TKE obtained without taking into account the bulk viscosity effects are again underestimated compared to those issued from the computations performed with κ6= 0 and this trend is slightly modified further downstream as the end of the computational domain is approached. It can be concluded that, in the absence of the second shock wave interaction and associated parasitic pressure waves issued from the top of the computational domain, the TKE levels would be underestimated if the effects of κ were not taken into account. In an attempt to better understand the behavior of the TKE, the analysis of the main terms involved in its transport equation is now carried out. The transport equation for the turbulent kinetic energy Kis given by ∂t(ρK) + ∇·(ρe uK) = P+ε+T+ Π + Σ,(12) In this equation, Pis the production term, εis the dissipation term, Tdenotes the turbulent transport term, Πis the pressure-strain term, and finally Σthe mass flux term. The budget (12) is deduced from the transport equation of the Reynolds tensor components: ∂(ρRij) ∂t +∂(ρeukRij) ∂xk =Pij +εij +Tij + Πij + Σij,(13) with                                        Pij =−ρRik ∂euj ∂xk +Rjk ∂eui ∂xk, εij =−τ0 ik ∂u00 j ∂xk −τ0 jk ∂u00 i ∂xk , Tij =−∂ ∂xkρu00 iu00 ju00 k+P0u00 iδjk +P0u00 jδik −τ0 jku00 i−τ0 iku00 j, Πij =P0∂u00 i ∂xj +P0∂u00 j ∂xi , Σij =u00 i ∂τjk ∂xk +u00 j ∂τik ∂xk−u00 i ∂P ∂xj +u00 j ∂P ∂xi, (14a) (14b) (14c) (14d) (14e) The analysis of the main terms involved in the TKE transport equation is carried out at two distinct locations to infer the impact of the volume viscosity. Figure 20 shows that the most important contributions are associated to the production and dissipation terms. Their amplitude is found to be slightly smaller when κis not taken into account. The turbulent transport term is positive at the periphery of the mixing layer while it tends to be negative within the mixing layer. This quantity, which is larger in the case featuring κ6= 0, removes energy from regions characterized by large fluctuations levels to deposit it in regions characterized by lower levels of TKE. Figure 20 also shows that the contributions due to pressure-strain and mass flux terms remain negligible compared to the others, for both cases. 17 0 50 100 150 200 250 0 2 4 6·10−2 x1/δω,0 max ( ρK)/ρ0∆U2 κ6= 0 κ= 0 Figure 19: Spatial evolution of the longitudinal maxima of normalized turbulent kinetic energy ρK. 4 6 8 10 12 14 0.0 5.0 ·10−4 x2/δω,0 x1/δω,0= 130 P(14a) ε(14b) (a) 4 6 8 10 12 14 16 18 −4.0 −2.0 0.0 2.0 4.0 6.0 ·10−4 x2/δω,0 x1/δω,0= 200 P(14a) ε(14b) (b) 4 6 8 10 12 14 −4.0 −2.0 0.0 2.0 ·10−4 x2/δω,0 x1/δω,0= 130 T(14c) Π(14d) Σ(14e) (c) 4 6 8 10 12 14 16 18 −2.0 −1.0 0.0 1.0 2.0 ·10−4 x2/δω,0 x1/δω,0= 200 T(14c) Π(14d) Σ(14e) (d) Figure 20: TKE budget with all terms normalized by ∆U3/δω,0, same symbols as those retained in Figure 17. Figure 21 reports the distribution of the Reynolds stress components as well as the variance of the passive scalar and the scalar to velocity correlations for both cases κ6= 0 and κ= 0. The three streamwise positions under consideration are representative of the variations observed on the TKE profile reported in Figure 19. The profiles of the Reynolds stress components show that the maxima of its diagonal components follow the trends reported in Figure 19.Figure 21(f), which displays the longitudinal evolution of the scalar flux component ] u00 1ξ00/(u1,RMSξRMS), reveals that the maximum value of the correlation between the longitudinal velocity fluctuation and the scalar fluctuation is slightly underestimated when the effects of bulk viscosity are not considered. Figure 22 reports the variance of the mass fractions of chemical species present in the 18 0 5 10 15 20 25 0.00 0.05 0.10 0.15 0.20 x2/δω,0 q] u00 1u00 1/∆U x1/δω,0= 140 x1/δω,0= 170 x1/δω,0= 200 (a) 0 5 10 15 20 25 0.00 0.05 0.10 0.15 x2/δω,0 q] u00 2u00 2/∆U x1/δω,0= 140 x1/δω,0= 170 x1/δω,0= 200 (b) 0 5 10 15 20 25 0.00 0.05 0.10 0.15 x2/δω,0 q] u00 3u00 3/∆U x1/δω,0= 140 x1/δω,0= 170 x1/δω,0= 200 (c) 0 5 10 15 20 25 0.00 0.05 0.10 x2/δω,0 q] u00 1u00 2/∆U x1/δω,0= 140 x1/δω,0= 170 x1/δω,0= 200 (d) 0 5 10 15 20 25 0.00 2.00 4.00 6.00 ·10−2 x2/δω,0 g ξ00ξ00 x1/δω,0= 140 x1/δω,0= 170 x1/δω,0= 200 (e) 0 5 10 15 20 25 −2.00 −1.00 0.00 ·10−2 x2/δω,0 ] u00 1ξ00/(u1,RMSξRMS) x1/δω,0= 140 x1/δω,0= 170 x1/δω,0= 200 (f) Figure 21: Profiles of the Reynolds stress tensor components normalized by ∆U together with the mixture fraction variance g ξ00ξ00 and longitudinal component of the scalar flux ] u00 1ξ00 at three abscissae, same symbols as those retained in Figure 17. mixture. The hydrogen, which is characterized by the highest ratio κ/µ is the one that displays the largest differences (up to approximately ten percent) between the two cases, i.e., κ= 0 and κ6= 0. The differences observed at the three locations concern both the shape and maximum levels, which depend on the species under consideration. Indeed, it is found that the distribution of the profiles for all chemical species is slightly wider – indicating that the fluid is incorporated more markedly – when the effects of bulk viscosity are taken into account, which leads to a reduction of fluctuations around the averaged value. A similar effect is observed when the convective Mach number values are increased [22,24]. For the present three-dimensional simulation, it is interesting to consider the evolution of 19 6 8 10 12 14 16 0.0 0.2 0.4 0.6 0.8 1.0 1.2·10−4 x2/δω,0 ^ Y00 H2Y00 H2 x1/δω,0= 140 x1/δω,0= 170 x1/δω,0= 200 (a) 6 8 10 12 14 16 0.0 2.0 4.0 6.0·10−3 x2/δω,0 ^ Y00 O2Y00 O2 x1/δω,0= 140 x1/δω,0= 170 x1/δω,0= 200 (b) 6 8 10 12 14 16 0.0 0.5 1.0 1.5 2.0 ·10−3 x2/δω,0 ^ Y00 H2OY00 H2O x1/δω,0= 140 x1/δω,0= 170 x1/δω,0= 200 (c) Figure 22: Profiles of the variances of chemical species mass fractions at three abscissae, same symbols as those retained in Figure 17. higher moments (thirdand fourth-order) of velocity and passive scalar fluctuations. Once properly normalized, these quantities provide the skewness and kurtosis of the statistical distribution, i.e., the probability density functions (PDF). They are defined by S=µ3/σ3and F=µ4/σ4, respectively, with µ3and µ4the third and fourth centered moments, σbeing the standard deviation. The skewness factor measures the symmetry of the fluctuations around the mean whilst the flatness factor characterizes if the PDF tends to be peaked or not. Their values can be compared to those associated to a Gaussian distribution with SG= 0 and FG= 3. Figure 23 displays the transverse profiles of Su1,Su2,Fu1, and Fu2obtained in three planes for κ= 0 and κ6= 0. In the free streams (i.e., outside the mixing layer), the skewness and flatness coefficients tend to 0.0and 3.0, respectively, which reflects a quasi-Gaussian behaviour of the residual turbulence outside the mixing zone. The consideration of bulk viscosity effects favors this trend. A sharp variation of these two coefficients is observed when approaching the mixing zone boundaries. These changes are associated to the intermittency between turbulent “puffs” (mixed fluid) in the freestream as well as fluid incursions into the mixing layer (fluid entrainment). The sign of the variation of the asymmetry coefficient depends upon the stream from which the mixing layer boundary is approached, either from the fast-cold stream side or from the slow-hot stream. On the one hand, on the high-speed and low-temperature side, intermittent events, which disrupt the fast and cold uniform flow, correspond to hotter and slower conditions and give rise to Su1<0and Su2<0. On the other hand, on the low speed and high temperature side, such events correspond to colder and faster fluid particles, and we 20 have Sξ<0,Su1>0, and Su2>0. These intermittent events influence the flatness coefficients as abruptly as the skewness coefficients, see Figure 23. Therefore, the different profiles approach an antisymmetric shape for the skewness coefficients and a symmetrical shape for the flatness coefficients. The positions where the asymmetry coefficients Su1and Su2cancel out correspond to the positions of the extrema on their respective variance profile. The bulk viscosity has the effect of amplifying the amplitude of the peaks observed for the asymmetry and flatness coefficients, which may be explained by the stabilizing action of molecular processes that induces a delayed development of the mixing layer and a subsequent delayed action of molecular processes. In addition, from the inner part of the mixing layer towards its boundaries (inlet stream conditions), it has to be noticed that the intermittent zone appears early in the presence of bulk viscosity effects. 0 5 10 15 20 25 −2 −1 0 1 x2/δω,0 Su1 x1/δω,0= 140 x1/δω,0= 170 x1/δω,0= 200 (a) 0 5 10 15 20 25 −2 −1 0 1 2 x2/δω,0 Su2 x1/δω,0= 140 x1/δω,0= 170 x1/δω,0= 200 (b) 0 5 10 15 20 25 2 4 6 8 10 12 14 x2/δω,0 Fu1 x1/δω,0= 140 x1/δω,0= 170 x1/δω,0= 200 (c) 0 5 10 15 20 25 5 10 15 20 x2/δω,0 Fu2 x1/δω,0= 140 x1/δω,0= 170 x1/δω,0= 200 (d) Figure 23: Skewness and flatness coefficients issued from the statistics of the two velocity components u1and u2, same symbols as those retained in Figure 17. Figure 24 reports the longitudinal evolution of the enstrophy maximum Ωnormalized by (∆U/δω,0)3. These two profiles display three distinct regions. Before the interaction of the second reflected shock with the mixing layer, the maximum value of the enstrophy obtained with κ= 0 seems to be underestimated in comparison to that obtained with κ6= 0 while beyond the second interaction, the two values become quite comparable. Two specific locations that are typical of each region are now considered to study the origin of observable differences in the light of Eq. (8). It has to be noted that the most important contributions to the enstrophy budget are associated to vortex stretching and viscous diffusion. There is an indirect effect of bulk viscosity that tends to promote the stretching term and, as a direct outcome of the increased molecular effects, the viscous diffusion term is also slightly 21 0 50 100 150 200 250 0 0.2 0.4 0.6 0.8 x1/δω,0 max ( Ω ) /( ∆U/δω,0)2 κ6= 0 κ= 0 Figure 24: Longitudinal evolution of the normalized enstrophy maxima. larger when κ6= 0. Finally, the baroclinic term is slightly larger in the absence of bulk viscosity, while the dilatation contribution does not seem to be significantly modified by the consideration of bulk viscosity effects. 4 6 8 10 12 14 −0.05 0.00 0.05 0.10 0.15 x2/δω,0 x1/δω,0= 130 E(8) V(8) (a) 4 6 8 10 12 14 16 18 0.00 0.05 0.10 x2/δω,0 x1/δω,0= 200 E(8) V(8) (b) 4 6 8 10 12 14 −1.50 −1.00 −0.50 0.00 0.50 ·10−2 x2/δω,0 x1/δω,0= 130 B(8) D(8) (c) 4 6 8 10 12 14 16 18 −6.00 −4.00 −2.00 0.00 2.00 ·10−3 x2/δω,0 x1/δω,0= 200 B(8) D(8) (d) Figure 25: Enstrophy budget normalized by (∆U/δω,0)3, same symbols as those retained in Figure 17. 5 Summary and conclusions In the present manuscript, twoand three-dimensional numerical simulations of spatiallydeveloping compressible mixing layers impacted by an oblique shock wave are conducted for a convective Mach number Mc= 0.48. The emphasis is placed on the possible influence of the bulk viscosity on the mixing processes. Thus, a mixture of hydrogen and air is considered 22 in conditions that are representative of experimental benchmarks relevant to high-speed flow combustion. In a first step of the analysis, two-dimensional computations of inert and reactive mixing layers are computed. A significant impact of the bulk viscosity is observed on the instantaneous flowfields while averaged quantities do not exhibit any remarkable modification. It is also worth noting that the reactive cases only display slight differences with respect to inert cases: this is especially true for the longitudinal evolutions of the vorticity thickness and turbulent kinetic energy. Three-dimensional simulations of inert mixing layers are subsequently conducted. The influence of the bulk viscosity is more visible in these three-dimensional cases: it tends to reduce the mixing layer growth rate compared to the case where it is not taken into account. The comparison is also performed in terms of higher-order statistical moments. This last part of the analysis shows that the bulk viscosity effects tend to amplify the velocity gradients at the boundaries of the mixing layer, and consequently favor the return to equilibrium. From the above synthesis of the obtained results, one may expect that refined large-eddy simulations (LES) may be rather sensitive to the consideration of bulk viscosity, while Reynolds-averaged Navier-Stokes (RANS) simulations, which are based on statistical averages, are not. Finally, from the present set of results, it is recommended to take the bulk viscosity effects into account especially when highly-resolved large-eddy simulations (LES) are considered. Acknowledgments The computations were performed using the High Performance Computing resources from the mésocentre de calcul poitevin and from genci under allocations x20142a0912 and x20142b7251. The first author also benefited from interesting discussions with Aimad Er-raiy. References [1] F. Bahmani and M. S. Cramer. Suppression of shock-induced separation in fluids having large bulk viscosities. Journal of Fluid Mechanics, 756, 2014. [2] U. Balucani and M. Zoppi. Dynamics of the liquid state, volume 10. Clarendon Press, 1995. [3] G. Billet, V. Giovangigli, and G. De Gassowski. Impact of volume viscosity on a shock– hydrogen-bubble interaction. Combustion Theory and Modelling, 12(2):221–248, 2008. [4] R. Boukharfane, Z. Bouali, and A. Mura. Evolution of scalar and velocity dynamics in planar shock-turbulence interaction. Shock Waves, 2018 (to appear). [5] M. Brouillette. The richtmyer-meshkov instability. Annual Review of Fluid Mechanics, 34(1):445–468, 2002. [6] D. Bruno and V. Giovangigli. Relaxation of internal temperature and volume viscosity. Physics of Fluids, 23:093104, 2011. [7] R. Buttay, G. Lehnasch, and A. Mura. Analysis of small-scale scalar mixing processes in highly under-expanded jets. Shock Waves, 26(2):93–212, 2016. 23 [8] S. Chapman and Cowling T.G. The mathematical theory of non-uniform gases. Cambridge University Press, 1970. [9] A. V. Chikitkin, B. V. Rogov, G. A. Tirsky, and S. V. Utyuzhnikov. Effect of bulk viscosity in supersonic flow past spacecraft. Applied Numerical Mathematics, 93:47–60, 2015. [10] M. S. Cramer. Numerical estimates for the bulk viscosity of ideal gases. Physics of Fluids, 24:066102, 2012. [11] M. S. Cramer and F. Bahmani. Effect of large bulk viscosity on large-reynolds-number flows. Journal of Fluid Mechanics, 751:142–163, 2014. [12] A. Ern and V. Giovangigli. Multicomponent transport algorithms, volume 24. Springer Science & Business Media, 1994. [13] A. Ern and V. Giovangigli. Fast and accurate multicomponent transport property evaluation. Journal of Computational Physics, 120:105–116, 1995. [14] A. Ern and V. Giovangigli. Volume viscosity of dilute polyatomic gas mixtures. European Journal of Mechanics. B, Fluids, 14(5):653–669, 1995. [15] A. Ern and V. Giovangigli. EGlib User’s Manual, 1996. [16] A. Ern and V. Giovangigli. EGLIB: A general-purpose fortran library for multicomponent transport property evaluation. Technical report, Manual of EGLIB Version, 2004. [17] G. Fru, G. Janiga, and D. Thévenin. Impact of volume viscosity on the structure of turbulent premixed flames in the thin reaction zone regime. Flow, Turbulence and Combustion, 88(4):451–478, 2012. [18] L. Gomet, V. Robin, and A. Mura. A multiple-inlet mixture fraction model for nonpremixed combustion. Combustion and Flame, 162:668–687, 2015. [19] H. Gonzalez and G. Emanuel. Effect of bulk viscosity on couette flow. Physics of Fluids A: Fluid Dynamics, 5(5):1267–1268, 1993. [20] A. C. Hindmarsh, P. N. Brown, K. E. Grant, S. L. Lee, R. Serban, D. E. Shumaker, and C. S. Woodward. Sundials: Suite of nonlinear and differential/algebraic equation solvers. ACM Transactions on Mathematical Software (TOMS), 31(3):363–396, 2005. [21] J. Lin, C. Scalo, and L. Hesselink. High-fidelity simulation of a standing-wave thermoacoustic–piezoelectric engine. Journal of Fluid Mechanics, 808:19–60, 2016. [22] I. Mahle. Direct and large-eddy simulation of inert and reacting compressible turbulent shear layers. PhD thesis, Technische Universität München, 2007. [23] P. J. Martínez Ferrer, R. Buttay, G. Lehnasch, and A. Mura. A detailed verification procedure for compressible reactive multicomponent navier–stokes solvers. Computers & Fluids, 89:88–110, 2014. 24 [24] P. J. Martínez Ferrer, G. Lehnasch, and A. Mura. Compressibility and heat release effects in high-speed reactive mixing layers. growth rates and turbulence characteristics. Combustion and Flame, 180:284–303, 2017. [25] R. E. Mitchell and R. J. Kee. General-purpose computer code for predicting chemicalkinetic behavior behind incident and reflected shocks. Technical report, Sandia National Labs., Livermore, CA (USA), 1982. [26] L. Monchick and E. A. Mason. Transport properties of polar gases. The Journal of Chemical Physics, 35(5):1676–1697, 1961. [27] M. Ó Conaire, H. J. Curran, J. M. Simmie, W. J. Pitz, and C. K. Westbrook. A comprehensive modeling study of hydrogen oxidation. International Journal of Chemical Kinetics, 36(11):603–622, 2004. [28] C. D. Pierce. Progress-variable approach for large-eddy simulation of turbulent combustion. PhD thesis, Stanford University, 2001. [29] G. J. Prangsma, A. H. Alberga, and J. J. M. Beenakker. Ultrasonic determination of the volume viscosity of n2,co,ch4, and cd2between 77 and 300 k. Physica, 64(2):278–288, 1973. [30] J. D. Ramshaw. Simple model for mixing at accelerated fluid interfaces with shear and compression. Physical Review E, 61(5):5339, 2000. [31] Y. Yan, C. Chen, P. Lu, and C. Liu. Study on shock wave-vortex ring interaction by the micro vortex generator controlled ramp flow with turbulent inflow. Aerospace Science and Technology, 30(1):226–231, 2013. [32] J. L. Ziegler, R. Deiterding, J. E. Shepherd, and D. I. Pullin. An adaptive high-order hybrid scheme for compressible, viscous flows with detailed chemistry. Journal of Computational Physics, 230(20):7598–7630, 2011. 25