Aurélien Vadrot

Alexis Giauque

Christophe Corre

Direct Numerical Simulations of temporal compressible mixing layers in a BZT Dense Gas: influence of the convective Mach number

). Several 3D direct numerical simulation (DNS) of compressible mixing layers are performed with FC-70 using the fifth order Martin-Hou thermodynamic equation of state (EoS) and air using the perfect gas (PG) EoS. After having carefully defined self-similar periods using the temporal evolution of the integrated streamwise production term, the evolutions of the mixing layer growth rate as a function of the convective Mach number are compared between perfect gas and dense gas flows. Results show major differences for the momentum thickness growth rate at M c = 2.2. The well-known compressibility-related decrease of the momentum thickness growth rate is reduced in the dense gas. Fluctuating thermodynamics quantities are strongly modified. In particular, temperature variations are suppressed leading to an almost isothermal evolution. The small scales dynamics is also influenced by dense gas effects, which calls for a specific sub-grid scale modelling when computing dense gas flows using large eddy simulation (LES). Additional dense gas DNS are performed at three others initial thermodynamic operating points. DNS performed outside and inside the BZT inversion region do not show major differences. BZT effects themselves therefore only have a small impact on the mixing layer growth.

Introduction

Dense gases (DG) are single-phase vapours characterized by long chains of atoms and by medium to large molecular weights. They have been widely used in the Organic Rankine Cycles (ORCs) industry over the past forty years. Their large heat capacity and their low boiling point temperature make them suitable working fluids for lowtemperature heat sources (solar, geothermal, biomass, heat recovery). The coupling with a turbine enables power generation. Recently, because of issues caused by carbon-based fossil fuels, there has been a strong research effort in developing this technology by improving ORC turbine efficiency. Rotating elements are a main source of losses for turbines. Their use in transonic and supersonic regimes are associated with shocks which generate entropy. However, for dense gases, entropy jumps through shocks are significantly reduced in specific thermodynamic regions [START_REF] Cinnella | Inviscid and viscous aerodynamics of dense gases[END_REF]. This feature could enable to increase ORC turbines efficiency, but the lack of knowledge about dense gases in these particular thermodynamic regions close to the vicinity of the critical point restrains ORC designers. This study seeks to widen knowledge about turbulence characteristics of these gases by comparing their behaviour to perfect gases on a classical configuration: the mixing layer.

A specific type of dense gas is used in these simulations: the Bethe Zel'dovich Thompson (BZT) gases, whose name was given at first by [START_REF] Cramer | Nonclassical dynamics of classical gases[END_REF] to acknowledge pioneering works of [START_REF] Bethe | The theory of shock waves for an arbitrary equation of state[END_REF], [START_REF] Zel'dovich | On the possibility of rarefaction shock waves[END_REF] and [START_REF] Thompson | A fundamental derivative in gasdynamics[END_REF]. Unlike other dense gases, they comprise an inversion thermodynamic region where the fundamental derivative of gas dynamics Γ becomes negative as shown in Fig. 1. [START_REF] Thompson | A fundamental derivative in gasdynamics[END_REF] defines Γ as:

Γ = v 3 2c 2 ∂ 2 p ∂v 2 s = c 4 2v 3 ∂ 2 v ∂p 2 s = 1 + ρ c ∂c ∂ρ s (1.1)
where v is the specific volume, ρ the density, c = ∂p/∂ρ| s the speed of sound, p the pressure and s the entropy. For thermally and calorically perfect gases, the fundamental derivative is equal to (γ + 1)/2, with γ the heat capacity ratio. In this case, its value is always greater than one, unlike dense gas flows, where Γ can become lower than one and even be negative for BZT dense gases. In that case, rarefaction shock-waves can occur, which is forbidden by the Second Law of Thermodynamics in usual gases, where only compression shock-waves are allowed. [START_REF] Bethe | The theory of shock waves for an arbitrary equation of state[END_REF] expressed the entropy jump expression across shock-waves as a function of the fundamental derivative.

∆s = s 2 -s 1 = - ∂ 2 p ∂v 2 s ∆v 3 12T + O ∆v 4 = - c 2 Γ v 3 ∆v 3 6T + O ∆v 4 (1.2)
with T the temperature. In the case of compression shock-waves, the specific volume variation is negative (∆v < 0), so that the fundamental derivative must be positive (Γ > 0) to ensure that the entropy jump remains positive (∆s > 0), thus satisfying the Second law of thermodynamics. Only compression shock-waves are physically admissible for classical ideal gases since Γ > 1. For BZT gases, the fundamental derivative being negative (Γ < 0), physically admissible shock-waves in the inversion region are expansion shock-waves such that the specific volume variation is positive (∆v > 0) to ensure the entropy jump remains positive. Moreover, since entropy jumps are proportional to the fundamental derivative Γ , which is of small amplitude in DG flows, intensity of shocks is significantly reduced [START_REF] Cramer | On the propagation of waves exhibiting both positive and negative nonlinearity[END_REF]. In addition to a peculiar thermodynamic behaviour, the sound speed is much lower in dense gases when compared to perfect gases, which makes compressibility regimes much more easily accessible.

Up to now, although thermodynamic features of dense gases are very different from ones of perfect gases, in absence of a better option, perfect gas turbulence closure models coupled with real-gas thermal and calorific EoS have been used for Reynolds-Averaged Navier Stokes (RANS) and Large Eddy Simulation (LES) to simulate dense gas flows [START_REF] Cinnella | Numerical solver for dense gas flows[END_REF][START_REF] Wheeler | A study of the three-dimensional unsteady real-gas flows within a transonic ORC turbine[END_REF][START_REF] Galiana | A study of trailing-edge losses in organic rankine cycle turbines[END_REF]. This choice implicitly assumes that turbulent structures are not affected by dense gas effects. This hypothesis is not yet verified and constitutes an open research field. There is currently no experimental data to verify this hypothesis because maintaining the flow in the vicinity of the critical point where physical quantities are experiencing strong variations is a very complex task.

Direct Numerical Simulation (DNS) is the tool of choice used in this study to assess this hypothesis. DNS enables to solve every turbulent scale down to the smallest one corresponding to the Kolmogorov length scale without resorting to any turbulence closure model. So far, few DNS of dense gas flows have been achieved. DNS of decaying Homogeneous Isotropic Turbulence (HIT) performed by [START_REF] Giauque | Direct numerical simulations of homogeneous isotropic turbulence in a dense gas[END_REF] shows that the dynamic Smagorinsky sub-grid scale model is not able to correctly capture the temporal decay of the turbulent kinetic energy. They extended their analysis by performing a forced HIT highlighting significant differences in the sub-grid scale (SGS) baropycnal work and the resolved pressure-dilatation, which is reduced by a factor of 2 in the dense gas (DG) when compared to the perfect gas (PG) [START_REF] Giauque | Direct numerical simulations of forced homogeneous isotropic turbulence in a dense gas[END_REF]. Sciacovelli et al. (2017b) performed DNS of decaying HIT and notice reduced levels of thermodynamic fluctuations in dense gas flows due to the decoupling of thermal and dynamic phenomena caused by the large heat capacity. The Eckert number, which quantifies the ratio between the kinetic energy and the internal energy, is indeed much smaller in dense gas flows. They also display a more symmetric probability density function (PDF) of the velocity divergence in BZT DG flows, explained by the presence of expansion shocklets and by the attenuation of compression shocklets. They show that turbulence structures are modified by expansion regions: the occurrence of non-focal convergent structures in DG flows diminishes the vorticity and counterbalances enstrophy destruction. Sciacovelli et al. (2017a) analyse DG flow behaviour in a turbulent channel flow. The initial thermodynamic state was this time chosen in a non-BZT DG region. They observe significant differences with respect to PG flows in thermodynamic variables.

A. Vadrot, A. Giauque and C. Corre Temperature variations are small in DG which leads to an almost isothermal evolution. The viscosity decreases from the wall towards the centreline unlike in PG flows. They also notice significant differences in the shape and rates of the fluctuating density and temperature distributions. It is also found that the structure of turbulence is not deeply affected in DG flows. An extent of this study to the BZT DG region and to a larger Mach number would help to conclude on BZT DG effect on turbulence development. [START_REF] Gloerfelt | Dense-gas effects on compressible boundary-layer stability[END_REF] performed the DNS of a dense gas compressible boundary layer at Mach numbers ranging from 0.5 to 6. They especially confirm the decoupling between dynamical and thermal effects, which leads to a suppression of friction heating. The most remarkable consequence is that the boundary layer thickness remains equal to its value in the incompressible regime as the Mach number increases.

Recently, [START_REF] Vadrot | Analysis of turbulence characteristics in a temporal dense gas compressible mixing layer using direct numerical simulation[END_REF] performed DNS of temporal compressible mixing layers for BZT DG flow and PG flow at a convective Mach number M c = 1.1, which is defined as:

M c = (u 1 -u 2 )/(c 1 + c 2 ) (1.3)
where u i and c i denote the flow speed and the sound speed of stream i (upper or lower) of the mixing layer. They show that the mixing layer is significantly affected by dense gas effects during the initial unstable growth phase, revealing a much faster unstable growth in the DG flow. However, only slight differences are observed during the self-similar period, which is the regime of interest when studying mixing layers. Self-similarity is thoroughly described in section 3.1. Results from this initial study at M c = 1.1 also show that the turbulent Mach number (Equation 1.4) is in the low-limit to get shocklets.

M t = u i u i c (1.4)
The authors expect that shocklets, which exhibit very different properties in DG flows when compared to PG flows, would have an impact on the mixing layer growth. In order to account for these additional effects, an extent of the study to larger convective Mach numbers is hereby considered. Since it is known that there are major differences between BZT DG flow and PG flow in shocklet generation, a study in a higher compressible regime would help to answer the following question: is the mixing layer growth rate modified in BZT dense gas flows ?

Since the past thirty years, many DNS of mixing layers have been achieved. The first ones were performed by [START_REF] Sandham | Compressible mixing layer -Linear theory and direct simulation[END_REF]; [START_REF] Luo | On the formation of small scales in a compressible mixing layer[END_REF]; [START_REF] Vreman | Compressible mixing layer growth rate and turbulence characteristics[END_REF]. These DNS use the perfect gas hypothesis. A common feature of compressible mixing layers, shown by experiments at first and DNS afterwards, is the reduction of the mixing layer growth rate with the increase of the convective Mach number. However, detailed mechanisms responsible for this trend are still under investigations.

