NUMERICAL ANALYSIS OF A REYNOLDS STRESS MODEL FOR TURBULENT MIXING: THE ONE-DIMENSIONAL CASE

A mixed hyperbolic-parabolic, non conservative, Reynolds Stress Model (RSM), is studied. It is based on an underlying set of Langevin equations, and allows to describe turbulent mixing, including transient demixing effects as well as incomplete mixing. Its mathematical structure is analysed, and specific regimes, related to acoustic-like, Riemann-type, or self-similar solutions, are identified. A second-order accurate numerical scheme is proposed in arbitrary curvilinear geometry. Its accuracy and convergence behaviour are tested by comparison with analytical solutions in the different regimes. The numerical scheme can be generalized to multi-dimensional configurations, with potentially cylindrical symmetry, on unstructured meshes. Mathematics Subject Classification. 65M12, 65M22, 76F25. Received May 9, 2020. Accepted July 13, 2021.


Introduction
Turbulent mixing at fluid interfaces plays an important role in a wide variety of domains, ranging from the study of astrophysical objects like supernovae, to engineering applications like Inertial Confinement Fusion (ICF). Depending on the configuration, turbulent mixing can be triggered by different instabilities, among which the Rayleigh-Taylor and Richtmyer-Meshkov instabilities. Reynolds Stress Models (RSMs), presented in Section 2, are common tools to predict the main effects and evolution of such fully developed turbulence. RSMs are one-point statistical models, which rely on a local-in-space decomposition of instantaneous fields into a mean and a turbulent fluctuating field. Still, they do not compute the fluctuating fields, but rather compute their double correlations. Among existing RSMs, the BHR [2] and GSG [19,20] models are widely used for engineering purposes related to variable-density turbulence and mixing at interfaces. They rely on a set of evolution equations for the covariances of velocity and specific volume (or other equivalent covariances, in the case of the BHR model). However, the mass fraction fluxes, which are defined as the covariances between velocity and mass fractions, are algebraically closed, and expressed in terms of the evolved variables of the model. Hence a first gradient closure is used, which implies that the material mass fractions evolve according to a nonlinear parabolic equation. This modelling is firmly established in the self-similar regime, assuming turbulent spectra at equilibrium, but it lacks some important phenomenology, as mentioned in [37,38]. More precisely, countergradient turbulent transport, with transient demixing effects, and incomplete mixing, with segregation effects, are not described. Yet, they are likely to appear in many configurations of interest, for instance just after a shock wave crosses a Turbulent Mixing Zone (TMZ) from the heavy-fluid side to the light-fluid side [8]. In order to take into account these effects, the closure on mass fraction fluxes should be postponed at a higher level by adding evolution equations for the mass fraction fluxes and covariances. This yields a coherent second-order model (in a sense that will be precised in Sect. 2).
In doing so, the nature of the model changes. The conservative parabolic equation for the material mass fractions is superseded by a non-conservative, mixed hyperbolic-parabolic subsystem. Physically relevant solutions rely on both hyperbolic and parabolic parts of the model, with comparable magnitude (see the analytical solutions of Sect. 3.2.2), at least when the characteristic time of the external forcing by the mean flow is not too small compared to the turbulent one. However, the regularizing closure may become unphysical when strong accelerations occur, for instance when a shock wave crosses a TMZ. In that case, the parabolic terms imply a supersonic transport of information from the shocked flow to the flow upstream from the shock front. This violation of a basic physical principle leads to a "modelling" instability as shown in [36]. Such failures of the closures must be cured through limitations which reduce or even discard the regularizing terms.
The purpose of this work is to propose a numerical scheme that can be applied to the hyperbolic subsystem, appart from the regularizing terms. The proposed scheme falls in the class of Godunov-type methods, that rely on a solution of the Riemann problem. While they seem appropriate in our case, at first sight (an exact solution of the Riemann problem is known), they should be employed with care and adapted. The proposed adaptation is rather specific to our turbulence RSM, in the sense that some non-conservative products, which are not defined, even in the weak sense, for the Riemann problem, require an appropriate treatment. Other adapations were introduced to enforce discrete invariant domains (realizability), entropy stability and stationary states, in 1D cuvilinear coordinates. These latter accuracy and stability requirements are used by a large class of numerical schemes devoted to non-conservative hyperbolic systems [7,41]. In the context of other turbulence RSM, discrete entropy and realizability have already been discussed [5,6], but from a different point of view (with emphasis on non-conservative products in presence of shocks). Our approximate Godunov scheme is extended to second-order accuracy, and can be generalized to multi-dimensional and non-planar configurations, on unstructured meshes. Little will be said here on the numerical treatment of parabolic operators, having a compact diffusion support, since the turbulent diffusion terms do not exhibit any particular stiffness. This is not the main concern of this work.
The article is organized as follows. In Section 2, the RSM is presented, its origins are detailed. In Section 3, a mathematical analysis of the hyperbolic-parabolic subsystem is proposed. Filiation with other RSMs is pointed out. The mathematical analysis supports a numerical analysis in 1D curvilinear coordinates, in Section 4. A first-order scheme is proposed, which constitutes a basis for a subsequent second-order extension in space and time. A special attention is paid to the simplicity of the algorithm and to stability at each step, having in view the generalization to an unstructured multi-D scheme. In Section 5, we investigate, in a series of numerical tests, the response of the proposed scheme with respect to the variety of possible regimes allowed by the RSM.