At first, additional terms in the turbulent kinetic energy equation due to compressibility effects: compressible dissipation d and pressure-dilatation Π ii were suspected to be responsible for the growth rate reduction. [START_REF] Zeman | Dilatation dissipation: the concept and application in modeling compressible mixing layers[END_REF] and Sarkar et al. (1991) especially proposed models for the dilatation dissipation. However, it was shown by [START_REF] Sarkar | The stabilizing effect of compressibility in turbulent shear flow[END_REF] that the growth rate diminution is primarily due to the reduction of turbulent production and not to dilatation terms. [START_REF] Vreman | Compressible mixing layer growth rate and turbulence characteristics[END_REF] confirmed that dilatation terms play a minor role in mixing layer growth and extended previous analysis, showing that pressure-strain terms Π ij diminution is responsible for the turbulent production decrease. They also noticed thanks to DNS that this diminution is mainly due to the decrease of pressure fluctuations normalised by the dynamic pressure (p rms /( 1 2 ρ 0 (∆u) 2 )). [START_REF] Pantano | A study of compressibility effects in the high-speed turbulent shear layer using direct simulation[END_REF] later demonstrated analytically the aforementioned observation. [START_REF] Hamba | Effects of pressure fluctuations on turbulence growth in compressible homogeneous shear flow[END_REF] performed the DNS of an homogeneous shear flow varying M t from 0.1 to 0.3. The author identifies a dissipative term, responsible for the normalised pressure fluctuations diminution, in the transport equation for p 2 called pressure-variance dissipation and which depends on the thermal conductivity. Several turbulence models were next proposed, based on the normalised pressure fluctuations reduction [START_REF] Fujiwara | A turbulence model for the pressure-strain correlation term accounting for compressibility effects[END_REF][START_REF] Park | A compressible turbulence model for the pressure-strain correlation[END_REF][START_REF] Huang | Modelling of pressure-strain correlation in compressible turbulent flow[END_REF].

However, few experiments and DNS have been achieved at high M c . [START_REF] Rossmann | Evolution and growth of large scale structures in high compressibility mixing layers[END_REF] have experimentally studied higher compressibility regimes up to M c = 2.25 and [START_REF] Matsuno | Compressibility effects in high speed turbulent shear layers-revisited[END_REF] recently performed DNS of temporal mixing layers up to M c = 2.0, but none of them is performed for real gas let alone for DG flows.

In the present article, several 3D DNS of compressible dense gas mixing layers are performed for the first time at M c = 2.2. A comparison is made between PG and DG flows. Evolution of the mixing layer growth rate as a function of the convective Mach number is compared between perfect gas and dense gas flows. This study extends previous analysis conducted at M c = 1.1 [START_REF] Vadrot | Analysis of turbulence characteristics in a temporal dense gas compressible mixing layer using direct numerical simulation[END_REF]).

An unusual behaviour is noticed, as the decrease of the mixing layer growth rate with the convective Mach number does not follow the same evolution between DG and PG flows. The discrepancy is not significant at lower Mach number M c = 1.1 [START_REF] Vadrot | Analysis of turbulence characteristics in a temporal dense gas compressible mixing layer using direct numerical simulation[END_REF]) but when the convective Mach number increases, DG mixing layer growth is influenced by modified thermodynamics behaviour. Differences are first analysed in the context of the peculiar shocklets properties in BZT DG flows. Finally, thermodynamics behaviour of DG flows is also investigated.

The first section is devoted to the problem description exposing the main physical and numerical parameters. Results are validated for the perfect gas flow in the second section with a comparison to available results in the literature. Comparison is made between dense gas and perfect gas in section 4. Finally, a physical analysis of discrepancies between DG and PG flows is conducted thanks to additional DNS performed at different thermodynamic operating points (Section 5). The aim of this analysis is to highlight and explain differences between BZT DG and PG flows at large convective Mach number.

Problem formulation

Initialisation

The problem consists in extending the analysis conducted at M c = 1.1 in [START_REF] Vadrot | Analysis of turbulence characteristics in a temporal dense gas compressible mixing layer using direct numerical simulation[END_REF] by performing a DNS of a 3D mixing layer at a convective Mach number M c = 2.2 for air considered as a perfect gas and for a BZT dense gas: the perfluorotripentylamine (FC-70, C 15 F 33 N ). Physical parameters associated to FC-70 and used in these DNS are given in table 1.

The initial thermodynamic state is chosen inside the inversion region in order to favour the occurrence of expansion shocklets, physically allowed in BZT dense gases. Figure 1 shows the initial state in the p -v diagram and its distribution during the beginning of the self-similar regime at τ = 4000 for DG flow. The initial value of the fundamental derivative is Γ initial = -0.284 which makes possible the appearance of expansion shocklets. The distribution spreads inside and slightly outside the inversion region. One can also note that the distribution does not perfectly follow the initial adiabatic curve. Mechanical dissipation and shocklets entropy losses are responsible for this discrepancy because their effect cannot be neglected at M c = 2.2. [START_REF] Cramer | Negative nonlinearity in selected fluorocarbons[END_REF]. The critical pressure pc, the critical temperature Tc, the boiling temperature T b and the compressibility factor Zc = pcvc/(RTc) are the input data for the Martin-Hou equation. The critical specific volume vc is deduced from the aforementioned parameters. The exponent n and the cv(Tc)/R ratio are used to compute the heat capacity cv(T ) (R = R/M being the specific gas constant computed from the universal gas constant R and M , the gas molar mass). For air, the same values of reduced specific volume and reduced pressure are selected for the initial thermodynamic state. Critical values used for air are the critical pressure p c = 3.7663 × 10 6 P a and the specific volume v c = 3.13 × 10 -3 m 3 .kg -1 [START_REF] Stephan | The thermal conductivity of fluid air[END_REF].

Tc (K) pc (atm) Zc T b (K) m (= cv(Tc)/R) n FC-

Mc ρ1/ρ2

Key non-dimensional parameters are the convective Mach number (Equation 1.3) and the Reynolds number based on the initial momentum thickness δ θ,0 :

Re δ θ,0 = ∆uδ θ,0 /ν (2.1)
where ν denotes the kinematic viscosity and the momentum thickness at time t is defined as:

δ θ (t) = 1 ρ 0 (∆u) 2 +∞ -∞ ρ (∆u) 2 4 -ũ2 x dy (2.2)
with ρ 0 = (ρ 1 +ρ 2 )/2 the averaged density and ũx the Favre averaged streamwise velocity defined in Eq. 2.9. The initial momentum thickness Reynolds number is set equal to 160 for all the DNS following [START_REF] Pantano | A study of compressibility effects in the high-speed turbulent shear layer using direct simulation[END_REF]. Table 2 summarizes the computational parameters of simulations performed for different M c (domain size, number of grid elements, dimensional values of velocity, initial momentum thickness and initial turbulent structures sizes). Additional DG simulations given in Appendix A have been performed for other domain sizes and resolutions to validate the current DNS. The impact on the selection of the self-similar period is also analysed in Appendix A.

The temporal mixing layer consists of two streams flowing in opposite directions. The velocity in the upper part of the domain U 1 is set equal to -∆u/2, whereas U 2 is set to ∆u/2. A representation of the computational domain is provided in figure 2. A snapshot of the velocity magnitude is also plotted for the DG DNS at M c = 2.2. A further analysis of the flow field visualisation is given in Appendix B. Periodic boundary conditions are imposed in the x and z directions and non-reflective conditions are set in the y directions using the NSCBC (Navier-Stokes Characteristic Boundary Conditions) model proposed by [START_REF] Poinsot | Boundary conditions for direct simulations of compressible viscous flows[END_REF].

The streamwise velocity field is initialised using an hyperbolic tangent profile:

ūx (y) = ∆u 2 tanh - y 2δ θ,0
(2.3)

The complete streamwise velocity field is obtained by adding fluctuations to the average velocity. For the y and z components, the average velocity is set equal to zero. A Passot-Pouquet spectrum is imposed for initial velocity fluctuations:

E(k) = (k/k 0 ) 4 exp(-2(k/k 0 ) 2 ) (2.4)
where k denotes the wavenumber. The peak wavenumber k 0 controls the size of the initial turbulent structures. Its influence on the mixing layer growth is investigated in Appendix A. Its value only influences the initial unstable growth regime. It has been noted that a larger value of k 0 accelerates the transition to the unstable growth. Its value for each DNS is given in table 2. The velocity field is then filtered to initialize turbulence only inside the initial momentum thickness.

Governing equations

In order to describe the temporally evolving mixing layer, the unsteady, threedimensional, compressible Navier-Stokes equations are solved:

∂ρ ∂t + ∂(ρu i ) ∂x i = 0 (2.5) A. Vadrot, A. Giauque and C. Corre ∂(ρu i ) ∂t + ∂(ρu i u j ) ∂x j = - ∂p ∂x i + ∂τ ij ∂x j (2.6) ∂(ρE) ∂t + ∂[(ρE + p)u j ] ∂x j = ∂(τ ij u i -q j ) ∂x j (2.7)
where

τ ij = µ( ∂ui ∂xj + ∂uj ∂xi -2 3 ∂u k ∂x k δ ij )
denotes the viscous stress tensor (µ the dynamic viscosity), E = e + 1 2 u i u i , the specific total energy (e, the specific internal energy), q j = -λ ∂T ∂xj , the heat flux given by Fourier's law (λ the thermal conductivity). Part of this study is conducted thanks to the analysis of the turbulent kinetic energy (TKE) equation terms. It requires to decompose density, pressure and velocity into mean and fluctuating components as follows:

   ρ = ρ + ρ p = p + p u i = ũi + u i (2.8)
where φ denotes the Reynolds average for a flow variable φ while the Favre average φ is defined as :

φ = ρφ ρ (2.9)
Reynolds fluctuations are noted φ while Favre fluctuations are noted φ . Reynolds averaging is equivalent to plane averaging along x and z directions because of the use of periodic boundary conditions. The TKE equation is obtained from the Navier-Stokes equation by applying the averaging process:

∂ ρk ∂t + ∂ ρk ũj ∂x j = -ρu i u j ∂ ũi ∂x j Production -τ ij ∂u i ∂x j Dissipation - 1 2 ∂ρu i u i u j ∂x j Turbulent transport - ∂p u i ∂x i Pressure transport + ∂u i τ ij ∂x j Viscous transport +p ∂u i ∂x i Pressure dilatation -u i ∂ p ∂x i - ∂ τij ∂x j Mass-flux term (2.10)
where k = 1 2 u i u i denotes the specific turbulent kinetic energy. The main terms of equation 2.10 are production, dissipation and transport terms. Pressure dilatation and mass-flux term (the later comprises the baropycnal work) are equal to zero in the incompressible case. The dissipation term can be decomposed into a solenoidal, a low-Reynolds number and a dilatational component. The latter is associated to losses occurring in eddy shocklets. [START_REF] Lee | Eddy shocklets in decaying compressible turbulence[END_REF] expressed the dilatational dissipation also called the compressible dissipation as:

d = - 4 3 ν ∂u k ∂x k 2 -2u k ∂ν ∂x k ∂u k ∂x k (2.11)
This expression comprises the effect of viscosity variations unlike [START_REF] Sarkar | Application of a reynolds stress turbulence model to the compressible shear layer[END_REF] and [START_REF] Zeman | Dilatation dissipation: the concept and application in modeling compressible mixing layers[END_REF] who expressed it as d = -4 3 ν