Model
The model studied in the present article is a building block of a larger-purpose turbulence model intended to predict turbulent mixing in multi-species compressible flows subject to strong accelerations or decelerations and shock crossing as encountered in ICF applications [24] or shock tube experiments [8,29]. In its previous form [20,21], this model focuses on the effect of pressure and velocity gradients but employs a simple approach for the species mass fractions. Its present extension introduces a more general treatment of concentrations leading to a more complex system of equations. Since the analysis of the full coupled model is a difficult task, this article makes a first step by considering the simplified but physically meaningful situations of freely evolving turbulent mixing zones where the mean pressure and velocity gradients can be discarded.
To be more specific, the model under consideration here pertains to the class of the Reynolds-averaged Navier-Stokes (RANS) turbulence models. RANS models rely on the decomposition of each flow variable into a mean¯and a fluctuating part ′ = −¯with respect to a well-suited Reynolds average 1 , for instance an ensemble average (as defined in [16], Appendix 1) as in the present case. Writing the equations for the mean values¯of the instantaneous variables involves unknown correlations which require to be either modelled from known quantities or evolved from their own equation which anew bring forth correlations of even larger order. This decomposition procedure can be illustrated with the instantaneous equation for the conservation of the mass of a material of index . The latter reads ( )+∇·( U) = 0 with the density, the mass fraction (i.e. mass of material per unit mass of the mixture, made of materials, such that ∑︀ =1 = 1) and U the velocity. When working with conservation equations, it is convenient to use Favre (i.e. mass-weighted) averages together with Reynolds averages. For any quantity , the Favre average is denoted̃︀ and is related to the Reynolds average, , according tõ︀ = / , with the instantaneous density. The Favre fluctuation is then denoted with ′′ = −̃︀. Summing over the indices, then averaging the mass conservation equation leads tō + ∇ · (¯Ũ) = 0 and does not involve any unclosed correlation, when the Favre average of the velocityŨ is used instead of its Reynolds averaged counter-part. On the opposite, the RANS procedure leads to the following unclosed equation for the averaged mass fraction¯˜+¯Ũ · ∇˜+ ∇ · (¯ ′′ u ′′ ) = 0 where the turbulent scalar flux ′′ u ′′ is an unknown quantity. A simple approach consists in writing a closure based on phenomenological considerations ′′ u ′′ = − ∇˜where the so-called turbulent diffusivity ≥ 0 has to be expressed from known correlations. Such a formula, reminiscent of the Fickian approximation for diffusion processes, is referred to as a first gradient closure. A more general approach consists in writing an equation for the flux ′′ u ′′ by combining equations for ′′ and for u ′′ but this implies that the unknown third order correlation ′′ u ′′ ⊗ u ′′ requires a closure. Assumptions and approximations allowing to retrieve the gradient closure of correlations from their unclosed equations can be found in [35], which serves as a guide for defining the turbulent diffusivity and having an idea of the domain of validity of the simple closure. RANS models are called first-order turbulence models when closures are applied to second-order correlations, whereas they are second-order turbulence models when only third-order correlations are closed, while second-order correlations are kept as unknown of the system. The second-order correlation tensor of the fluctuating velocity, u ′′ ⊗ u ′′ , is an important variable of RANS models. It is referred to as the Reynolds stress tensor. Its half trace˜= 1 2 tr( u ′′ ⊗ u ′′ ) is the turbulent kinetic energy. The rate of molecular destruction of˜is described by the variable˜. The ratio˜=˜/˜is a characteristic turbulent frequency. Many models decompose the Reynolds stress tensor into an isotropic and a deviatoric part and use a closure for the latter, thereby needing only one equation for˜, rather than equations for each component of the tensor. On the opposite, models which evolve the full Reynolds stress tensor are called Reynolds stress models (RSM). The present model is of that kind. Some RSMs can be obtained following two distinct but coherent derivation routes. The first route is determined by the procedure described above. It starts from the hydrodynamics equations for the instantaneous fields, introduces a decomposition into a mean and a turbulent fluctuating component, isolates and averages the equations for the products of fluctuations and finally models the unclosed correlations at second order. The second route relies on an underlying model, which directly computes the Probability Density Function (PDF) of the fluctuating fields [31]. This ensures that the tensor built from the double correlations of the fluctuating fields is semi definite positive. Hence, the resulting RSM stands as realizable. It is this second route that is followed to present the RSM under study.
In order to give more insights into the origins of the RSM of interest here, we introduce its underlying PDF model. The RSM shares some important features with RSM equations presented in [39] (see the turbulent flux Eq. (223), terms I to IV), while the underlying PDF model is adapted from [30,32]. It is able to account for a number of important identified mechanisms that occur, for instance, but not only, in turbulence induced by the Richtmyer-Meshkov instability. For the sake of simplicity, the underlying PDF model will be described for this particular type of flows. Thus, before presenting the PDF model, we will first give a few details about Richtmyer-Meshkov turbulent mixing zones. They result from the impulsional energy deposit by a shock wave crossing an interface between two fluids. The unavoidable small scale perturbations of this interface grow in time, proportionally to the Atwood number, = ℎ − ℎ + , which compares the densities, and ℎ , of the lighter and heavier fluids, respectively. The instability ultimately reaches a self-similar, fully turbulent state. Before reaching such a state, compressibility effects progressively fade: the shock gets far from the mixing zone and the Boussinesq approximation becomes valid [37,38,44]. The mean acceleration, mean pressure gradient, ∇ , and mean velocity gradient, ∇̃︀ U, become negligible. The only turbulence production comes from the mean mass fraction gradients, which drive, through the fluctuating fields, both turbulent mass fluxes and a modification of the mixing heterogeneity. This low compressibility limit regime embeds a variety of non-trivial aspects of fully developed turbulence: fluctuating velocity field showing incomplete return to isotropy, transient demixing effects, and incomplete mixing (defined bỹ︂ ′′2 ̸ = 0) between materials [37,38].
To model this type of flows, it has been shown that the Langevin model got satisfactory results [30,32], notwithstanding the issue of anisotropy, that will be left aside, for the sake of simplicity. To deal with mixing, we also propose to use Langevin models, and their corresponding PDF model. Namely, a Langevin model inspired from the Interaction Exchange with the Mean (IEM) model for the instantaneous mass fractions [14], together with a Simplified Langevin Model (SLM) [31] for the fluctuating velocities, with only dissipation and complete redistribution (i.e. complete return to isotropy) terms. The IEM model has been enriched with a Brownian noise, in order to include, by construction, a class of self-similar solutions in the subsequent RSM. Self-similarity is indeed an important, generic feature of a wide class of fully developed turbulent regimes, that cannot be disregarded. Even though it does not allow the respect of the mass fraction boundedness, this model is only meant to be used as a basis for a RSM. For this purpose, this inherent defect of the Langevin model applied to a turbulent field is not prohibitive. This bare Langevin model can be enriched with mean acceleration and compressibility effects, for instance to account for the large scale interaction with a shock [20]. Such an extension will not be considered in this work and we only focus on these reduced Langevin models.
The flow under consideration involves materials, also referred to as "species". For each material ∈ [1, ] composing the fluid, the mass fraction is defined as where is the local mass of material per unit volume, so that = ∑︀ =1 is the fluid density. Each material is itself composed of different constituents among the ones over the whole set of materials. The mass fraction of constituent , in material , is defined as , = , / , where , is the local mass of constituent belonging to the material , per unit volume. As an illustration, ICF capsules [24] involve materials such as the filling gas and solid shells, whereas the deuterium and tritium ions which should fusion are constituents that can be present in different materials. Mixing between materials can be triggered by molecular or turbulent diffusion. The turbulent material mass fraction flux of each material , ′′ u ′′ , is a cornerstone to the turbulent mixing modeling. It can be obtained by integration of the Favre PDF,̃︀, over the fluctuating field q ′′ = (u ′′ , ′′ , ′′ , · · · ):︂ ′′ ′′ = ∫︀ dq ′′ ′′ ′′̃︀ (q ′′ ) = X . The so-called correlation tensor, X, writes The Favre PDF satisfies the following equation The left-hand side gathers terms related to turbulent transport, that arise from the decomposition of the Navier-Stokes equations into the mean and fluctuating components. The procedure that leads to these terms is described in [32], for instance. The right-hand side is a Fokker-Planck model. It is commonly used in turbulence modelling [32], as it can handle with Gaussian PDF with appropriate mean and variance of fluctuating quantities. The coefficients G and H are defined by where the diagonal components of G, that are and , model the fluctuation dissipation -or return to the mean -, while the contribution HH accounts for the dispersion from the mean due to the turbulent motion. I is the identity matrix of dimension . The coefficients that govern the fluctuating velocity dissipation and dispersion are = − 1 2 1̃︀ and 2 = 2 3 ( 1 − 1)̃︀̃︀, respectively. Here, 1 > 1 is a positive constant related to the redistribution processes in the fluctuating velocity components, and determines the isotropization rate of the Reynolds stress tensor. On the other hand, the coefficients that govern the fluctuating material mass fraction dissipation and dispersion are ≡ − 2̃︀ and H , respectively. The off-diagonal coefficient where x is the position of a Lagrangian fluid particle. Let be the corresponding Eulerian PDF, as defined in [30,32]. Then, from [30,32], it can be shown formally that is solution to the Fokker-Planck equation (2.1). Here, and W are independent Brownian noises, and is implicitely defined by the relation ∑︀ =1 = ′′ ′′ , is the Lagrangian time derivative, x the position of a fluid particle. The Langevin equations (2.3), for the Favre velocity fluctuations, u ′′ , is the SLM (12.109) from [32]. The first term in the right-hand side is a turbulent pressure contribution, that arises from the decomposition between the mean and fluctuating components of the Navier-Stokes equations. Equation (2.3) can also be viewed as a simplification of equation (B3) of [20], for the limited scope of the present framework without pressure nor velocity gradients. The Langevin equation (2.2) is a multi-species generalization of the Langevin model (5.52) of [30] (note that Eq. (2.2) is expressed as a Langevin equation for a fluctuating quantity, while Eq. (5.52) of [30] is expressed as a Langevin equation for an instantaneous quantity. The passage from one to another can be made using the equation for the averaged material mass fraction,̃︀ ). It satisfies the mass conservation property ∑︀ =1 ′′ = 0. The two first terms in the right-hand side of equation (2.2) also arise from the decomposition between the mean and fluctuating components of the Navier-Stokes equations.
To many uses, the instantaneous mass fractions are attached to a flow made of well-mixed constituents. In this work, it is assumed that there is no separation effect between constituents inside a material. Thus, the instantaneous mass fractions of the material constituents satisfy the relation This amounts to a so-called "carrying material" approximation. The equations for the second-order correlations can be obtained by integrating the PDF equation (2.1). A compact tensorial formulation of the evolution equations for the second-order correlations can be written as . The triple correlation flux, yet to be closed, is denoted by ℱ . The latter tensorial expression may also be developed component-wise, say for a material with index , having a constituent with index , in presence of another material with index . Also, the equations for the mean material mass fractions and their constituent mass fractions can be deduced from the Favre average of the instantaneous equation ( ) + ∇ · ( U) = 0, and from equation (2.5). This leads to the following system Finally, an equation for the turbulent mean dissipation rate,̃︀, completes the system [15] where 2 > 1 is introduced as a dissipation constant. The last term in the right-hand side of equation (2.12) is a redistribution term among the components of the Reynolds stress tensor, which reflects its isotropization. All other terms in the right-hand side of equations (2.10)-(2.13) are interpreted as dissipation terms, as they carry a minus sign.
In order to complete the modelling at second-order, a closure on the triple correlations needs to be introduced. These are closed with a first gradient model. Note that such modelling relies on a presumed asymptotic regime where external sollicitations can be discarded, which is valid at long time scales. Such physical considerations hold for both first-order and second-order RSMs [34]. The second-order closure reads where > 0 and > 0 are positive diffusion constants related to the transport processes. The model (2.7)-(2.13), supplemented with these second-order closures, may be viewed either as an asymptotic limit, in the so-called Boussinesq regime, with constant-in-time mean density (Eq. (2.7)), or simply as a sub-system of larger purpose RSMs dedicated to time-dependent, variable-density, compressible flows (by considering a general Reynolds averaged Navier-Stokes equations, that embeds the present RSM, without Eq. (2.7)). It is intended to cure several flaws of these RSMs [2,19,20,28]. First, the model becomes fully second-order. This ensures a better consistency between the various fields, since only triple correlations are closed. In particular, the turbulent material mass fraction fluxes, ′′ u ′′ , are now computed by equation (2.10), instead of being closed, in equation (2.8), with the first-order closure as was classically done in [2,19]. Hence, demixing and counter-gradient fluxes (i.e. fluxes that are not oriented according to the direction imposed by the minus sign in Eq. (2.14)), with coherent material mass fraction covariances, are now allowed by the modelling, owing to equations (2.10) and (2.11). Moreover, the first-order model is an asymptotic approximation of the second-order model in the absence of an external sollicitation, and at long time scales, which is detailed in [35]. In this limit regime, the first-order and second-order models share common self-similar solutions (see Sects. 3.2.2 and 3.2.3). The closure modification amounts to introduce terms with a finite relaxation time which affects the transitory behaviour toward self-similar solutions. In situations where short time scales are at stake, this transient modelling is indeed crucial for reproducing physical phenomena under study. Second, the material mass fraction covariances, ′′ ′′ , are also part of the evolved-intime correlations of the model. It is an important data as it can serve as a measure of the heterogeneity of the TMZ [33]. When turbulent combustion modelling comes into consideration [5,39], the model may thus provide with a modification of the combustion rates by the heterogeneity inside the TMZ [33]. Similarly, radiation transfer across a TMZ is modified by the heterogeneity level [11].