∂u k ∂x k 2
, neglecting viscosity variations. For decaying compressible turbulence, [START_REF] Lee | Eddy shocklets in decaying compressible turbulence[END_REF] found that [START_REF] Sarkar | Application of a reynolds stress turbulence model to the compressible shear layer[END_REF] and [START_REF] Zeman | Dilatation dissipation: the concept and application in modeling compressible mixing layers[END_REF]'s expression overestimates by about 15% the compressible dissipation.

In addition to equations 2.5, 2.6 and 2.7, the thermal perfect gas and the following calorific equations of state (EoS) are used for air.

       p = ρRT e = e ref + T T ref c v (T )dT (2.12)
where R is the specific gas constant, c v the specific heat capacity, p the pressure, T the temperature, ρ the density.

For FC-70, the Martin-Hou EoS (referred as MH) will be retained to provide an accurate representation of BZT dense gas thermodynamic behaviour [START_REF] Guardone | Assessment of thermodynamic models for dense gas dynamics[END_REF].

               p = RT v -b + 5 i=2 A i + B i T + C i e -kT /Tc (v -b) i e = e ref + T T ref c v (T )dT + 5 i=2 A i + C i (1 + kT /T c )e -kT /Tc (i -1)(v -b) i-1
(2.13) where (.

) ref denotes a reference state, b = v c (1 -(-31, 883Z c + 20.533)/15), k = 5
.475 and the coefficients A i , B i and C i are numerical constants determined by [START_REF] Martin | Development of an Equation of State for Gases[END_REF] and [START_REF] Martin | An improved equation of state for gases[END_REF] from physical parameters summarized in table 1.

To complete the thermodynamic description of the BZT dense gas, Chung's model is used to compute dynamic viscosity and thermal conductivity [START_REF] Chung | Generalized multiparameter correlation for nonpolar and polar fluid transport properties[END_REF]. FC-70 is assumed to behave as a nonpolar gas, its dipole moment is neglected [START_REF] Shuely | Model liquid selection based on extreme values of liquid state properties in a factor analysis[END_REF]. For PG transport coefficients, the Sutherland's model is used associated to a constant Prandtl number set equal to 0.71. Values of the initial Prandtl number are given in Appendix C for DG flows. The selected constants for Sutherland's law are the ones given by White (1998).

Numerical Setup

DNS are performed using the explicit and unstructured numerical solver AVBP. It solves the 3D unsteady compressible Navier-Stokes equations coupled with the perfect gas EoS (Equation 2.12) for Air and the MH EoS for FC-70 (Equation 2.13) using a two-step time-explicit Taylor Galerkin scheme (TTGC) for the hyperbolic terms based on a cell vertex formulation [START_REF] Colin | Development of high-order taylor-galerkin schemes for LES[END_REF]. The scheme provides high spectral resolution and low numerical dissipation ensuring a third-order accuracy in space and in time. AVBP is designed for massively parallel computation and can be used to perform LES as well as DNS simulations [START_REF] Desoutter | DNS and modeling of the turbulent boundary layer over an evaporating liquid film[END_REF][START_REF] Cadieux | DNS and LES of separated flows at moderate reynolds numbers[END_REF]. The scheme is completed with a shock capturing method. In regions where strong gradients exist, an additional dissipation term is added following the approach of [START_REF] Cook | A high-wavenumber viscosity for highresolution numerical methods[END_REF]. Its impact on the resolution of the smallest scales has been analysed in a previous article [START_REF] Giauque | Direct numerical simulations of forced homogeneous isotropic turbulence in a dense gas[END_REF].

DNS of Perfect gas mixing layer: verification and validation

This section is devoted to the selection of self-similar periods and the assessment of the quality of perfect gas DNS performed for air at three different convective Mach numbers (M c = 0.1/1.1/2.2).

Temporal evolution and self-similarity

Figure 3 shows the temporal evolution of the momentum thickness normalised by its initial value. This key quantity characterizes the development of mixing layers. Time is non-dimensional (τ = t∆u/δ θ,0 ). The evolution is plotted for three different convective Mach numbers (M c = 0.1/1.1/2.2). Results at M c = 1.1 are extracted from [START_REF] Vadrot | Analysis of turbulence characteristics in a temporal dense gas compressible mixing layer using direct numerical simulation[END_REF]. The same Reynolds number (Re δ θ,0 = 160) based on the initial momentum thickness is used for the three different DNS. Simulation parameters are given in table 2. At M c = 2.2, the size of initial turbulent structures has been enlarged in order to speed up the development of the mixing layer.

One can identify three main phases: an initial delay caused by a transition of modes from the modes in which turbulent kinetic energy is initially injected to the most unstable ones; an unstable over-linear growth; and the self-similar period, during which the mixing layer evolves linearly with time. The procedure used to select the self-similar period is detailed in subsequent paragraphs.

At M c = 2.2, one can notice that the mixing layer takes a much longer time to develop. This is consistent with observations of [START_REF] Pantano | A study of compressibility effects in the high-speed turbulent shear layer using direct simulation[END_REF] who noticed that the time necessarily to reach self-similar regime increases with compressibility. Self-similarity is reached around τ ≈ 11500 after a long unstable growth phase. As a comparison, at M c = 0.1 and M c = 1.1, self-similarity is reached respectively at τ = 700 and τ = 1700. Moreover, the self-similar period is also stretched as the convective Mach number increases.

A long time delay is observed at the beginning of the simulation. That delay is associated to the transition of modes. Turbulent kinetic energy is initially injected at a given integral length set equal to L x /8. Afterwards, energy is distributed over the whole spectrum and some unstable modes are amplified leading to the unstable growth phase. In order to reduce this time delay, initial turbulent structures have been chosen larger in proportion to the initial momentum thickness at M c = 2.2 (Table 2). This modification of initial turbulent structures size does not impact the growth rate over the self-similar regime. This has been carefully verified for DG flows in Appendix A.

In addition, domain lengths are doubled in x and z directions and multiplied by four in the y direction when compared to DNS at M c = 1.1 relatively to initial momentum thicknesses. This enables the mixing layer to develop until larger values of δ θ (t)/δ θ,0 and to obtain a long enough self-similar period without reaching domain boundaries. Other simulations performed with smaller domains did not allow the flow to reach self-similarity.

Slopes and standard deviations mentioned in figure 3 are computed over the self-similar period. One can observe that the growth rate is divided by a factor of about two between DNS at M c = 2.2 and at M c = 1.1. Indeed, compressibility effects tend to reduce mixing layer development as the convective Mach number increases.

DNS performed at M c = 0.1 constitutes our reference incompressible case used to plot δθ / δθ,inc = f (M c ). The computed growth rate is about 0.0131 which is relatively close to the empirical value of 0.016 given by [START_REF] Pantano | A study of compressibility effects in the high-speed turbulent shear layer using direct simulation[END_REF]. One can notice a very short unstable growth phase when compared to larger convective Mach number cases.

Self-similarity is a major characteristic of mixing layers: during the self-similar period, flow development can be described using single length and velocity scales. The momentum thickness linearly evolves with time. This particular state in the development of mixing layers is widely used to extract key features of mixing layers. The well known chart giving the evolution of the mixing layer growth rate as a function of the convective Mach number [START_REF] Papamoschou | The compressible turbulent shear layer: an experimental study[END_REF]) is plotted during the self-similar regime. This period is also used to investigate the balance of the TKE equation, because temporal solutions can be averaged during self-similarity since the flow is in a statistically stable state.

The selection of the self-similar period is thus a key point in the study of turbulent mixing layers, but this choice is difficult especially at high compressible regimes which require lengthy simulations. One can note that in our case the time required to achieve self-similarity is multiplied by a factor of about five when the convective Mach number increases from

M c = 1.1 to M c = 2.2.
Lots of authors evoke difficulties in reaching self-similarity [START_REF] Pantano | A study of compressibility effects in the high-speed turbulent shear layer using direct simulation[END_REF][START_REF] Pirozzoli | Early evolution of the compressible mixing layer issued from two turbulent streams[END_REF] particularly because of computational domain lengths. Moreover, criteria to define self-similarity are not standardised. Superposition of the mean velocity profiles, linear evolution of the momentum thickness, collapse of the Reynolds stress profiles are three different ways to define the self-similar period.

The same methodology used in [START_REF] Vadrot | Analysis of turbulence characteristics in a temporal dense gas compressible mixing layer using direct numerical simulation[END_REF] is applied here to select the selfsimilar period: it relies on the stabilisation of the streamwise production term integrated over the whole domain. The underlying reason for using this criterion comes from [START_REF] Vreman | Compressible mixing layer growth rate and turbulence characteristics[END_REF] who demonstrated the following relation between the mixing layer growth rate and the production power (ρP xx = -ρu x u y ∂ ũx ∂y ):

δ θ = dδ θ dt = 2 ρ 0 (∆u) 2 ρP xx dy (3.1)
Figure 4 shows the temporal evolution of the non-dimensional streamwise production integrated over the whole domain for the three DNS at M c ranging from 0.1 to 2.2 performed for air using the PG EoS. A constant integrated production is directly related to a self-similar regime according to Eq. 3.1. Selected self-similar periods are indicated on each plot. As the convective Mach number increases, the maximum peak of integrated turbulent production decreases which is consistent with the decrease of the momentum thickness growth rate. Time required to achieve self-similarity lengthens but self-similar periods last longer.

Difficulties can be encountered to get a fully stable plateau with an almost constant integrated turbulent production. Domain lengths have a major influence on self-similarity. The evolution of the turbulent production follows a piecewise decrease, reaching several plateaus. It is observed that these piecewise plateaus are directly related to integral lengths scales. When some turbulent structures grow and become too large for the computational domain, the integrated turbulent production decreases and reaches another plateau lower than the previous one. The mixing layer therefore adapts its growth to domain lengths when the computational box is not large enough. Since the integrated turbulent production is related to the mixing layer growth rate, a lower plateau leads to a smaller mixing layer growth rate. Great care therefore needs to be taken selecting the size of the computational domain as well as a good stabilization of the integrated turbulent production in order to precisely select the self-similar period. Influence of the domain size on self-similarity is thoroughly investigated in Appendix A for dense gas flows and correlations with integral length scales are analysed.

Validation over the self-similar period

Since self-similar periods are now well defined for each DNS, it is possible to plot the evolution of the mixing layer growth rate with respect to the convective Mach number. Figure 5 shows a comparison between current PG results and available numerical [START_REF] Freund | Compressibility effects in a turbulent annular mixing layer. Part 1. Turbulence and growth rate[END_REF][START_REF] Pantano | A study of compressibility effects in the high-speed turbulent shear layer using direct simulation[END_REF][START_REF] Kourta | Computation of supersonic mixing layers[END_REF][START_REF] Fu | Numerical simulation of compressible mixing layers[END_REF][START_REF] Zhou | Direct numerical simulation of a spatially developing compressible plane mixing layer: flow structures and mean flow properties[END_REF][START_REF] Ferrer | Compressibility and heat release effects in high-speed reactive mixing layers I.: Growth rates and turbulence characteristics[END_REF][START_REF] Matsuno | Compressibility effects in high speed turbulent shear layers-revisited[END_REF] and experimental results [START_REF] Rossmann | Evolution and growth of large scale structures in high compressibility mixing layers[END_REF]) from the literature. Current DNS follow the tendency observed and described in the literature: the well-known compressibility-related reduction of the momentum thickness growth rate as M c increases. From the incompressible case to M c = 2.2, the mixing layer growth rate is divided by a factor of about five. Standard deviations have also been computed and are reported on the plot. It represents about 5% of the computed growth rates. It is rather difficult to reduce this uncertainty because of difficulties encountered in reaching perfect self-similarity. This is also illustrated by the scattering of literature results, which might be a consequence of this phenomenon. Moreover, the lack of numerical results at highly compressible regimes makes the validation process more complex.