Analysis of the continuous 1D model
In this Section, the continuous model is analysed, its important properties highlighted, having in view its discretization in the remainder of this article.

The hyperbolic sub-system
In this section, the closure with diffusion terms, having a superscript , are dropped in equation (2.6) and system (2.10)-(2.13).

Realizability
Proposition 3.1. Let ℱ = 0 in equation (2.6). Consider X( ) a solution of this system, and assume that︀ U, X, G, H are smooth enough so that all terms have classical meaning. Then, if X(0) ≥ 0, in the sense of symmetric matrices, X( ) ≥ 0, ∀ > 0.
Proof. We first consider the casẽ︀ U = 0. Then, defining one easily checks that where I is the identity matrix. Therefore, the solution X of (2.6) reads In the above sum, both terms are semi-definite positive, hence X(x, ) ≥ 0. Next, we turn to the case wherẽ︀ U ̸ = 0. Equation (3.1) is then valid if one follows the characteristics of the field̃︀ U. More precisely, let us define = (x, ) the solution to )︂ .
Then we have Therefore, we deduce that X( (x, ), ) ≥ 0. In order to infer that X(x, ) ≥ 0 for all x, we only need to check that any x satisfies x = (x 0 , ), for some x 0 . In order to prove this, we consider the backward characteristics Sincẽ︀ is smooth enough, the map ↦ → (x, ) is defined for all . Setting x 0 = (x, − ) ensures that x = (x 0 , ), which concludes the proof.
Remark 3.2. The smoothness assumption in Proposition 3.1 is important for two things: first, smoothness on X, G, H allows to write down the equation in a strong form and use Duhamel formula (3.1). Second, smoothness of̃︀ allows to follow the characteristics backward, since it exists for all time. We do not know how to prove such a result in a non-smooth case. However, as we will see below (Prop. 3.6), the particular solution of the Riemann problem we use to build the scheme does satisfy the realizability property.
This multi-species system is hyperbolic if and only if the mono-species, quasi-linear system (3.16), is hyperbolic, and all its characteristic fields are linearly degenerated (LD). Its eigenvalues arẽ︀ U,̃︀ U + ± , and̃︀ U + , while the corresponding set of Riemann invariants, , ± , and , can be expressed as Proof. We find it useful to write the 1D system in mass coordinates, using the variable change d = d , where is one cartesian coordinate, or the radial coordinate. Then the system writes We rewrite the system (3.3)-(3.8) under the following compact, quasi-linear formulation from which we deduce straightforwardly the eigenvalues of the system (3.3)-(3.8):̃︀ U,̃︀ U + ± , and̃︀ U + .
Here, the subscript ⊥∈ { , } refers to any transverse component of the mass fraction flux.
Since the model is realizable, all the eigenvalues are real-valued. The model reformulation (3.9)-(3.14) makes it also clear that hyperbolicity can be studied for each material independently, discarding equation (3.11). The latter can indeed be viewed as a diagnostic (one-way coupled) equation, whose eigenvalue is̃︀ U. We can also isolate the components of the Reynolds tensor, in equation (3.13), that are independent from the other variables of the quasi-linear system (the eigenvalue associated to these components is alsõ︀ U). The equation for the remaining component, + = − − , can be written The dissipation rate of the turbulent kinetic energy,̃︀, obtained from equation (3.14) of the quasi-linear system, is independent from the other variables of this quasi-linear system as well. The associated eigenvalue is alsõ︀ U. Therefore, hyperbolicity of the model (3.3)-(3.8) is equivalent to the hyperbolicity of the reduced, mono-species system of quasi-linear equations (3.9)-(3.10)-(3.12)-(3.15). The latter can be rewritten as where I is the identity matrix whose size is equal to the number of constituents in the material . The eigenvalues of the matrix arẽ︀ U,̃︀ U + ± , and̃︀ U + , with respective multiplicity equal to 2, 1, and to the number of material constituents in the material . Their right eigenvectors are for the eigenvaluẽ︀ U, respectively, for the eigenvalues̃︀ U + + and̃︀ U + − , for the eigenvaluẽ︀ U + , with e = ( 1 , 2 , · · · ) , using the Kronecker symbol . The determinant of the matrix composed by these right eigenvectors (the change-of-frame matrix to the eigenspace) can be computed, and is equal to + − . Therefore, this determinant is positive provided̃︂ ′′2 is non-null, that is inside the TMZ. This completes the proof of hyperbolicity.
Finally, the field associated to the eigenvaluẽ︀ U + is LD because it is a function of all field but̃︀ , /̃︀ , and due to the structure (non-null components) of the corresponding eigenvector. Moreover, the field associated to the eigenvaluẽ︀ U + ± is LD since ∇ (̃︀ U + ± ) = (0, 0, 0, 1, 0) is orthogonal to the eigenvector of̃︀ U + ± . Therefore all the fields of the system are LD.