Yet, numerical parameters given in table 3 confirm the validation of the current DNS. The integral lengths l x and l z are computed using the streamwise velocity field:

l x = 1 2u x 2 Lx/2 -Lx/2 u x (x)u x (x + re x )dr (3.2) l z = 1 2u x 2 Lz/2 -Lz/2 u x (x)u x (x + re z )dr (3.3)
Integral length scales show that the domain is chosen sufficiently large. The largest value 0.20 is obtained at the end of the self-similar period for DG flow at M c = 1.1. Otherwise, values do not exceed 0.16 in the streamwise direction and 0.13 in the z direction. As a comparison, [START_REF] Pantano | A study of compressibility effects in the high-speed turbulent shear layer using direct simulation[END_REF] 0.178 in the streamwise direction for a configuration with M c = 0.7 and a density ratio of 4. Appendix A also confirms that domain lengths have been properly chosen for DG mixing layer at M c = 2.2.

The ratio r = L η /∆x characterizes the resolution of simulations. The larger the ratio, the better the resolution. Minimum value is about 0.52 computed for DNS at M c = 2.2. For other simulations, values are larger than 0.6 and the maximum value is 1.64 for PG at M c = 2.2 because of small dissipation in high compressible regimes. As a comparison [START_REF] Pantano | A study of compressibility effects in the high-speed turbulent shear layer using direct simulation[END_REF]'s ratio is about 0.38 for the most resolved simulation and recently [START_REF] Matsuno | Compressibility effects in high speed turbulent shear layers-revisited[END_REF] performed a DNS at M c = 2.0 with a L η /dx ratio equal to 0.41. One can thus consider that turbulent scales are adequately resolved for all simulations presented in this paper since in addition the turbulent kinetic energy is very low close to the Kolmogorov scale [START_REF] Moin | Direct numerical simulation: a tool in turbulence research[END_REF].

Dense gas effect on mixing layer growth

Temporal evolution

As previously done for the perfect gas mixing layer, it is required to precisely define the self-similar range for the dense gas flow. This is done through both figures 6 and 7. Figure 6 enables the comparison of normalised DG momentum thickness over time at three different convective Mach numbers : M c = 0.1 -1.1 -2.2. The three DNS are performed at the same initial Reynolds number Re δ θ,0 = 160. Additional simulation parameters are given in table 2. At M c = 0.1, similarly to PG mixing layer, the domain length is doubled in the y direction to get a long enough self-similar period. At M c = 2.2, the domain length is divided by two in the y direction when compared to PG flow. The domain is therefore large enough to reach a self-similar period which lasts 4000τ . Initial turbulent structures are chosen six times larger at M c = 2.2 when compared to other M c to be consistent with PG simulation. It is nevertheless shown in Appendix A that the size of initial turbulent structures does not influence the growth rate during self-similarity. This choice was motivated by the will to shorten the simulation. Enlarging the size of initial turbulent structures accelerates the unstable growth phase. As a consequence, in figure 6, M c = 1.1 and M c = 2.2 curves overlap after τ ≈ 2500. Slopes and standard deviation computed over the self-similar range are given in figure 6. At M c = 0.1, because of the suppression of compressibility effects, growth rate is very close to that of PG flow: the difference is about 1.5% and is below the standard deviation range. Like for PG, the DNS at M c = 0.1 is considered as the reference incompressible case and is used to plot the dependence of the normalised momentum thickness growth rate with respect to M c . At M c = 1.1, comparison between DG and PG flows is detailed in [START_REF] Vadrot | Analysis of turbulence characteristics in a temporal dense gas compressible mixing layer using direct numerical simulation[END_REF] during unstable growth and self-similar phases.

Figure 6 shows that the momentum thickness growth rates are very close between M c = 2.2 and M c = 1.1 unlike the perfect gas case. The well-known decrease of the growth rate with the convective Mach number is modified by dense gas effects. Despite being a highly compressible fluid, compressibility effects decrease in F C -70. Explanations for this effect are given in section 5.

Slopes provided in figure 6 are determined using the same methodology used for PG in section 3.1. For each convective Mach number, the non-dimensional integrated turbulent production term P * int is plotted over time. The three main phases described for the PG flow can also be identified for DG. One can notice that, at M c = 2.2, the initial phase corresponding to an energy transfer to the most unstable modes is much shorter for DG flow, likely because unstable modes are different between the two types of gas. After this phase, turbulent production reaches a maximum which decreases as M c increases. Finally, self-similar periods are defined selecting the range during which turbulent production is almost constant. As observed for PG flow, the self-similar period extends as M c increases. One can also notice that integrated production terms in DG flows are consistent with momentum thickness growth rates: the values of P * int are very close between M c = 2.2 and M c = 1.1 and the value of P * int at M c = 0.1 is twice larger than the one at M c = 1.1. This observation confirms the relevance of [START_REF] Vreman | Compressible mixing layer growth rate and turbulence characteristics[END_REF]

Comparison with perfect gas over the self-similar period

Self-similar periods have been selected for both types of gas. It is thus possible to plot the evolution of self-similar growth rates as a function of the convective Mach number. Slopes are usually normalised using an incompressible reference case at very low convective Mach number for which compressibility effects can be neglected. DNS at M c = 0.1 is considered here as the reference incompressible case. For example, [START_REF] Pantano | A study of compressibility effects in the high-speed turbulent shear layer using direct simulation[END_REF] use a simulation at M c = 0.3 as a reference case. There is no consensus on this choice, which can partly explain the spreading of PG results observed in figure 8 -where the same literature results used in figure 5 are reported. DG mixing layer results are plotted with error bars coloured in black. They represent the standard deviation of the normalised growth rate over the self-similar range. Unlike the PG mixing layer which shows a fairly abrupt decrease of its growth rate as M c increases, the DG mixing layer seems to be much less influenced by compressibility effects as M c becomes larger than 1.1. Differences between DG and PG mixing layers are large enough when compared to standard deviations to reveal that turbulence development is actually modified by dense gas effects in mixing layer flows.

In order to analyse the impact of compressibility effects, [START_REF] Pantano | A study of compressibility effects in the high-speed turbulent shear layer using direct simulation[END_REF] study the TKE equation and particularly the importance of the turbulent production term. They find that this term is decreasing in consistent proportion with the growth rate as the convective Mach number increases. The computation of TKE equation terms requires to statistically average the terms. This can only be done during the self-similar A. Vadrot, A. Giauque and C. Corre period during which both mixing layers are in a statistically stable state. Figure 9 shows the comparison between DG and PG mixing layers of the normalised main terms of the TKE equation over the non-dimensional cross-stream direction y/δ θ (t). Production, dissipation and transport terms are averaged during corresponding self-similar ranges. The production term (denoted P ) is always positive and is responsible for the growth of the mixing layer. Viscous dissipation (denoted D) is always negative and counterbalances the production term. The transport term (denoted T ) enables the propagation of TKE from the center to the edges of the mixing layer. It is thus negative at the center and positive close to the edges. Consistently with the comparison of slopes between DG and PG flows, all main terms and particularly the production term are two to three times larger for DG.

Another noticeable feature which was highlighted in the previous analysis at M c = 1.1 [START_REF] Vadrot | Analysis of turbulence characteristics in a temporal dense gas compressible mixing layer using direct numerical simulation[END_REF]) is confirmed here: curves are wider for the PG mixing layer, when compared to the DG mixing layer. For the DG mixing layer, TKE is more localised at the center. This is directly linked to the thermodynamic profiles, which are wider for PG mixing layer (see Figure 19 in section 5.3).

Other terms of the TKE equation, namely the compressible dissipation, the mass-flux coupling term, the convective derivative of the TKE and even the pressure dilatation are negligible for both types of gas. The pressure dilatation term which is directly linked to shocklet effects is carefully analysed in section 5.1 to quantify shocklet effects on the mixing layer growth.

As mentioned in the introduction, [START_REF] Pantano | A study of compressibility effects in the high-speed turbulent shear layer using direct simulation[END_REF] demonstrate that the compressibility-related reduction of the momentum thickness growth rate is induced by the reduction of pressure-strain terms Π ij , which causes a reduction of turbulent production. In the TKE equation, which is obtained from the sum R ii , the pressurestrain terms do not appear. Their sum Π ii , which constitutes the pressure-dilatation term, appears in the TKE equation but is negligible. In order to study pressure-strain terms, one needs to plot turbulent stress tensor equations terms. Figure 10 shows the main terms of the x-and y-components of the turbulent stress tensor equations. In the streamwise direction, the pressure-strain term counterbalances the streamwise production, whereas in the cross-stream directions, pressure-strain term is and is balanced by viscous dissipation. In the cross-stream direction, turbulent production term can be neglected unlike in the streamwise direction for which it is maximal.

One can notice that pressure-strain terms are significantly reduced for PG flows compared to DG flows at M c = 2.2: streamwise pressure strain term is twice larger for DG when compared to PG. This is consistent with the comparison of momentum thickness growth rates. For both types of gas, growth rates are identically linked to their pressure-strain terms. Compressibility effects impact the same terms for both DG and PG.

It remains to verify the last step in [START_REF] Pantano | A study of compressibility effects in the high-speed turbulent shear layer using direct simulation[END_REF]'s explanation, which is that the reduction of pressure-strain terms is caused by a reduction of normalised pressure fluctuations. Figure 11 shows the cross-stream distribution of the root-mean squared value of pressure normalised by the dynamical pressure 1 2 ρ 0 (∆u) 2 . Comparison is made between DG and PG flows at M c = 1.1 and M c = 2.2.

At M c = 1.1, DG and PG distributions are very close as are their corresponding momentum thickness growth rates. As the convective Mach number increases, DG non-dimensional pressure fluctuations experience a 20% decrease also consistent with the observed decrease in the growth rate. This decrease is yet much smaller than that of the PG mixing layer, in which normalised pressure fluctuations are approximately divided by a factor of two. To sum up, although the same mechanism is responsible for the growth rate decrease in both types of gas (i.e. the reduction of non-dimensional pressure fluctuations), its effect is significantly different between the two types of gas. For DG flows, the well-known compressibility-related reduction of the momentum thickness growth rate is almost suppressed by dense gas effects at convective Mach numbers above M c = 1.1.