Entropies
for each material, denoted by the subscript . The entropy expressions read Both share the same entropy flux, Proof. The Reynolds tensor, u ′′ ⊗ u ′′ , is nonsingular inside the TMZ. Moreover, its derivatives can be obtained with the relation and its derivatives are well-defined quantities. The next step is to prove that both functions and are convex functions. Concerning , it is indeed convex, because it is the sum of two convex functions. As such, and also because is non-positive, it is an entropy for the system (3.3)-(3.8). Next, the convexity of can be deduced from the convexity of the following function: where ( ) is the positive definite cone for matrices of size . The proof starts by noticing that the following function is convex because its Hessian matrix has only non-negative eigenvalues. Then, the spectral theorem can be used to show that, given , ∈ ( ), there exists a matrix which satisfies = and = , where = diag( 1 , · · · , ) is a diagonal matrix. Let x,y ∈ R . We obtain . Introducing ∈ (0, 1), we finally get which has been obtained from the convexity of the function applied on each index . Hence, the function is convex, and so is .
is also a couple of entropy and entropy flux.  .7), with null right-hand side, in mass coordinates. This system admits a four-state piecewise constant solution of a Riemann problem in the mass coordinates, that is:

21)
where the intermediated states, are defined as the unique solution of the algebraic system , and can be either , , * or * .
This solution is realizable, in the sense that if It is also isentropic, in the sense that it is solution of system  ]︁ > 0, where 0 is the position that initially makes the separation between the and states. This can be shown classically by first integrating equation (3.32) by parts, then inserting the solution of System (3.22)-(3.31) and integrating the derivatives of the test function. Such procedure cannot be applied a priori for equation (3.5), because it has a non-conservative product. However, this difficulty can be circumvented, noticing that equation (3.5) with null right-hand side is satisfied in a strong sense, at = 0 , and for all times, by the solution of System (3.22)-(3.31). Therefore, it is sufficient to show that equation (3.5) is satisfied in a weak sense on [0, 0 [ and ] 0 , ∞[ separately, which writes ]︁̃︀ ]︁ ( , ) = 0, (3.34) Figure 1. Riemann fan associated to the system (3.9)-(3.14), in a 1D geometry, with transverse symmetry in cartesian coordinates, assuming̃︀ U = 0.
for any test function smooth enough. This is true because the integration by parts is allowed, even with the presence of a non-conservative product, because As a prelude to the discretization step, several remarks can be done.
Remark 3.7. The particular solution of the Riemann problem that is exhibited in Proposition 3.6 can be built using the Linearly Degenerated (LD) stucture of the system of interest. The continuity of Riemann invariants is expressed across the waves having speed̃︀ U+ ± and̃︀ U, as illustrated in Figure 1. Even in this LD case, we cannot rely, to our knowledge, on a well-established theoretical result, that would grant a unique solution of the Riemann problem (for such a non-strictly hyperbolic, non-conservative system). However, it has been shown a posteriori that the solution of (3.22)-(3.31) is indeed a solution of the Riemann problem. Conversely, every four-state piecewise constant solution of the Riemann problem must satisfy the Riemann invariant continuity relations (3.22)-(3.31).
Remark 3.8. The algebraic system (3.22)-(3.31) can be solved in two steps. The first steps amounts to solve System (3.24)-(3.29), that can be put in terms of a two-equation system and interpreted as an "acoustic" solver. The second step can be obtained by inserting the first step solution in equations (3.30) and (3.31). This decomposition reflects that equation (2.11) is a diagnostic (one-way coupled) equation of the system of interest.
Remark 3.9. The Riemann invariants have been expressed as functions of the stationary states.

Analysis of the non-conservative products
Based on a particular solution of the Riemann problem, made explicit in Proposition 3.6, and taking advantage of the LD nature of the system, we can propose an analysis of the non-conservative products. First, from the following proposition, either u ′′ ⊗ u ′′ or̃︀ is smooth across each wave of the Riemann problem.
Owing to Proposition 3.10, the non-conservative product, u ′′ ⊗ u ′′ ∇̃︀ , in equation (2.10), is always defined in a weak sense for the Riemann problem. Therefore, it cannot be considered as a pathological non-conservative product. A similar argument does not hold for the other non-conservative product, ′′ u ′′ · ∇̃︀ , in equation (2.11), owing to Remark 3.11.
Remark 3.11. Both variables ′′ u ′′ and̃︀ may be discontinuous fields across the extremal waves of the Riemann problem.
This apparent pathology can however be adressed, relying on the following proposition.
The left-hand side of equation (3.37) is defined in a weak sense, provided the temporal derivative is kept Lagrangian. The hyperbolic system has been analysed through the prism of the Riemann problem, and the identified pathological non-conservative products can be recast under a well-posed formulation. The system can be viewed as "marginally non-conservative".
Remark 3.13. We do not have much theoretical result at disposal for discontinuous solutions: uniqueness of the solution, as well as realizability, have not been proven. Our analysis of non-conservative products is based on a particular realizable solution of the Riemann problem. Therefore, our approach will consists in encoding as much as possible the known theoretical results in the numerical scheme, at the discretization step, using a Godunov approach.

The full hyperbolic-parabolic system
In this section, we consider the additional contribution of the diffusion-like closures that were not considered in Section 3.1.