Figure 12 shows the comparison between PG and DG streamwise specific turbulent kinetic energy spectra computed over the centreline. Spectra are normalised by (∆u) 2 δ θ (t) in the same way as [START_REF] Pirozzoli | Early evolution of the compressible mixing layer issued from two turbulent streams[END_REF] and averaged over the self-similar period. The longitudinal Taylor microscale λ x is also indicated for each gas in figure 12. Its value is much larger for DG flow consistently with Reynolds numbers computed from Taylor microscales given in table 3. The inertial phase is thus significantly reduced for the PG flow. Dissipation occurs at much larger scales making the comparison between the two inertial phase slopes difficult. Spectra yet confirm previous results observed at M c = 1.1 [START_REF] Vadrot | Analysis of turbulence characteristics in a temporal dense gas compressible mixing layer using direct numerical simulation[END_REF]): dense gas effects tend to increase small scales energy. The dissipation term, which is the main term at these scales, is significantly reduced. These results suggest a need for a specific sub-grid scale modelling of dense gas flows the small-scale dynamics of which is significantly modified with respect to perfect gas flows.

Analysis of discrepancies between DG and PG

First hypothesis: Effect of shocklets

The previous analysis conducted at M c = 1.1 shows that the growth rate is not influenced by the dense gas effect during the self-similar period [START_REF] Vadrot | Analysis of turbulence characteristics in a temporal dense gas compressible mixing layer using direct numerical simulation[END_REF]). However, significant differences are observed during the unstable growth phase. At M c = 1.1, the evolution of the turbulent Mach number shows that shocklets might be detected during the unstable growth phase but not during the self-similar range, during which M t decreases well below the range of values for which shocklets are expected. It is known that the generation of shocklets is different between BZT DG flow and PG flow [START_REF] Giauque | Direct numerical simulations of forced homogeneous isotropic turbulence in a dense gas[END_REF]), yet can shocklets alone explain discrepancies between DG and PG flows ?

In the current analysis, we increase the convective Mach number in order to reach larger turbulent Mach numbers during the self-similar period and to analyse the influence of shocklets. Figure 13 shows the temporal evolution of the turbulent Mach M t (see equation 1.4). Turbulent Mach numbers increase the initial phase up to 0.9 respectively for DG and flows. Then M t decreases and reaches a rather plateau corresponding to the self-similar period. During this phase, average values of turbulent Mach numbers are respectively equal to 0.67 and 0.49 for DG and PG flows. Shocklets can thus be observed during both DG and PG self-similar periods.

In order to study their effect on the growth rate, one can analyse the compressible component of the dissipation given in Equation 2.11. [START_REF] Zeman | Dilatation dissipation: the concept and application in modeling compressible mixing layers[END_REF] and Sarkar et al. (1991) show that the dilatational part of the dissipation increases with the turbulent Mach number because of the occurrence of eddy shocklets in the compressible regime. [START_REF] Wang | Effect of flow topology on the kinetic energy flux in compressible isotropic turbulence[END_REF] perform compressible isotropic turbulence simulations and observe that shocklets act as kinetic energy sinks which absorb large-scale kinetic energy. Shocklets are thus an additional source of dissipation. The dilatational dissipation is computed over the self-similar period. Figure 14 shows the ratio between the compressible and the total dissipation rate over the cross-stream direction. Around y/δ θ (t) ≈ 3.5, one can note an increase of the ratio. It corresponds to the borders of the mixing layer, outside of which the dissipation drops to zero (see Figure 9). Except for these regions, at M c = 1.1, the compressible dissipation represents less than 0.5% of the total dissipation for both DG and PG flows. At M c = 2.2, the ratio increases consistently with the increase of turbulent Mach numbers. The ratio is thus larger for DG flow compared to PG flow. However, the rate of dilatational dissipation with respect to the total dissipation remains below 4% for DG and below 1% for PG. Compressible dissipation can therefore be neglected with respect to the total dissipation. Shocklets have a limited influence on the TKE equation. Since the TKE equation governs the mixing layer dynamics, one cannot explain discrepancies observed between DG and PG flows with shocklets effect.

Additional simulations varying the initial thermodynamic operating point

In order to explain discrepancies observed between DG and PG flows, we perform additional DNS varying the initial thermodynamic operating point. Figure 15 shows the four selected operating points. DGA corresponds to the reference simulation analysed in section 4. DGA's initial operating point is located inside the inversion zone also called BZT region. The operating point of the second simulation DGB is chosen outside the inversion region and inside the dense gas zone. This enables us to investigate the impact of BZT effects on the mixing layer growth. Finally, for DGC and DGD, initial operating points are chosen on the same adiabatic curves as respectively DGB and DGA but outside the dense gas zone. The diversity of targeted thermodynamic regions aims at providing a proper insight into the effects of dense gas on the shear layer growth.

At first, one needs to validate DNS named DGB, DGC and DGD. Table 4 gives simulations parameters including r, l x , and l z for the four different simulations. Achieved values are very close to DGA and since DGA has been validated previously (see section 4 and Appendix A), one can consider that DGB, DGC and DGD are adequately resolved. The size of computational domains have been enlarged for DGC and DGD in the y direction in order to provide the mixing layer with more space in order to reach self-similarity. Self-similar periods are defined for each DNS using the same methodology previously presented in section 3.1. Plateaus showing constant integrated turbulent production correspond to self-similar periods. They are identified with vertical lines in figure 16. In addition, beginning and ending times are given in the caption for each case. Although all the DNS are performed at the same convective Mach number M c = 2.2, results are quite different. The initial evolution is similar, but after τ ≈ 1100, discrepancies appear especially for DGD. Maximum values and self-similar regimes are influenced by the initial thermodynamic operating point.

The comparison of mixing layer momentum thickness evolutions is done in figure 17. Slopes with standard deviations computed during self-similar regimes are indicated on the plot. From these results, one can deduce that BZT region does not have a major influence on the mixing layer growth. DGC's growth rate is indeed very close to DGA's although initial thermodynamic operating points are located respectively outside and inside DG and BZT regions. The relation between the mixing layer growth and the initial thermodynamic operating point is not obvious: operating points located on the same adiabatic curve (respectively DGA, DGD and DGB, DGC) are far away in terms of growth rate. Looking at the growth rate, simulations can be classified by pairs: DGA goes with DGC and DGB goes with DGD. One can observe that slopes are all below the M c = 1.1 growth rate. It means that the well-known compressibility-related reduction of the momentum thickness growth rate is still verified. Yet there is an additional effect due to the initial thermodynamic operating point. At the end of section 4, the physical explanation provided by [START_REF] Pantano | A study of compressibility effects in the high-speed turbulent shear layer using direct simulation[END_REF] was assessed on DGA: the reduction of the momentum thickness is due to a reduction of normalised pressure fluctuations. It remains to check whether this reduction of normalised pressure fluctuations is also observed for DGB, DGC and DGD. Figure 18 shows the A. Vadrot, A. Giauque and C. Corre normalised growth rate as a function of the normalised pressure fluctuations computed at the center of the mixing layer. For PG flow, the reduction is significant. Between M c = 1.1 and M c = 2.2, growth rate and normalised pressure fluctuations are divided by a factor of two. For DG, the decrease of the normalised growth rate is also correlated with a decrease of pressure fluctuations. Among cases at M c = 2.2, the ranking purely based on the level of pressure fluctuations is not entirely satisfactory but this could be explained by standard deviations caused by variations of the plateaus of integrated turbulent production. Moreover, other effects must also be taken into account for dense gases: this is the topic of the next section.

Analysis of discrepancies between DG and PG flows

There is a significant effect of dense gas on the well-known compressibility-related reduction of the momentum thickness growth rate. Dense gas effects modify the decrease at convective Mach numbers larger than M c = 1.1. Between M c = 1.1 and M c = 2.2, the growth rate slope does not vary much for DG. Several factors can be identified, which contribute to explain the observed discrepancies between DG and PG mixing layers. The first main difference between DG and PG flows is the ratio between the enthalpy and the kinetic energy. It is associated to the Eckert number, which is defined for the mixing layer as:

Ec = (∆u) 2 c p0 T 0 (5.1)
where c p0 denotes the initial specific heat capacity at constant pressure and T 0 , the initial temperature. Initial Eckert numbers are computed for each DNS and results are gathered in table 5. For DG flows, values are about two orders of magnitude lower than PG flows. Two features of DG mixing layers are responsible for these significant differences: the large heat capacity of FC-70 and the small differential speed ∆u. The differential speed is defined in order to get the same initial convective Mach number between DG and PG mixing layers. Since the sound speed is much lower in dense gases, a much lower differential speed is obtained for a given value of the convective Mach number, which mechanically reduces the Eckert number. With small Eckert numbers, kinetic energy becomes negligible when compared to the enthalpy or to the internal energy †. It is the case for all DG flows in this study even though the convective Mach number is large. As shown by the present results, kinetic energy also decouples from thermodynamics compressibility effects and the growth rate of the momentum thickness is allowed to reach larger values. It can be observed that the close values of the momentum thickness growth rates for DGA / DGC on one hand and DGB / DGD on the other hand are well correlated with the values of the initial Eckert number reported in table 5. The lower Eckert numbers for DGA / DGC correspond to higher growth rates for these shear layer configurations, induced by an even stronger decoupling between internal and kinetic energy for DGA/DGC with respect to DGB/DGD. However, the Eckert number can not be the only factor explaining dense gas effect on the growth rate since DGC displays a slightly lower growth rate with respect to DGA, with a slightly lower value of the initial Eckert number.

For DG flows, the amount of internal energy is much larger when compared to kinetic energy. Internal and kinetic energies are decoupled in that case. In the equation of † At the initial conditions γ is about 1.3 and internal energy and enthalpy are of the same order of magnitude. The evolution of the average density confirms this reduction. The PG air density suffers a 40% decrease at the center between M c = 1.1 and M c = 2.2. In the PG, friction heating is important and leads to an increase of the temperature, which induces a decrease of the density. The mechanism is significantly reduced in dense gas flows. For DG, the temperature is almost constant and averaged density displays very limited variations. At M c = 2.2, the averaged density decrease at the center of the mixing layer represents about 8% of the initial density compared to 45% for air. Equation 3.1 shows that this effect influences the mixing layer growth rate which depends on the density. As the mixing layer develops in PG, strong friction occurs at the center, which decreases the density. The momentum thickness growth rate is thus significantly reduced for PG when compared to DG.

Mc

Figure 19c displays the root mean square value of density fluctuations. Between PG and DG flows, the distribution across the mixing layer changes shape. For PG, it consists in two symmetric peaks with respect to the center of the mixing layer. Peaks are located at the borders of the mixing layer, where the cross-stream gradient of averaged density is maximal. In this region, the mixing layer flow experiences strong dynamic and thermal variations with an important coupling between internal and kinetic energy. For DG, the distribution is composed of a single peak located at the center of the mixing layer. The distribution is much less affected by the variation of the averaged density. For DG, thermal quantities are less influenced by flow dynamics because of the decoupling of internal energy and kinetic energy. The root mean square value of density fluctuations diffuses from the center of the mixing layer.