A class of self-similar solutions, with the diffusion contribution
Let us assume a fixed background fluid having zero mean velocity, in 1D geometry with transverse symmetry in cartesian coordinates, homogeneous turbulent frequency, homogeneous density, and isotropic Reynolds tensor across the TMZ. Under these asumptions, System (2.8)-(2.13) can be simplified. The evolution of the turbulent kinetic energy,̃︀, and its dissipation rate,̃︀, reduce to the following̃︀ −̃︀ system [15] . The remaining equations of our model becomẽ︀ The existence of self-similar solutions to the problem of freely developing TMZ is suggested by dimensional analysis, in the absence of external dimensioning parameters. The TMZ width Λ( ) itself can then be taken as the scale length of the problem. It means that possible self-similar solutions may have the form̃︀ =ˇ( /Λ( )), where is an integer andˇis an arbitrary function. Both the time exponent and the shape functionˇdepend on the considered quantity, . In addition, one can assume TMZ having compact support, meaning zero-valued turbulent quantities outside the TMZ of increasing width. Under these considerations, simple profiles can be assumed, namely parabolic and linear profiles, for the turbulent quantities and mass fractions, respectively. Following closely the computation of [10], we obtain some self-similar solutions for this system︀ The following compatibility relations must be satisfied for this self-similar solution to hold Finally, all the above discussion may be summarized in the following result:    [3,10]. The stability of these solutions, of great interest, has not been studied here, as it is beyond the scope of this work. See [18], for instance, for examples of such stability analysis led on simpler RSMs. Remark 3.19. Proposition 3.18 shows that these self-similar solutions establish a connection between the second-order RSM (2.8)-(2.13) and the GSG and BHR RSMs [2,19]. The latter indeed rely on equation (3.46), for which the material mass fraction flux, ′′ u ′′ , is closed at first-order instead of being computed.

Numerical analysis in 1D curvilinear coordinates
We present, in this section, a first-order numerical scheme that maintains, at the discrete level, the main properties of the model (3.3)-(3.8), namely mass conservation, realizability, preservation of stationary states, and non-increasing entropies. This scheme is also able to handle correctly the non-conservative products, because it relies on a reformulation of the hyperbolic model (3.3)-(3.8) with equation (3.37). The proposed scheme can be viewed as an approximate Godunov scheme. As a preambular statement, let us notice that schemes having collocated components of the correlation tensor are good potential candidates to obtain the discrete realizability. Moreover, locally non-increasing discrete entropies can be obtained from collocated scheme, see [9,13,25] for Euler hydrodynamics.

Preliminaries: notations at the discrete level
Let us consider a space domain Ω = ∪ ∈ Ω , paved with cells. The volume of the th cell Ω of a given set of cells , is denoted by |Ω |. Given a 1D geometry along the unit vector e , each cell Ω is an interval is the set of node indexes that are located at the boundaries of Ω , while is the set of cells sharing the node with index .

A splitting approach
We propose a splitting scheme to solve equations (3.3)-(3.8), during ∆ = ( +1) − ( ) . The left-hand side of this system is first evolved from an initial time, labelled with ( ), to an intermediate diamond state, labelled with ◇, during the time step ∆ . The remaining contributions, namely the dissipation and redistribution contributions, are then evolved from the diamond state, to a final state, labelled with ( + 1), during the time step ∆ . Note that stationary states do not prevent a priori this splitting approach, as they do not involve a balance between terms located from both sides of equations (3.3)-(3.8).
The discretization leading to the diamond state is described from Sections 4.3 to 4.10, while the discretization leading to the dissipation and redistribution contributions is described in Section 4.11.

An approximate Godunov scheme for the mass fractions and their fluxes
The mass fractions and their fluxes evolve according to equations (2.8)-(2.10). The resulting sub-system stands as a non-conservative sub-system of our model, whose non-conservative product, u ′′ ⊗ u ′′ ∇̃︀ , can be analysed with the help of the Riemann invariants. Therefore, a Godunov-type scheme is a good candidate for its discretization. Having the knowledge of the stationary states of the model, we propose to introduce such a scheme for the variables = (︁̃︀ , u ′′ )︁ , which are exactly the variables that can become stationary. We will show later on that this choice does lead to the discrete preservation of stationary states. We also discard the dissipation terms in the mass fraction flux equation. We can do so safely because stationary states indeed occur for a null dissipation contribution. These dissipation terms will be reintroduced later on.
Let then ℛ (︀ / , , )︀ be the solution of the Riemann problem in mass coordinates, following [26]. The weak solution of the Cauchy problem, built from the initial piecewise constant data The projection of this solution on piecewise constant function writes with the choice ∆ = |Ω |. Inserting the solution of the Riemann problem and linearizing the characteristic curves as straight lines (required for > 1) yields The following proposition is then obtained from Lemmas 4.1, 4.3 and 4.5.
Finally, a quadrature formula is chosen to relate ′′ ′′ and ′′ ′′ , so as to deduce a discrete equation for

A scheme for the mass fraction covariances that handles with non-conservative products
The discretization of the mass fraction covariance equation relies on the ODE (3.37). Discarding the dissipation terms in the right-hand side, that will be reintroduced later on, we propose the following discretization , (4.10) Remark 4.8. Because both variables, ′′ ′′ and̃︀ , can jump across the extremal waves, the non-conservative product ′′ ′′̃︀ is not straightforwardly defined in a weak sense as such. The mass fraction covariance equation expressed in the primary variables is ambiguous in this respect. The ODE (3.37) restores back a weak sense, and thus controls the non-conservative product, as explained in Section 3.1.6, Remark 4.9. The discretization (4.10) can be reformulated in terms of discrete entropy variations, as (4.11)

A conservative, realizable and entropy scheme
Proposition 4.10. The approximate Godunov scheme (4.2)-(4.9) gives a discrete equation for the material mass fractions that can be recast under its conservative form︀ Proposition 4.11. The approximate Godunov scheme (4.2)-(4.10) satisfies unconditional discrete realizability.
Proof. The discrete equations for the correlation tensor components can be gathered in the tensorial expression Remark 4.12. The discrete tensorial formulation (4.13) can be identified with the continuous tensorial equation (3.2), provided the resolvent matrix, K, is approximated with a first-order Taylor expansion in the time variable: K = I + ∆ G + (∆ 2 ). Our scheme simply brings an additional, unique prescription for a space discretization of the mass fraction gradients, | ( ) , in the turbulent double correlation discrete equations. satisfies the following discrete entropy inequalities 14) where both discrete entropies have the same dissipation rate.
Proof. Relying on equation (4.11), we can only focus on , and establish its discrete entropy inequality. Equation (4.12) can be combined with equations (4.4)-(4.6) (the discrete Riemann invariant relations), to give Next, rewritting component-wise the tensorial equation (4.13), we obtain the discrete equation for ′′ u ′′ . The latter can be rearranged as (4.16) Summing equations (4.15) and (4.16), and noticing that The discrete Riemann invariant relations (4.4)-(4.6) can then be inserted in the second term of the right-hand side, which gives where Enforcing inequality , ≤ 0 yields the stability condition.
Remark 4.14. We observe that the first-order scheme is isentropic at the discrete level, for both entropies, in the very special case of uniform mesh and constant in space Reynolds tensor, with ∆ = ∆ / √︁̃︂ ′′2 ( ) .