The amplitudes of the distributions are also quite different between DG and PG flows. For DG, the maximum root mean square value of density fluctuations is multiplied by a A. Vadrot, A. Giauque and C. Corre factor of three from M c = 1.1 to M c = 2.2. In the PG case, it is multiplied by a factor of about two. Compressible flows are more subject to root mean square density fluctuations which increase as the Mach number grows. An explanation can be found in the definition of the isentropic compressibility coefficient, which is large for DG flows:

χ s = 1 ρ ∂ρ ∂p s (5.2)
For flows with large values of χ s , small variations of pressure lead to large variations of density. The sound speed is directly linked to the isothermal compressibility since:

c = 1 √ ρχ s (5.3)
For DG flows, the large isentropic compressibility factor strongly diminishes the sound speed. As a result, the initial sound speed in the computed DG flows is about six times smaller when compared to its initial value for the PG shear layers. Figure 20 shows the normalised momentum growth rate at M c = 2.2 as a function of the normalised sound speed. A rather clear correlation appears between the momentum thickness growth rate and the initial sound speed: the growth rate decreases with increasing sound speed.

The main conclusion that can be drawn from these observations is that the smaller Eckert number in DG flows causes a decoupling between internal and kinetic energy and induces less friction heating. Both phenomena influence the mean and fluctuating thermal physical quantities, which consequently limits the compressibility-related reduction of the momentum thickness growth rate.

remarks

The present work extends the previous analysis of a temporal compressible shear layer conducted at c = 1.1 [START_REF] Vadrot | Analysis of turbulence characteristics in a temporal dense gas compressible mixing layer using direct numerical simulation[END_REF]) a larger convective Mach number M c = 2.2 air described as perfect gas and FC-70 (BZT gas) described using Martin Hou EoS. A reference incompressible DNS is also performed at M c 0.1 to provide the incompressible growth rate δθ,inc used to normalise the growth rate δθ . The computed evolution of the mixing layer growth rate with respect to the convective Mach number is compared with available results from the literature for perfect gas. The PG results are found consistent with the literature and establish the accuracy of the present simulations.

The choice of the domain size is paramount in this study. The domain is enlarged at M c = 2.2 for both DG and PG DNS when compared to DNS at M c = 1.1 in order to ensure mixing layers reach self-similarity. An analysis presented in Appendix A is performed to thoroughly investigate the sensitivity of the DG mixing layer to domain extent and to the size of initial turbulent structures. Results establish the relevance of the choices of domain extent and initial structures size made in the present study.

The selection of the self-similar period is a key point in the study of mixing layers: this choice is complex and the diversity of criteria used for the selection process contributes to the scattering of δθ / δθ,inc = f (M c ) plots reported in the literature. In the present work, self-similar periods are selected using the integrated streamwise production over time, which is proportional to the momentum thickness growth rate under certain conditions [START_REF] Vreman | Compressible mixing layer growth rate and turbulence characteristics[END_REF].

The comparison between perfect and dense gases shows major differences for the momentum thickness growth rates at M c = 2.2. The dense gas flow limits the wellknown compressibility-related reduction of the momentum thickness growth rate. At M c = 2.2, the growth rate is twice as large for dense gas when compared to perfect gas. [START_REF] Pantano | A study of compressibility effects in the high-speed turbulent shear layer using direct simulation[END_REF] demonstrate that for perfect gas flows the growth rate Figure 21 shows temporal evolutions of momentum thickness for the simulations listed in Table 6. DG1 is performed with the same domain lengths and size of initial turbulent structures (relatively to the initial momentum thickness) as in the previous M c = 1.1 study DG0. At τ = 4000, self-similarity is not yet achieved but flow field visualisations indicate that the y boundaries of the domain are reached. DG2 is then conducted with a domain size doubled in the y direction and with smaller initial turbulent structures corresponding to L x /4 = 86δ θ,0 , in order to speed up the mixing layer development. Simulations show that the modification of initial structures size only modifies the time necessarily to reach the unstable growth phase but not the growth rate itself. where i denotes the direction of the velocity.

Figure 23 shows the evolution of the two-point correlation over the streamwise direction for the three velocity components. In the PG case, the correlation increases significantly for the x-and y-velocity components as M c increases. As noticed in [START_REF] Freund | Compressibility effects in a turbulent annular mixing layer. Part 1. Turbulence and growth rate[END_REF] and [START_REF] Matsuno | Compressibility effects in high speed turbulent shear layers-revisited[END_REF], at high compressible regimes, eddies are stretched in the streamwise direction. In the DG case, the correlation stays approximately the same between M c = 1.1 and M c = 2.2 for the three components, except for the x-component which is slightly larger for M c = 2.2 when compared with M c = 1.1. Consistently with Figure 8, which shows that the mixing layer growth rate is slightly affected by the convective Mach number from M c = 1.1 to M c = 2.2, the structure of eddies stays approximately the same in the streamwise direction unlike for PG flows. One can also notice that for all cases, the correlation drops to a low value at r x /L x = 0.5 which confirms that the streamwise domain length is sufficiently large.

Figure 24 shows some snapshots of the velocity magnitude. As noticed in Figure 23, the size of turbulent structures increase from M c = 1.1 to M c = 2.2 in PG flow unlike in DG flow, where the size of turbulent structures remains stable between M c = 1.1 to M c = 2.2. At M c = 1.1, there is no difference between DG and PG flow field visualisation. Consistently with the evolution of the normalised momentum thickness growth rate as a function of the convective Mach number (Figure 8), differences appear at M c = 2.2 in highly compressible regime. The Prandtl number is defined as the ratio between the kinematic viscosity and the thermal diffusivity (α):

P r = ν α (C 1)
A large Prandtl number indicates that the viscous diffusivity is faster than the thermal diffusivity. It would thus affect the temperature distribution. Table 7 gives the values of Prandtl numbers and mixing layer growth rates for each DNS. The growth rate of the DGC DNS is approximately the same as the one of the DGA DNS whereas the Prandtl number of the DGC DNS is twice larger than the one of the DGA DNS. Results show that there is no correlation between the Prandtl number and the mixing layer growth rate. However, the influence of the Prandtl number can be seen in Figure 25. Although the rate of normalized temperature fluctuations is very low in the DG DNS (below 6.10 -3 ) compared with the PG DNS (about 1.7), one can notice differences when varying the initial thermodynamic operating point. The larger the Prandtl number, the larger the temperature fluctuations.

Figure 1 .

 1 Figure 1. The initial thermodynamic state is represented in the non-dimensional p -v diagram for BZT dense gas FC-70 at Mc = 2.2. The dense gas zone (Γ < 1) and the inversion zone (Γ < 0) are plotted for the Martin-Hou equation of state. pc and vc are respectively the critical pressure and the critical specific volume. The initial value of the fundamental derivative of gas dynamics is equal to Γ initial = -0.284. The normalised distribution of the thermodynamic states is plotted at the beginning of the self-similar period (τ = 4000) along the curve where the distribution of thermodynamic states is the largest.

Figure 2 .

 2 Figure 2. Configuration of the temporal mixing layer. The velocity magnitude is plotted for the DG DNS at Mc = 2.2 at τ = 4000.

Figure 3 .

 3 Figure 3. Temporal evolution of the mixing layer momentum thickness for Mc = 0.1/1.1/2.2 using air with PG EoS. Slopes are non-dimensional and standard deviations computed over the self-similar period are indicated on the plot.

Figure 4 .

 4 Figure 4. Temporal evolution of the non-dimensional streamwise turbulent production term integrated over the whole domain P * int = (1/(ρ0(∆u) 3 )) Ly ρPxxdy (with

Figure 5 .

 5 Figure 5. Evolution of the mixing layer growth rate with respect to the convective Mach number for air using PG EoS. Comparison is made with available DNS results in literature and experimental results by Rossmann et al. (2001). Standard deviations are indicated on the plot.

Figure 6 .

 6 Figure 6. Temporal evolution of the mixing layer momentum thickness for DG at Mc = 0.1 -1.1 -2.2.

Figure 7 .

 7 Figure 7. Temporal evolution of the non-dimensional streamwise turbulent production term integrated over the whole domain P * int = (1/(ρ0(∆u) 3 )) Ly ρPxxdV (with

Figure 8 .

 8 Figure 8. Evolution of the mixing layer growth rate over the convective Mach number for air and for FC-70. Comparison is made with available DNS results in literature and experimental results in Rossmann et al. (2001).

Figure 9 .

 9 Figure 9. Distribution of the volumetric normalised powers over the non-dimensional cross-stream direction y/δ θ (t) at Mc = 2.2. P: Production, D: Dissipation and T: Transport are normalised by ρ0(∆u) 3 /δ θ (t). Distributions have been averaged between the upper and the lower stream to get perfectly symmetrical distributions.

Figure 10 .

 10 Figure 10. Distribution of the main non-dimensional volumetric power terms of the x-(top) and y-(bottom) turbulent stress tensor (Rxx and Ryy) equations over the non-dimensional cross-stream direction y/δ θ (t). Pxx and Pyy: Streamwise and cross-stream production, Πxx and Πyy: Streamwise and cross-stream pressure-strain and Dxx and Dyy: Streamwise and cross-stream dissipation terms are normalised by ρ0(∆u) 3 /δ θ (t). Results are computed at Mc = 2.2. Distributions have been averaged between the upper and the lower stream to get perfectly symmetrical distributions.

Figure 11 .

 11 Figure 11. Distributions of the root mean square value of pressure averaged over the self-similar period, plotted along the y direction and compared between FC-70 and Air at Mc = 1.1 and Mc = 2.2. Distributions have been averaged between the upper and the lower stream to get perfectly symmetrical distributions.

Figure 12 .

 12 Figure 12. Streamwise specific TKE spectra computed at the centreline.

Figure 13 .

 13 Figure 13. Temporal evolution of the turbulent Mach number Mt.

Figure 14 .

 14 Figure 14. Distributions of the ratio between the compressible dissipation ( d ) and the total dissipation ( ) (see details in Equations 2.10 and 2.11). Results are averaged over the self-similar period. Comparison is made between FC-70 and Air at Mc = 1.1 and Mc = 2.2. Distributions have been averaged between the upper and the lower stream to get perfectly symmetrical distributions.

Figure 15 .

 15 Figure 15. Four different initial thermodynamic states used to perform additional DNS are represented in the non-dimensional p -v diagram for BZT dense gas FC-70 at Mc = 2.2. The dense gas zone (Γ < 1) and the inversion zone (Γ < 0) are plotted for the Martin-Hou equation of state. pc and vc are respectively the critical pressure and the critical specific volume.

Figure 16 .

 16 Figure 16. Temporal evolution of the non-dimensional streamwise turbulent production terms integrated over the whole domain P * int = (1/(ρ0(∆u) 3 )) Ly ρPxxdV (with

Figure 17 .

 17 Figure 17. Temporal evolution of the mixing layer momentum thickness for DG at Mc = 2.2. Results are shown for the FC-70 for four different DNS: DGA, DGB, DGC and DGD.

Figure 18 .

 18 Figure 18. Evolution of the non-dimensional mixing layer growth rate over the center root-mean squared value of pressure normalised by 1 2 ρ0(∆u) 2 . Results are given for DG and PG at Mc = 1.1 and Mc = 2.2.

Figure 19 .

 19 Figure 19. The non-dimensional Reynolds averaged temperature (a) and density (b); and root mean squared value of the density (c) are averaged over the self-similar regime and plotted along the y direction. Comparison is made between FC-70 and Air at Mc = 1.1 and Mc = 2.2.

Figure 20 .

 20 Figure 20. of the non-dimensional mixing layer growth rate as a function of the sound speed normalised with pc/ρc. Results are given for DG and PG = 2.2.

Figure 21 .

 21 Figure 21. Temporal evolution of the mixing layer momentum thickness.

Figure 22 .

 22 Figure 22. Temporal evolution of the integral length scale lz.

Figure 23 .Figure 24 .

 2324 Figure 23. The streamwise two-point correlations of the (a) x-, (b) y-and (c) z-velocity component at the beginning of the self-similar period. Comparison is made between FC-70 and Air at Mc = 1.1 and Mc = 2.2.

Figure 25 .

 25 Figure 25. The root mean squared value of the temperature are averaged over the self-similar regime and plotted along the y direction. Comparison is made between DG simulations at Mc = 2.2.

Table 1 .

 1 Physical parameters of FC-70

	70 608.2	10.2	0.270 488.2	118.7	0.493

Table 2 .

 2 Simulation parameters. Lx, Ly and Lz denote computational domain lengths measured in terms of initial momentum thickness. Nx, Ny and Nz denote the number of grid points. L0 denotes the size of initial turbulent structures (k0 = 2π/L0) measured in terms of initial momentum thickness. All grids are uniform.

		Re δ θ,0 Lx × Ly × Lz	Nx × Ny × Nz	∆u (m.s -1 ) δ θ,0 (nm) L0
	Air 0.1 1.0	160	344 × 344 × 86 1024 × 1024 × 256	34.11	135.8	Lx/48
	Air 1.1 1.0	160	344 × 172 × 86 1024 × 512 × 256	375.2	12.35	Lx/48
	Air 2.2 1.0	160 688 × 688 × 172 1024 × 1024 × 256	753.0	6.153	Lx/8
	FC-70 0.1 1.0	160	344 × 344 × 86 1024 × 1024 × 256	5.665	2070	Lx/48
	FC-70 1.1 1.0	160	344 × 172 × 86 1024 × 512 × 256	62.32	188.2	Lx/48
	FC-70 2.2 1.0	160 688 × 344 × 172 1024 × 512 × 256	125.1	93.77	Lx/8

Table 3 .

 3 's integral length scale reaches Mc Re δ θ Re λx r = Lη/∆x lx/Lx lz/Lz Non-dimensional parameters computed at the beginning and at the end of the self-similar period for Mc = 2.2 simulations. Re λx denotes the Reynolds number based on the longitudinal Taylor microscale λx = 2u Lη denotes the Kolmogorov length scale computed at the centreline.

	Air (τ = 700)	0.1 1879 209	0.63	0.10 0.04
	Air (τ = 1450) 0.1 3444 194	0.81	0.11 0.13
	FC-70 (τ = 550) 0.1 1448 135	0.58	0.04 0.05
	FC-70 (τ = 900) 0.1 2176 201	0.7	0.07 0.06
	Air (τ = 1700) 1.1 1874 143	0.97	0.07 0.06
	Air (τ = 2550) 1.1 2413 156	1.09	0.12 0.08
	FC-70 (τ = 1700) 1.1 2469 176	0.80	0.09 0.05
	FC-70 (τ = 2550) 1.1 3304 241	0.87	0.20 0.05
	Air (τ = 11500) 2.2 3487 146	1.44	0.12 0.07
	Air (τ = 14100) 2.2 3700 191	1.64	0.11 0.10
	FC-70 (τ = 4000) 2.2 4663 263	0.52	0.10 0.06
	FC-70 (τ = 6000) 2.2 6259 390	0.57	0.16 0.05

2

x /(∂u x /∂x) 2 computed at the centreline.

Table 5 .

 5 Eckert numbers and normalised momentum thickness growth rates are given for each simulation.energy conservation (Equation 2.7), all the terms can be neglected with respect to the temporal and convective internal energy terms. Since the Eckert number quantifies the friction heating, it is significantly reduced in DG flows as previously shown by[START_REF] Gloerfelt | Dense-gas effects on compressible boundary-layer stability[END_REF]. Figure19shows the distribution of the Reynolds averaged temperature, density and the root mean square value of density fluctuations over the cross-stream direction of the shear layer. Results are averaged over the self-similar period. It can be observed in figure19that temperature variations are almost suppressed for DG.Sciacovelli et al. (2017a) confirm this remark in supersonic turbulent channel flows and state that dense gas flow are less subject to friction losses associated with Mach number effects. For the mixing layer, above M c = 1.1, compressibility effects associated with the increase of convective Mach number have less influence on DG flows in part because of the reduction of friction heating.

	Ec	δθ / δθ,inc
	DG 1.1 0.0040 0.484
	DGA 2.2 0.0162 0.395
	DGB 2.2 0.0226 0.352
	DGC 2.2 0.0147 0.389
	DGD 2.2 0.0203 0.342
	PG 1.1 1.94	0.450
	PG 2.2 7.74	0.188

Table 7 .

 7 Prandtl numbers and normalized momentum thickness growth rates are given for each DNS.

	Mc P r δθ / δθ,inc
	DG 1.1 1.52 0.484
	DGA 2.2 1.52 0.395
	DGB 2.2 1.41 0.352
	DGC 2.2 3.33 0.389
	DGD 2.2 6.41 0.342
	PG 1.1 0.71 0.450
	PG 2.2 0.71 0.188

† Email address for correspondence: aurelien.vadrot@ec-lyon.fr

Yet, a large decrease of the growth rate is observed for DG2 around τ = 4000; selfsimilarity cannot be reached. Figure22displays the time evolution of the integral length scale in the z direction l z for DG2 and DG3 simulations. Around τ = 4000, the integral length scales l z /L z suddenly decreases for DG2 after having reached a value of 0.2. The domain is thus not large enough to account for spanwise turbulent structures, which causes a growth rate decrease and prevents the transition to self-similarity.Because of the aforementioned observations, domain sizes have been doubled in x and z directions when compared to DG1. This corresponds to the DG3 simulation, which is the reference DNS used in Section 4 to compare results between DG and PG. For DG3, the momentum thickness evolution reaches a perfectly linear stage and self-similarity is well achieved as confirmed by Figures 7 (right) and 21.

Acknowledgements -This work is supported by the JCJC ANR EDGES project, grant #ANR-17-CE06-0014-01 of the French Agence Nationale de la Recherche. Simulation have been carried out using HPC resources at CINES under the project grant #A0062A07564.

Declaration of Interests -The authors report no conflict of interest.

A. Vadrot reduction is due to the reduction of pressure fluctuations leading to the reduction of pressure-strain terms. We show that growth rate is also correlated with pressure fluctuations in dense gas flows. Yet, the small scales dynamics is very different. A much larger dissipation is also observed for perfect gas mixing layer. These results call for a specific sub-grid scale modelling for dense gas flows when simulated using Large Eddy Simulation. Additional dense gas DNS have been performed at three others initial thermodynamic operating points. Results show that BZT effects have only a small impact on the mixing layer growth. Shocklets indeed produce only a limited effect on mixing layer growth. The compressible dissipation is negligible when compared with the total dissipation. For dense gas mixing layers, several physical factors tend to reduce compressibility effects: the decoupling of kinetic and internal energy reduces the effect of increasing M c ; reduced friction losses in dense gas flows modify the distribution of the averaged density, which therefore favours the momentum thickness growth rate. Finally, it is found that increasing the initial isothermal compressibility also increases the momentum thickness growth rate in dense gas flows. Initial sound speed could therefore be an appropriate indicator when forecasting the mixing layer growth rate in real gas flows. Note that the perfect gas results are restricted to air, with heat capacity ratio equal to γ = 1.4. Further exploration could investigate the effect of γ close to unity over the perfect gas results and provide a comparison with the dense gas results in order to separate possible γ-effects from dense gas effects.

Appendix A. DG mixing layer: Influence of domain size, resolution and initial turbulent structures size

Additional simulations have been performed for DG mixing layer with Re δ θ,0 = 160 and M c = 2.2 in order to confirm proper resolution and domain size. The computational parameters corresponding to these simulations are summarized in Table 6 along with the parameters used in the previous study at M c = 1.1.

Appendix B. Analysis of spatial correlations.

This section is devoted to the analysis of two-point spatial correlations of the velocity components. Both the PG and DG flows are analysed and compared. In the streamwise


Khouloud Cheour

Mustapha Assarar

Daniel Scida

Rezak Ayad

M Asssarrar

X-L Gong

Xiao-Lu Gong

Effect of stacking sequence on the mechanical and damping properties of flax-glass fibre hybrid composites

Keywords: Damping, flax fibre, glass fibre, hybrid composite, mechanical properties

à la diffusion de documents scientifiques de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

Introduction

Synthetic fibres reinforced composites are used nowadays in various fields, since they offer many advantages. But their use raises environmental problems at the end of life cycle. The aim of this work is to investigate the effect of adding flax fibre reinforcement to the glass fibre composites. Several works showed the interesting specific mechanical and dynamical properties, several properties of the natural fibre reinforcement [START_REF] Wambua | Natural fibres: can they replace glass in fibre reinforced plastics?[END_REF][START_REF] Murali Mohan Rao | Extraction and tensile properties of natural fibers: Vakka, date and bamboo[END_REF][START_REF] Ren | Effects of fluctuation of fibre orientation on tensile properties of flax sliver-reinforced green composites[END_REF][START_REF] Genc | Comparison of the mechanical properties of flax and glass fiber composite materials[END_REF]. But they have an hydrophilic behavior. So to reduce the significant lost of mechanical properties, several studies proposed the hybridization as a solution of both problems [START_REF] Nunna | A review on mechanical behavior of natural fiber based hybrid composites[END_REF][START_REF] Jawaid | Cellulosic/synthetic fibre reinforced polymer hybrid composites: a review[END_REF][START_REF] Sathishkumar | Hybrid fiber reinforced polymer compositesa review[END_REF][START_REF] Swolfs | Fibre hybridisation in polymer composites: A review[END_REF]. The aim of this work is to show the intrest of glass-flax hybridization, since it can provide good mechanical and dynamical properties, mass gain and also environmental advantages.

Experimental procedure 2.1 Materials

In this study, non-hybrid and hybrid composite materials with glass/flax fibres reinforcement and epoxy matrix were elaborated. All materials were manufactured by press platen process in the IFTS of Charleville-Mézières. The layers of dry fabrics were impregnated one by one with an epoxy resin SR1500 associated with hardener SD2503 (Sicomin Company). Once impregnated, fabrics were put on heat seal press under 6 bars pressure and 40°C temperature for three hours and a half. The nonhybrid plates were cut to obtain specimens of different lengths (240, 260 and 280 mm) and orientations (0°, 90° and 45°) and hybrid plates in only one orientation 0°. The width of both nonhybrid and hybrid specimens was 26 mm for all orientations and lengths.

Experimental protocol

The aim of the present work is to characterize dynamic properties of the studied materials using a free flexural vibrations analysis. The same equipment and method were used in a previous work of the authors [START_REF] Cheour | Effect of water ageing on mechanical and damping properties flax-fibre reinforced composite materials[END_REF]. The specimen is supported vertically by a rubber in order to have a free-free configuration. It is excited using an impulse hammer (PCB 086C03) in different points and the flexural vibrations are detected by an accelerometer (PCB 352C23) pasted in the specimen using wax. Both signals of excitation and response are digitalized by LMS SCADAS dynamic analyser. The frequency response enables us to identify the natural frequencies and loss factors from LMS PolyMAX method. More details can be found in [START_REF] Peeters | The PolyMAX frequency domain method: a new standard for modal parameter estimation[END_REF]. The results for each length, orientation and material showed in this work are the average of the five specimens.

Finite element analysis

To identify dynamic properties of the studied materials, specimens were modelled by finite element method. The element used in this work is shown in Figure 1. It is a four-node multilayered shell finite element, named S4, of the Abaqus software. This element is based on the first-order laminate theory. Transverse shear stresses could be neglected since thicknesses of specimens are widely lower than their lengths. Once natural frequencies are calculated by finite element analysis, corresponding values of stress and strain on the lower (l) and the upper (u) faces of each layer k of the element e (Figure 1), , ,

)

k e k e k e k e i      (1)
The energy

         (3)
S e is the surface of the finite element e. The evaluation of the loss factor deduced from the finite element analysis can be done with an energetic approach [START_REF] Assarar | Damping Analysis of Sandwich Composite Materials[END_REF]. The total energy stored in the structure can be expressed as follow:

11 22 66 U U U U    (4) With: elements, layers, , , ( 1, 2,6) ek ke ii ii U U i   (5)
Then, the dissipated energy with the damping phenomena in the layer k of the element e can be expressed as a function of loss factors as follow: The dissipated energy in the finite element e can be expressed by:

k k k U U U U        (6)
1 , n ee k k UU      (7)
The total energy dissipation U  in the whole structure is then:

Finally, loss factor estimated by finite element analysis can be deduced from:

U U    (9)
The previously cited formulations, are used to estimate the bending moduli and the loss factors of all the composites of the study.

Results and discussion 4.1 Elastic coefficient evaluation Elastic coefficient (E

L , E T and G LT ) of non-hybrid composites are deduced from free vibration tests of the specimens with different fibre orientations 0°, 45° and 90°. It's worthy to note that natural frequencies of the specimens with 0° and 90° fibre orientation depend on the longitudinal modulus E L and transverse modulus E T respectively. based on the natural frequencies obtained by experimental modal analysis, all moduli can be determined. Thereby for specimens with fibre oriented to 0°, the identification of the longitudinal modulus E L is done using an iterative method. A variation of modulus value in each iteration is done, until equality between natural frequencies obtained by experimental method and finite element analysis. This method is also used for 90° fibre orientation so transverse modulus E T can be identified. Then, shear modulus G LT is determined in the same way using results of the specimens with 45° fibre orientation and taking into account variations of both moduli E L and E T as a function of frequency. Regarding hybrid composites, the same procedure was used introducing for each layer the properties deduced for non-hybrid materials. The obtained results were then compared to experimental results. As explained in a previous work of the authors [START_REF] Cheour | Effect of water ageing on mechanical and damping properties flax-fibre reinforced composite materials[END_REF], it is possible to identify elastic moduli using natural frequencies. Table 1 resume the obtained values of the bending modulus of all hybrid and nonhybrid composites. The obtained results prove that introducing two internal flax layers to a glass fibre reinforced composite doesn't considerably affect the bending modulus of the stratified [G 2 /F 2 ] S since the difference is only about 3%. Natural fibre reinforced composites have a lightweight compared to glass fibre composites. That's why a comparison between specific modules (E/ρ) is done in this work. Figure 2 shows the evolution of the specific bending moduli as a function of a volume fraction of flax fibre in the hybrid composite. Unlike the tensile modulus that doesn't depend on the position of the layers in the hybrid [START_REF] Zhang | Tensile and interfacial properties of unidirectional flax/glass fiber reinforced hybrid composites[END_REF], the obtained results show that the stiffness of the composites depends on the stacking sequence. Regarding for example a composite with a flax fibre volume fraction about 25%, the specific modulus of the hybrid [G 2 /F 2 ] S increases by 32% and 61% compared to non-hybrid glass fibre and flax fibre reinforced composites respectively. On the other hand, the hybrid composite [F 2 /G 2 ] S has a lower specific modulus (27% and 11% lower than non-hybrid glass fibre and flax fibre materials). This result is explained by the fact that glass layers that have a better bending modulus are placed further of the medium plan (glass fibre composites E=43 GPa and flax fibre composite E=21 GPa).

Loss factor evaluation

Figure 3 compares the experimental results and those obtained by finite element method for the studied composites. Both methods show the same results, proving that the finite element method allows us to identify the loss factor of the hybrid materials. Figure 3 shows also that the variation of loss factor as a function of the frequency has two different behaviours. This difference is associated to the evolution of the distribution of the dissipated energies on the glass and flax layers as shown in Loss factor 0,6

[G 3 /F] s [G 2 /F 2 ] s [G/F 3 ] s Flax [F/G 3 ] s [F 2 /G 2 ] s [F 3 /G] s 0,8 1,0 1,2 1,4 1,6 1,8 2,0 2,2 Glass

Experimental result Modelisation

Volumique fraction of flax fibre (%) Flax layers

[L/V 3 ] s [L 2 /V 2 ] s [L/V 3 ] s [L 2 /V 2 ] s [L 3 /V] s a) b)
Glass layers