A non Godunov-type scheme for the material constituents
At the mixing zone boundary, there is at least one material mass fraction, denoted by index for instance, that vanishes. This may lead to a divergent speed , at the locations where the material mass fraction becomes close to zero. This speed is associated to the material constituents, while the speed associated to the material, ± , do not diverge. The system partially degenerates at such boundary because the wave speed blows up. Thus it is no more hyperbolic, if the material constituent evolution is considered. We can nevertheless mitigate this point, because the Riemann invariant,̃︀ , /̃︀ , associated to the wave speed , has no physical significance (as it is not linked with any physical matter). Its transport at infinite speed cannot be forbidden, contrary to the transport of the constituent mass fractioñ︀ , . Therefore a Godunov scheme relying on the Riemann invariant︀ , /̃︀ is clearly not an appropriate choice, as the CFL condition would force the time step toward zero. We propose an alternative scheme, that can be written as︀ Proof. The discrete equation (4.12) for the material mass fraction,̃︀ = ∑︀ =1̃︀ , , can be obtained from the discrete equation (4.19), by summing over the material constituents.
Remark 4.16. Initial uniform profiles of̃︀ , /̃︀ remain uniform profiles, since this mass fraction ratio is a Riemann invariant. This stationarity property holds true at the discrete level also, owing to the discrete equations (4.12) and (4.19).

Second-order extension in space
We propose the following second-order extension in space. It relies on the reconstruction of a variable, , which can be either̃︀ or ′′ ′′ , at a node ∈ where | refers to an a priori limiter at cells, while refers to an a posteriori limiter at nodes. The a priori limiter is defined as ≥ 0, and | = 0 otherwise. The a posteriori limiter is set to its a priori value, = 1. We shall set = 0 only if an a posteriori entropy criterion, yet to be defined, would be violated. This reconstruction is then inserted in the expressions for the discrete Riemann invariant at nodes which can be viewed as a nodal solver, whose solution is These expressions can be inserted in the discrete equations (4.12), (4.13) and (4.19), that rely on a simple forward Euler time discretization. It could also be combined with a Runge-Kutta higher-order extension, as we shall show in the next Sections.
Proposition 4.17. The second-order a priori reconstruction preserves linear solutions, and degenerates at first-order at extrema, by construction.
Proposition 4.18. The second-order scheme preserves the stationary states and the mass conservation relations (3.38), at the discrete level, by construction.
In the following, we investigate the stability of this high-order in space extension, in relation with the time discretization.

Runge-Kutta 2 stabilization: a Von Neumann analysis
The extension of the space discretization to second-order may affect stability. To evaluate its impact, a standard Von Neumann analysis is proposed for the discrete equations satisfied by the variable = (︁̃︀ , u ′′ )︁ . Strong hypotheses are made to do so. Hence, all subsequent stability constraints should be understood as conditions related to linear stability, distinct from nonlinear stability.
and is set to 0 or 1, respectively for the first-order or the second-order case. The spectral radius of is bounded by one, under the CFL criterion where -the first-order discretization in space ( = 0) yields: sin 2 -the second-order discretization in space ( = 1) yields: The scheme will be deemed as linearly stable, with respect to the exact harmonic solution, if all eigenvalues of are, in modulus, less than or equal to one. From Proposition 4.21, we observe that the constant mode, = 0, is unconditionally stable for both first-order and second-order discretization in space. Note also that the smallest wavelength, with wavenumber = /∆ , is stable under the condition ≤ 1 in either cases. The first-order scheme in space and time is 2 stable if ≤ 1. On the contrary, the second-order in space, first-order in time scheme is unconditionally unstable. Figure 2 indeed shows the modulus of the eigenvalues (amplification factor) related to the first wavenumbers, against the Courant number. This suggests the existence, for each Courant number, of at least one wavenumber having an amplification factor greater than one.
Let us now investigate whether a higher-order time discretization might bring stabilization. Note that the Taylor series should be pursued for higher-order time discretization. At the third and fourth order, this series corresponds to the amplification matrices of the third and fourth order, low storage Runge-Kutta scheme of Appendix A. Figure 3 shows the amplification factor of the first wavenumber, for the high-order time discretizations, up to fourth-order, against the Courant number. We observe that the stability of our scheme only depends on the stability of the 2∆ -wave, because the amplification factor from this particular wave dominates the other   Figure 2, now related to the high order time discretization: second-order (left), third-order (center) and fourth-order (right). Increasing the wavelength now broadens the stability domains.
ones, as soon as one wave becomes unstable (its amplification factor exeeds unity). The CFL criteria are given analytically in Table 1 for the first and higher-order scheme.
As a conclusion, we have shown that a stabilization is required for the second-order in space discretization. If combined with a first-order forward Euler time discretization, it is indeed unconditionally unstable. We have also shown that high-order Runge-Kutta schemes provide the desired stabilizing effect, and that a second-order Runge-Kutta time discretization is sufficient. In the next section, we shall show how the second-order spatial discretization can be appropriately combined with such second-order Runge-Kutta time discretization, having in view an efficient and simple implementation.
is Semi Definite Positive (SDP), then is also SDP.
The proof of Lemma 4.23 is postponed to Appendix B, and can be easily extended to the case where is a matrix and a vector, using the Sylvester characterization for the SDP property. Proof of Proposition 4.24. Dropping the cell index, the discrete correlation tensor can be written as The first term of the previous equation for X ◇ is clearly Semi Definite Positive (SDP). Therefore, requiring the second term to be SDP is a sufficient condition for X ◇ to be also SDP. This is the case because the matrix ∆ = (∆ ) is also SDP, according to Lemma 4.23.
In the following section, we rely on Proposition 4.24 to propose a second-order Runge-Kutta stabilization on only, followed by the discrete equation for the mass fraction covariances (4.10). This minimizes the required storage, while ensuring discrete realizability.

An entropy-based a posteriori treatment
We wish to obtain a second-order in space scheme for smooth solutions, while ensuring a nonlinear entropy stability. Recall that linear stability has already been obtained with a second-order in time Runge-Kutta scheme, which controls smooth solutions of the linearized hyperbolic system. Note that such entropy stabilization may also bring an additional control over spurious oscillations, for non-smooth solutions. To achieve so, let us investigate the modification of the first-order discrete entropy inequality (4.17), by the second-order scheme. Noticing that the discrete entropy dissipation rates are the same for both entropies and , at first and also second-order, we drop safely and focus on , which is related to the variable . The explicit Euler stages, = 1, 2, of our RK2 method, can be written where ℒ and ℒ are symbolic operators refering to the first-order in time, second-order in space discretizations of Sections (4.3), (4.6) and (4.7). The final convex combination stage of our RK2 method writes, ∀ ∈ , .

(4.28)
This RK2 method brings the following discrete entropy increments where the superscript (*, ) denotes the second-order solution of the nodal solver at the Runge-Kutta stage having index , and are the (unsigned) second-order entropy dissipation rates. The last Runge-Kutta stage to the diamond state is a convex combination stage. Proof.
by convexity of the function .
From Proposition 4.25, we obtain straightforwardly . Gathering all Runge-Kutta stages, and defining we finally obtain the second-order discrete entropy inequality This latter expression gives a non-ambiguous discrete expression of the entropy flux through the Runge-Kutta stages. Now equipped with the first-order and second-order discrete entropy inequalities, respectively (4.17) and (4.29), we can step toward a discrete entropy-based a posteriori treatment for the high-order approximation of . Relying on the isentropic nature of the Riemann problem in the absence of diffusion and dissipation, we impose, a posteriori, the left-hand side of the second-order entropy inequality (4.29), ℰ 2nd , ℎ , to be bounded by the first-order discrete entropy dissipation rate, defined by equation (4.18), and denoted by (1st) , . Would this criteria be violated, then the scheme would locally reduce to first-order. The order degradation is performed following the e-MOOD approach [4], where the local reconstruction order is determined uniquely at each node, for the sake of coherence in the nodal solver, which gets the same reconstruction order from all the neighbour meshes. The reconstruction is therefore not unique for cells which have both "first-order" and "higher-order nodes". Let us define -and the set of cell indices and nodes indices, respectively, )︁ a boolean that locally points out the need for an order degradation, -= { ∈ / or +1 or −1 }, the domain that requires degradation to first-order, -RK2 the set of 3 stages of our second-order Runge-Kutta scheme. We have chosen a specific Runge-Kutta scheme, see Appendix A and references [17,43], so that each stage, having index , and denoted by RK2( ), can be either an "Explicit Euler" stage, or a "Convex Combination" stage.
The algorithm that makes the transition from the initial state, , to the diamond state, The second step of the splitting computes a final state, from the intermediate diamond state, owing to the discrete equations (4.30)- (4.33).
The numerical solution satisfies, for arbitrary geometry, the following properties, at the discrete level: -entropy stability, -preservation of stationary states, -realizability, -conservativity of the mass fraction and material constituent equations, -mass conservation.
Moreover, the stability restriction on the time step is not stringent and ensures that the time step does not vanish.

Numerical tests
We propose four different test cases, three of them having analytical solutions. They highlight the variety of possible regimes that the RSM under study is able to describe. The first one focuses on wave-like behaviour in an acoustic-like regime, while the second focuses on transitory discontinuous solutions, occuring typically at the initialization stage of the model. The third test case deals with the important self-similar regime. The fourth test case deals with the transient demixing regime. In the first three cases, we investigate the convergence rate, showing the benefit of the second-order scheme. For the sake of clarity, a mixture of only two materials, having index and , will be considered in each test. The mass fractions and fluxes of material is easily deduced from that of material , hence we present results only on material . The mean velocitỹ︀ U is assumed to be zero.

Test case in the acoustic-like regime
In the regime considered here, the settings are compatible with the Von Neumann analysis of Section 4.8, with periodic boundary conditions. The stability analysis is confirmed numerically for the variety of wavenumbers, and time discretizations shown in Figures 2 and 3. Let us now investigate the convergence properties and the behaviour of the a posteriori procedure. We recall that the diffusion and dissipation terms are discarded. The density is set to one, and the Reynolds tensor is chosen as a constant, in a planar geometry. Thus the "acoustic" speed reduces to = √︀ ′′2 . The initial mass fraction fluxes and covariances are set to zero, while the initial mass fraction is chosen as̃︀ Here, the spatial domain, Ω = [−10; 10], is chosen periodic, with length |Ω| = 20. It is splitted uniformly in cells, so that |Ω | ≡ ∆ . The wavenumber is set to = 20. The analytical solution reads︀ For the purpose of convergence study, the spatial resolution varies in the range ∆ ∈ [10 −4 ; |Ω|], while the time step is fixed at ∆ = 10 −4 . Hence, the stability condition is satisfied for any value of ∆ considered. We show in Figure 4, on the right-hand side, the convergence curves. The logarithm of the 1 -error is plotted against log , for each variables, and for both the first-order (solid curves) and second-order (dashed curves) schemes. The convergence rates are given in Table 2. As expected, the second-order scheme is second-order  convergent. The convergence rate of the first-order scheme is slightly smaller for the mass fraction covariances. This is due to the fact that these variables depend nonlinearly on the numerical dissipation of both the mass fractions and their fluxes.
As an illustration, Figure 4 also shows, on the left-hand panel, the numerical and analytical solutions at time = 2.5 with = 300 and ∆ = 10 −4 . We clearly observe that the second-order scheme is closer to the analytical solution, because it is less damped than the first-order scheme. For larger time steps, with the same spatial resolution, the improvement from the first-order to the second-order scheme can still be observed at = 1/2, where = √ ′′2 ∆ /∆ . Above = 0.7, the second-order solution is degraded to first order, by the a posteriori entropy constraint, which we consider acceptable.

Initialization test case: the Riemann problem
Let us now investigate the behaviour of our scheme with respect to the discontinuous, analytical solutions of a Riemann problem. We wish, in particular, to evaluate the treatment of the non-conservative terms by our scheme. Diffusion and dissipation terms are again discarded, while the initial mass fractions are chosen as Heaviside functions. An initial constant heterogeneity level is prescribed via the mass fraction covariances. The analytical solution is resumed in Table 4, where the initial condition is chosen realizable. A mixing zone evolves from the initial location of the mass fraction discontinuity, and expands at a constant speed, equal to the diagonal elements of the Reynolds tensor. For each material, two constituents have been considered, having initial linear profiles. The constituent ratio was initially 99/1% at the domain edges, for the material with index . Table 3. Convergence rates corresponding to Figure 5, for the Riemann problem. The constituent ratio for the material with index was initially chosen 70/30% at the domain edges. The resolution domain, Ω = [−10; 10], is splitted uniformly in cells, so that |Ω | ≡ ∆ . The spatial resolution varies in the range ∆ ∈ [0.02; |Ω|], while the time step is uniquely set at ∆ = 10 −3 , to satisfy the stability conditions at the highest spatial resolution, that corresponds to ∆ = 0.02. Figure 5 (with Table 3), in the right-hand panel, and Figure 6, show good convergence results. The logarithm of the 1 -error is plotted against log , for each variable, for both the first-order (solid curves) and secondorder (dashed curves) scheme. As expected, see [12,27], the first-order scheme converges with a 1/2 rate. The high-order scheme proves to converge at a 2/3 rate. These convergence rates are maintained for the material constituents, despite the fact that the high-order correction of the scheme only occurs in the nodal solver.
On the left-hand panel of Figure 5, the numerical and analytical solutions are plotted at = 10 with a resolution of = 500 cells, with ∆ = 10 −3 . We do not observe any pathology due to the non-conservative product, which would occur at the initial discontinuity position for a standard Godunov scheme applied for the mass fraction covariances. We again observe the improvement of the second-order over the first-order scheme.
Note that the second-order to first-order degradation by the a posteriori loop is activated at the discontinuities only, whatever the space and time resolution.