Flax layers

Glass layers

Flax layers

[G3/F]s [G3/F]s [G2/F2]s [G/F3]s [G2/F2]s [G/F3]s [F/G3]s [F/G3]s [F2/G2]s [F2/G2]s [F3/G]s [F3/G]s
For the hybrid with flax external layers, the effect of the dissipated energy by the flax layers predominates (Figure 4.b). Thus, the observed decrease of the loss factors as a function of frequency of the flax fibre composites leads to a decrease of the loss factors of the hybrid with flax external layers.

In order to compare the loss factors of all studied hybrids, the evolution of these latter is plotted according to flax fibre volume fraction for a 500Hz frequency (Figure 5).

Discussion

The choice of a composite material in a structure needs to be proved by a good mechanical performance and damping properties. Tables 2 and3 show the properties of all hybrid composites of the study. These results show that replacing flax internal layers by glass layers leads to a lowering of both specific bending moduli and loss factors of the composites (table 2). Otherwise, the stratified [G 2 /F 2 ] S , made of internel layers of flax, has very interessant mechanical and dynamicals properties. In fact, this material has a specific bending modulus and loss factor which exceed those of the glass fibre composite by about 31% and 39% respectively (table 3).

Conclusion :

Compared to glass fibre and carbon fibre composites, flax fibre showed a very interesting damping properties. Those composites have also a specific bending modulus similar to glass fibre composites one. This flax-glass hybridation has been proposed to have a material with better combined properties.

The results of this study proved that hybrid materials have better dynamic properties with, in some cases, better specific moduli. Results also showed that hybridation can provide better properties with an important mass gain.

Figure 1 .

 1 Figure 1. Multilayred finite element S4.

  Loss factors are evaluated in the axis   , LT of the material in each layer: 11 L and T direction respectively, and 66e k is the loss factor in the plan   , LT .

Figure 2 .

 2 Figure 2. Evaluation of specific modulus as a function of the volume fraction of flax fibre of hybrid materials. Volumique fraction of flax fibre (%) 0 5 10 15 20 25 30 35 Specific bending modulus (GPa.cm 3 /g)

Figure 3 .Figure 4 .

 34 Figure 3. Loss factor as a function of frequency (experimental and modelisation results)Figure 4. Concerning the [G 3 /F] S laminate where the external layers are glass layers, the loss factor slightly increases with the frequency. In fact, the stored energy in the laminates is essentially dissipated by the glass layers (Figure 4.a), that's why the evolution of the loss factor as a function of the frequency is directed by that of the glass fibre.

Figure 4 .

 4 Figure 4. Dissipated energy in flax layers and glass layers for a 500 Hz frequency as a function of volume fraction of flax fibre in hybrid with (a) internal flax layers and (b) external flax layers.

Frequency

Figure 5 .

 5 Figure 5. Loss factor as a function of flax fibre volume fraction of hybrid composites for 500Hz frequency A significant increase of loss factor is observed especially when the external layers are flax layers. For example, the loss factors of the [G 2 /F 2 ] S and [F 2 /G 2 ] S laminates are 40% and 146% greater than glass fibre composite, respectively. This growth of damping properties is due to the presence of flax layers which have a better loss factor than glass ones. The loss factor induced by the flax layers in the [G 2 /F 2 ] S laminate only represents 30% of the global loss factor, whereas it represents 96% in the hybrid [F 2 /G 2 ] S . This difference between the [G n /F m ] S and [F m /G n ] S stacking sequences is essentially attributed to the repartition of the stored energies in the different layers. In fact, according to the chosen stacking sequence, between 19% and 88% of the stored energy is dissipated by glass layers in the case of [G n /F m ] S hybrids. In the other hand, in the case of [F m /G n ] S , the energy dissipation is directed by flax layers which represent about 85% and 98% of the global energy.

Flax

  , M. Asssarrar, D. Scida, R. Ayad, X-L. Gong

Table 1 .

 1 bending modili of flax-glass hybrid composites

	composite	Flax	[F 3 /G]	[F 2 /G 2 ]	[F/G 3 ] Glass	[G 3 /F] [G 2 /F 2 ] [G/F 3 ]
	E fx (GPa)	21,38	22,08	22,75	25,11	43,21	42,64	41,63	33,19

Table 2 .

 2 Mechanical and damping properties of hybrids with flax external layers compared to flax composite

	Composite	Flax	[F 3 /G] s	[F 2 /G 2 ] s	[F/G 3 ] s
	E fx (GPa)	21,39	+2,95 %	+6,36 %	+17,40 %
	E fx /ρ	18,01	-5,56 %	-11,11 %	-16,67 %
	η (%)	2,01	-5,01 %	-9,01 %	-25,01 %

Table 3 .

 3 Mechanical and damping properties of hybrids with flax internal layers compared to glass composite

	Composite	Glass	[G 3 /F] s	[G 2 /F 2 ] s	[G/F 3 ] s
	E fx (GPa)	43,37	-1,68 %	-4,01 %	-23,46 %
	E fx /ρ	22,01	+18,18 %	+31,82 %	+13,64 %
	η (%)	0,74	+5,41 %	+39,19 %	+85,13 %