Test case in the self-similar regime
We finally propose a test case where diffusion and dissipation terms are considered. The dissipation and diffusion constants are given in Table 5. We make use of the analytical self-similar solution of Section 3.2.2, with̃︀ = 1 and̃︀ = 1, to obtain a reference solution.
At this point, let us give a few comments about the discretization of the tensorial diffusion operator, which can be written component-wise, under the generic form (in the present conditions, assuming constant density) where is the turbulent diffusivity. It is splitted from the hyperbolic part, with a Strang splitting, to get a global second-order time discretization if a second-order global accuracy is required. A centered discretization in space is used. The semi-discrete scheme is  Table 4 (dotted curves), and numerical solution from the first-order (solid curves) and second-order (dashed curves) schemes. The results for the mass fractions (blue), mass fraction fluxes (red) and mass fraction covariances (green) are plotted at = 5 for = 500 cells. Right-hand panel: logarithm of the 1 -error against log . The dotted lines correspond to the observed convergence rates, 1/2 and 2/3 ( Table 3). The symbols correspond to the different mesh refinements. The associated best-fit lines are also shown. Figure 6. Riemann problem: convergence study for the material constituents. Logarithm of the 1 -error against log . The dotted lines correspond to a 1/2 and 2/3 slope, close to the observed convergence rates 0.497 and 0.65, for the first-order and the second-order schemes, respectively. The symbols correspond to the different mesh refinements. The associated best-fit lines are also shown.
with the arithmetic average | +1/2 = (∆ +1 | +∆ | +1 )/(∆ +∆ +1 ), and the definition ∆ +1/2 = (∆ + ∆ +1 )/2. Using this arithmetic average makes the TMZ support (i.e. the domain of non-zero Reynolds tensor) able to expand, which is not the case if one uses harmonic averaging. The time discretization of the variable is chosen implicit, while the discretization of the coefficients are considered explicit. The resolution domain, Ω = [−10; 10], is splitted uniformly in cells, so that |Ω | ≡ ∆ . The spatial resolution varies in the range ∆ ∈ [20/510; |Ω|], while the time step is fixed at ∆ = 10 −4 , for the purpose of a convergence study. Hence, the stability condition is satisfied for any value of ∆ considered. Figure 8 (with Table 6) shows the good behaviour of the scheme with respect to the convergence study. The logarithm of the 1 -error is plotted against log , for each variable, and for both the first-order (with solid lines)  and second-order (with dashed lines) scheme. For the sake of illustration, Figure 7 shows the evolution of the analytical solution of several variables of interest, including the fast computed variables̃︀ , ′′ ′′ , and̃︂ ′′2 , as well as the slow rebuilt variable, ′′ ′′ / √︁̃︂ ′′2̃︂ ′′2 . We observe that the second-order scheme brings a subsequent accuracy improvement for all variables, especially for the slow quantity, see Figure 9 for an illustration. This is remarquable, because the analytical reference remains constant across the expanding TMZ, and abruptly falls to zero at its edge, see Figures 7 and 9. We also observe that the a posteriori, second-order to first-order degradation, is not activated, whatever the spatial resolution, and time step ranging from ∆ = 10 −4 to the maximum stable time step, satisfying = √ ′′2 ∆ /∆ = 1.
We finally investigate the behaviour of the proposed scheme on non-uniform meshes, again in the case of the self-similar solution. We introduce a geometrical progression across the TMZ, so that boundary cells are 100 times larger than the center cell. Figure 10 gives a comparison between the analytical and numerical solutions at = 10, for several variables of interest. As far as mass fractions are concerned, the first-order scheme generates an artificial inflexion point at the TMZ center, while the second-order scheme does not. This illustrates, once again, the beneficial effect of a second-order spatial discretization. The a posteriori correction never activates. We have checked so for time steps ranging from ∆ = 10 −4 to the maximum stable time step, satisfying = 1.  Table 6). In both panels, the dotted lines correspond to the theoretical convergence rates, 1 and 2. The symbols correspond to the different mesh refinements. The associated best-fit lines are also shown. Figure 9. Self-similar test case: the analytical solution (dotted line) and the numerical solutions, namely the first-order (solid curve) and second-order (dashed curve) schemes, are shown, at = 10, for the slow variable, on a uniform mesh with = 300.

Test case with transient demixing
Shock tube experiments can induce transient regimes where demixing occurs, see Figures 4 and 5 in [22]. Indeed, after strong sudden accelerations or shock front crossing, the orientation of mass fraction fluxes can transiently be opposite to the diffusive gradient closed flux of equation (2.14). This so-called counter-gradient situation causes a decrease of the TMZ width, which can be measured as = ∫︀ d̃︀ (1 −̃︀ ). In applications where the chronometry is important (e.g. ICF where only events occurring before maximum compression matter), this non-classical behaviour is important because it introduces a delay in the evolution of the TMZ.
In order to show that our model is able to handle transient demixing, we consider counter-gradient mass fraction flux as an initial condition, and analyse the evolution of resulting from the full hyperbolic-parabolic Figure 10. Self-similar test case: numerical and analytical solutions on a mesh having a geometric progression, at = 10, with = 50. The results from the first-order (solid curves), second-order (dashed curves) schemes and the analytical solution (dotted curves) are given for the mass fraction (blue), mass fraction fluxes (red), and mass fraction covariances (green), and the slow variable (purple). Table 6. Convergence rates corresponding to Figure 8, for the self-similar problem.
First-order scheme Second-order schemẽ︀ system. We choose the model constants according to Table 5, and propose two computations, at first order, differing only by the initial conditions: the first begins with a diffusive turbulent mass fraction flux (i.e. of the same sign as the first-gradient closed flux in Eq. (2.14)), and the second with an anti-diffusive flux, here taken as the opposite of the previous one. Let us consider one material, denoted by , of a two-material configuration.  We show, in Figure 11, the evolution of the integral mixing width , together with the evolution of the turbulent mass fraction flux at the center of the TMZ, from both initial conditions. As expected, a transient demixing is observed for the counter-gradient (anti-diffusive) configuration: the TMZ extent is decreasing for small times, see Figure 11, left-hand panel. Next, due to turbulent production related to the gradient of the mean mass fraction, second term in equation (2.10), the diffusive nature of the turbulent flux is restored with time and both computations should tend to a similar asymptotic state, as made clear by Figure 11, right-hand panel.
Since the turbulent flux drives the width of the TMZ, the latter reaches its minimal extent at the time when the sign of the mass fraction flux changes. This phenomenology cannot be accounted for when only diffusive closures are used.

Conclusion
The mathematical analysis of Reynolds Stress Models (RSM) dedicated to compressible turbulence is delicate: the wave-structure of wide-purpose RSMs, such as BHR [2] and GSG [19,20], is unknown, while they embed, due to the Reynolds and Favre averaging procedure, non-conservative products. We have proposed, in this article, a mathematical and numerical analysis of these models, and focused on a particular extension of them, though representative of a specific limit regime. This extension has been studied independently, but it can be used in the frame of a splitting of the comprehensive RSMs. We have justified, in terms of stability and accuracy, the use of an approximate, Godunov-type scheme, in 1D curvilinear coordinates. A broadband set of test-cases has been presented, that aims at covering all the identified relevant regimes. Hence, this work contributes to the mathematical and numerical analysis of RSMs, which deserves better understanding.
The multi-D extension of stable, high-order, and realizable 1D RSM discretization is a natural development, but it is not trivial. The proposed 1D scheme is a good candidate. It can indeed be interpreted as a scheme relying on nodal solvers. In an unstructured multi-D context, nodal solvers indeed enjoy a larger stencil than face-based solvers, encompassing neighbour meshes via nodes. In this respect, the resulting schemes can be considered closer to the genuine multi-D Riemann problem, especially for meshes having a large aspect ratio. We have paved the way towards a multi-D adaptation of the proposed scheme, into the GLACE or EUCCLHYD framework [13,25], as we intend to show in a future companion paper. In this adaptation, the material mass fraction and their fluxes would play respectively the role of the mean pressure and velocity of the Euler equation, in the acoustic nodal solver. Further, the mathematical entropies would play the role of the Gibbs entropy. We have also introduced an a posteriori limitation at high order, that may also bring an essential, entropy-compliant backstop, in multi-D flows with highly distored meshes.
Among other possible perspectives and applications, this work could be pursued by introducing a molecular inter-diffusion operator for the material mass fraction equations, see [1,42]. Such operator can be stiff, depending on the specific stages of an Inertial Confinement Fusion (ICF) capsule implosion, for instance, as it was shown in [40]. In the numerical analysis point of view, the interplay between the turbulent hyperbolic and the nonturbulent, molecular inter-diffusion part of the model will then have to be considered.  Let us show that if ≥ 0, then ≥ 0.