On Regularized Systems of Equations for Gas Mixture Dynamics with New Regularizing Velocities and Diffusion Fluxes

We deal with multidimensional regularized systems of equations for the one-velocity and one-temperature inert gas mixture dynamics consisting of the balance equations for the mass of components and the momentum and total energy of the mixture, with diffusion fluxes between the components as well as the viscosity and heat conductivity terms. The regularizations are kinetically motivated and aimed at constructing conditionally stable symmetric in space discretizations without limiters. We consider a new combined form of regularizing velocities containing the total pressure of the mixture. To confirm the physical correctness of the regularized systems, we derive the balance equation for the mixture entropy with the non-negative entropy production, under generalized assumptions on the diffusion fluxes. To confirm nice regularizing properties, we derive the systems of equations linearized at constant solutions and provide the existence, uniqueness and L2-dissipativity of weak solutions to an initial-boundary problem for them. For the original systems, we also discuss the related Petrovskii parabolicity property and its important corollaries. In addition, in the one-dimensional case, we also present the special three-point and symmetric finite-difference discretization in space of the regularized systems and prove that it inherits the entropy correctness property. We also give results of numerical experiments confirming that the discretization is able to simulate well various dynamic problems of contact between two different gases.


Introduction
Multicomponent compressible gas mixture dynamics is an important field in science and engineering, and a number of systems of partial differential equations (PDEs) were developed to describe phenomena of such type, see, in particular, references [1-3] and references therein.
Numerical methods serve as the most powerful tool to solve and simulate such systems of quasilinear PDEs. Originally, various numerical methods were designed to solve the compressible single-component gas dynamics systems of PDEs, and vast literature is devoted to this subject, see, in particular, references [4][5][6] and references therein.
Preliminary regularization of equations is an important and frequently used approach in constructing numerical methods for solving various scientific problems. In computational physics, those regularizations that have a physical basis are usually preferred. Such numerical methods are also used in computational gas dynamics. These include explicit in time conditionally stable and symmetric in space finite-difference and finite volume methods without limiters based on the discretization of regularized, or the so-called quasi-gas-dynamic (QGD), equations of gas dynamics. It is well known that, without regularization, methods of such type are unstable. These QGD equations were originally constructed on the basis of the Bhatnagar-Gross-Krook model kinetic equations, see monographs [7][8][9]. They can be rewritten in the form akin to compressible Navier-Stokes-Fourier equations with artificial coefficients of viscosity and heat conductivity and additional second order terms in space representing the regularizing velocity, viscous stress and heat flux, with a small parameter τ > 0. These equations can also be obtained on the basis of compressible Navier-Stokes-Fourier equations using formal procedures of time averaging and expansion [10][11][12]. Numerical methods based on the QGD equations have been successfully tested in practice for almost 40 years, including complex applied problems; see an extensive bibliography in the above monographs and numerous subsequent works, among which we highlight only a few works devoted to 3D turbulence and magnetohydrodynamics problems and the inclusion of such methods in the well known open source software package OpenFOAM [13][14][15][16]. Note that τ is taken proportional to the characteristic spatial mesh step in numerical methods. The QGD equations were proved to be physically correct, in the sense that they imply the correct entropy balance equation, i.e., with a non-negative entropy production. Some mathematical regularizing properties of the QGD equations were also confirmed, including their Petrovskii parabolicity, in contrast to the Euler equations and compressible Navier-Stokes-Fourier equations, which have hyperbolic and composite hyperbolic-parabolic types, respectively, as well as L 2 -dissipativity of the QGD equations linearized on constant and equilibrium solutions [17,18]. Recall that the important question of correct setting of boundary conditions is closely related to the type of a system of PDEs, and this setting is usually the most complicated in the hyperbolic case and the simplest in the parabolic case. In addition, conditional stability theorems were proved in the linearized statement for the above mentioned difference methods based on the QGD equations [19][20][21]. Notice also a quasi-hydrodynamic (QHD) regularization which can be considered as a simplification of the QGD one applicable to some subsonic or transonic flows [8,9].
There are other regularizations of the gas dynamics equations, which were also studied mathematically and aimed at constructing new numerical methods; in particular, see three approaches [22][23][24][25] and [26,27]. In all the approaches, much attention is paid to the entropy correctness of regularized equations. Among the listed approaches, the last one based on the so-called bi-velocity hydrodynamics [28,29], is closest to the QGD approach in structure of equations, although they are far from being the same. Alternative approaches also demonstrate success, but so far they have not yet undergone such extensive multi-year testing as the QGD approach. This paper is related to further development of the QGD and QHD regularizations and the corresponding numerical methods in the case of multicomponent gas mixture dynamics, and we deal with the so-called one-velocity and one-temperature homogeneous gas mixture model, with the perfect polytropic components. The reason is that such type models are widespread in practice including the computational design of aircraft and rocket engines.
For binary mixtures, the original regularized QGD multi-velocity and multi-temperature homogeneous gas mixture model was constructed on the basis of kinetic equations for mixtures in ( [8], Chapter 9). It was rewritten in [30] in the form of the compressible Navier-Stokes-Fourier equations with exchange terms for components and additional regularizing velocities, viscous stresses and heat fluxes, and a justification of its entropy correctness was given. Concerning applications, in particular, see [8,31]. For multicomponent mixtures, see similar QGD model and its applications in [32,33]. We do not touch this model here. The transition to the QGD one-velocity and one-temperature model was accomplished in ( [34], Section 1) by aggregating the PDEs of the original model. The aggregation procedure is simple and consists in using the balance PDEs for the mass of components and the momentum and total energy of the mixture, followed by taking the common velocities and temperatures of the components in them. The main advantage of this procedure is that the entropy correctness of the resulting QGD model is guaranteed. A mathematical analysis of such a multicomponent QGD regularization, as well as its QHD simplification, with additional allowance for diffusion fluxes between components, has recently been given in [35].
In [36], another approach to the construction of regularized QHD Navier-Stokes-Fourier-Cahn-Hilliard equations at low Mach numbers was given, based on the wellknown Coleman-Noll procedure and also ensuring the entropy correctness of the QHD model. The regularizing velocity w α of the component α, which play an important role in the QGD and QHD regularizations, turned out to be different in [34,36] depending on the partial pressure p α and the total pressure p, respectively. An additional full or partial averaging of w α from [34] can be applied (for the QHD regularization, they are the same), which also makes the result depending on p, see the next Section. A mathematical analysis of such a multicomponent QHD regularization has recently been given in [37]; among other things, it turned out that, in contrast to the single-component case, in the absence of diffusion fluxes, the QHD system of PDEs acquires the composite hyperbolicparabolic type, i.e., the regularization becomes incomplete. Theoretical constructions were accompanied by experiments with corresponding difference schemes, see [34,[38][39][40][41], etc. The full averaging of w α seems to be unsuccessful since the entropy correctness of the QGD regularization with it was established in ( [34], Section 2) only by passing to a non-conservative modification of the balance PDE for the total energy of the mixture, and the corresponding difference schemes gave satisfactory results only in simple 1D tests, in particular, see Section 5 below.
The partial averaging of w α corresponds to the combined w α , and the first sufficiently successful experiments with corresponding difference schemes are presented in [42]. An incomplete attempt to derive the combined w α and the entire regularized system using the approach from [10] is also made there. In this case, the regularizing viscous stress tensor and heat flux turn out to be different than in ([34], Section 1) in contrast to the single-component case, and a significant drawback of such a system in [42] is the loss of entropy correctness.
In this paper, we analyze the effect of these new combined regularizing velocities of components w α depending on the densities of the components and the total pressure. However, we still apply the same aggregated regularized balance PDEs for the momentum and total energy of the mixture as in ( [34], Section 1.2) and [35] in the case of binary and general multicomponent mixtures, respectively. In addition, we involve a new generalized form of the diffusion fluxes between the mixture components. Following [35], we study both the QGD and QHD regularizations in a unified manner by introducing a parameter in the corresponding PDEs. The first main theoretical result of the paper is the derivation of the balance equation for the mixture entropy with non-negative entropy production for our essentially modified system of equations. The second result concerns the derivation and study of the linearized system of PDEs: we justify the existence, uniqueness and L 2 -dissipativity of weak solutions to an initial-boundary value problem for this system. We also discuss that our results imply the Petrovskii parabolicity of the original quasilinear system of PDEs which allows one to obtain the local-in-time classical unique solvability of the Cauchy problem for this system identical to ([35], Theorem 3.3) and simplify the statement of correct boundary conditions for it. We emphasize that the presence of the diffusion fluxes is crucial for validity of the second and related results, for, without them, the original system of PDEs becomes a more complicated composite hyperbolic-parabolic, as in [37], not parabolic. This discovered regularizing role of the diffusion fluxes is nontrivial and even somewhat surprising. Notice that important mathematical results on the properties of other PDEs for compressible heat-conducting gas mixtures were proved, in particular, in [2, [43][44][45][46].
In the one-dimensional (1D) case, we also consider the new special three-point and symmetric finite-difference discretization which modifies one suggested in [41] in the case of the new regularizing velocities and more general form of the diffusion fluxes; for the single-component gas dynamics PDEs, this discretization was suggested, generalized and computationally tested in [47][48][49]. The discretization uses some non-standard nonlinear averages of the densities of the components and temperature and is conservative in the mass of the components and the momentum and total energy of the mixture. Our main theoretical result relating to this discretization is the semi-discrete balance equation for the mixture entropy with the non-negative entropy production; it means that the constructed discretization is entropy correct. In addition, results of our numerical experiments demonstrate better (sometimes, much better) or not worse behaviour, depending on the test, for the new combined regularizing velocities compared to those used previously. Now we can hope that the entropy correct discretizations of the considered type will be further developed for the general multidimensional gas mixture dynamics PDEs in the spirit of [48].
Vast literature is devoted to other numerical methods for solving multicomponent gas dynamics PDEs. We refer the reader to the brief review and a collection of references in the recent paper [50]. Note that only a few of the papers touch the entropy correct methods [51]. This is an important but complicated subject even in the case of singlecomponent gas dynamics PDEs, see, in particular, reviews: Tadmor The structure and the results of the paper in more detail are as follows. In Section 2, we present the aggregated regularized systems of PDEs describing the multidimensional one-velocity and one-temperature homogeneous gas mixture model, define the collection of all the involved functions and pass to the combined regularizing velocities. Proposition 1 concerns properties of the average gas mixture parameters, and Proposition 2 establishes a useful particular connection between solutions to the regularized systems of PDEs for the gas mixture dynamics and single-component gas dynamics. The main result is Theorem 1 about the balance equation for the mixture entropy with the non-negative entropy production. In Section 3, we derive and study the linearized system of PDEs. The key role belongs to the properties of symmetry/skew symmetry and positive definiteness of the related bilinear forms considered in Lemma 2. Theorem 2 states the existence, uniqueness and L 2 -dissipativity of weak solutions to an initial-boundary value problem for the linearized system. We also discuss the Petrovskii parabolicity of the original quasilinear system of PDEs and a local-in-time classical unique solvability of the Cauchy problem for this system. In Section 4, we pass to the 1D case of the regularized system of PDEs, introduce the mesh notation and present a special three-point and symmetric discretization in space for 1D regularized systems. Theorem 3 contains a semi-discrete balance equation for the entropy of the gas mixture, with a non-negative entropy production, and serves as a counterpart of Theorem 1. Section 5 is devoted to 1D numerical experiments. Applying the constructed discretization, we solve four known tests from [52][53][54][55]. The results confirm that the discretization is able to simulate well various dynamic problems of contact between two different gases, including the case of high initial pressure drops, and have some advantages over other choices of the regularizing velocities from [34] and especially from [37,39].
The paper also contains four appendices. Appendix A is devoted to derivation of the combined regularizing velocities and the full regularized system of PDEs from [39] based on the Euler-type system of PDEs for multicomponent gas mixture dynamics, by applying a formal procedure suggested in [11]. In Appendix B, we accomplish the scaling of the regularized system of PDEs from Section 2 that is often used to solve practical problems. In Appendix C, for the 1D case of the Euler-type system of PDEs from Appendix A, the Rankine-Hugoniot relations on the shock wave are given, and conditions for the existence of a stationary shock wave and the relationship between the values of the sought functions to the left and right of it are derived. Finally, in Appendix D, the 1D finitedifference counterpart of Proposition 2 is given.

A Regularized System of Equations for the Multicomponent Gas Mixture Dynamics with New Regularizing Velocities in the Presence of Diffusion Fluxes
The aggregated regularized system of PDEs for one-velocity and one-temperature multicomponent homogeneous gas mixture dynamics consists of the following balance PDEs for the mass of components, total momentum and total energy of the mixture: Here, the main sought functions are the densities of the mixture components ρ 1 > 0, . . . , ρ K > 0 (K 2 is their amount), their common velocity u = (u 1 , . . . , u n ) and absolute temperature θ > 0. These functions depend on x = (x 1 , . . . , x n ) ∈ Ω and t 0, where Ω is a domain in R n , n = 1, 2, 3, and α = 1, K means that α = 1, . . . , K. Vector-functions are written in bold. The operators div = ∇· and ∇ = (∂ 1 , . . . , ∂ n ) are taken in x, ∂ t = ∂ ∂t and ∂ i = ∂ ∂x i . In this section, the symbols ⊗ and · denote the tensor and scalar products of vectors, the tensor divergence is taken with respect to its first index, and · is the operation of summation over index α = 1, K.
This regularized system of PDEs was derived in [34] for K = 2 and d 1 = d 2 = 0 by aggregating the regularized multi-velocity and multi-temperature gas mixture PDEs [30]. In general case, it has recently been studied mathematically in [35] in general case. Now, we sequentially define a number of functions involved in these PDEs. We assume that the mixture components are perfect polytropic gases and exploit the following expressions for the pressure, specific internal energy, the total energy and the specific enthalpy of the component α: with physical constants γ α = R α c Vα + 1 > 1, R α > 0, c Vα > 0 and c pα = c Vα + R α = γ α c Vα , the last two of which are the specific heat capacities at constant volume and pressure, α = 1, K. One can consider any two of four constants γ α > 1, R α > 0, c Vα > 0 and c pα > 0 as the main independent ones; below, in computations in Section 5, such role is played by γ α and c Vα .
The total density and pressure, average specific internal energy and total energy of the mixture are expressed by the formulas with the average gas mixture parameters The second Formula (5) is the Dalton law for mixtures. The function ρ α ρ =: C α is the mass concentration of the mixture component α. Consequently, the important formula of the standard form for the total pressure holds as well In contrast to the single-component case, R, c V and γ are functions, not constants, except for the particular cases R 1 = . . . = R K , c V1 = . . . = c VK and γ 1 = . . . γ K , respectively.
In the above PDEs, the following regularizing velocities for the component α and average ones were originally used see [34], where τ = τ(ρ, u, θ) > 0 is a regularization (relaxation) parameter which is usually a function, not constant, with ρ := (ρ 1 , . . . , ρ K ). Here, = 0 or 1, and the regularization is of the so-called quasi-gasdynamic (QGD) type for = 1 or essentially simpler quasihydrodynamic (QHD) type for = 0, so actually we consider two different, albeit related, systems in a unified manner similarly to [35]. The Formula (8) mean takes the average of w α and w α , in other words, the full and partial averaging of w α . In this paper, we replace w α by w making w α combined and dependent on the total pressure p instead of the partial pressure p α : and analyze the effect of this replacement discussed above in Introduction. Notice that the replacement does not affect the validity of Formula (8) for w . The viscosity tensor and heat flux are expressed, respectively, by the formulas and contain the standard-type terms and the regularizing ones with the superscript τ. The classical Navier-Stokes viscosity tensor and the Fourier heat flux are given by the formulas where µ > 0, λ 0 and κ > 0 are the total coefficients of dynamic and bulk viscosities and heat conductivity (which can depend on the sought functions (ρ, u, θ)), ∇u = {∂ i u j } n i,j=1 and I is the unit tensor of order n. Next, the regularizing viscosity tensor and heat flux are given by the formulas The density of body force f and intensities of heat sources Q α 0 (acting on the component α) are given functions, and Q := Q α 0. Finally, we consider the diffusion fluxes and additional respective heat flux of the form where · β means the summation over index β = 1, K. The functions G α and s α are the usual Gibbs potential and specific entropy of the component α, and s α0 , ρ α0 > 0 and θ 0 > 0 are constant reference values for s α , ρ α and θ, α = 1, K. The functions-coefficients d αβ and e α can depend on the sought functions. Their specific form is not essential below, and we only assume the symmetry property d αβ = d βα for any α = β.
Let · α,β := · β mean the summation over α, β = 1, K. Using permutations of indices α and β and then this symmetry property, we obtain two identities for any numbers ϕ 1 , . . . , ϕ K and ψ 1 , . . . , ψ K . The first identity implies the important physical property d α = 0, and the second one will also be essential below.
We can avoid the explicit usage of s α due to the formulas In particular, for e α = s α , we obtain the simplest formulas The general multicomponent case (K 2) for = 0, 1 in the presence of d α and q d has recently been studied mathematically in [35] but only in the particular case d αβ ≡ d 0 and Ke α − e β β = b α , that is the same for K = 2 but much less general for K 3. The above quantities d α and q d generalize those proposed in [1] in the case K = 2; in this case, the formulas are transformed and discussed in more detail in [35]. Notice also that a much more general approach for introducing these quantities is known, for example, see [2].
In [34], the above total coefficients µ, λ and κ are defined simply as i.e., the sums of the corresponding coefficients of the components. These coefficients can be artificial depending on τ in order to ensure stability of symmetric in space discretizations for computations, or physical, or sums of them. In the first case, the typical formulas for τ and them are as follows in accordance with Formula (18). Here, 0 < a 1 is a parameter, a Sα > 0 and a Prα > 0 are the Schmidt and inverse Prandtl numbers for the component α; a 1Sα 0 is a counterpart of a Sα (in particular, a 1Sα = 0), which can be also used as adjusting numerical parameters, is the sound speed of the mixture, i τ = 0 or 1, and h is a characteristic size of the spatial mesh. In the case of a Sα = a S , a 1Sα = a 1S and a Prα = a Pr independent of α, the formulas for µ, λ and κ are simplified: For the single-component gas dynamics, see such formulas, in particular, in [8,9,19].
Recall that R α = R 0 m α , where R 0 is the universal gas constant and m α > 0 is the molecular mass of gas α. In some cases, γ α and m α are taken as the two main gas constants, and the average molecular mass of the mixture m is defined by 1 m = ρ α ρ 1 m α . Then the other gas constants can be expressed in the form We first give some inequalities for R, c V , m and γ.
Proposition 1. 1. The two-sided bounds hold Moreover, all the inequalities are strict, except for the particular cases R 1 = . . . = R K , c V1 = . . . = c VK and m 1 = . . . = m K , respectively. 2. The formula and two-sided bounds also hold Moreover, both the bounds are strict excluding the particular case γ 1 = . . . = γ K ; in that case, we have γ = γ 1 even if c Vα are not identical for all α = 1, K.
3. The following relations hold Proof. Items 1 and 2 are elementary. Item 3 is valid since owing to the Cauchy inequality The inequality becomes an equality only for (22), we can accomplish the following transformations

Remark 1. Starting from Formula
Permuting indexes α and β, similarly to identity (16), we derive the representation This formula also implies Item 3. (20) and (21) show that both of them are averages of γ 1 , . . . , γ K , and clearly c Vα and R α can be scaled in these formulas. Forγ, the same bounds as in Item 2 for γ are valid. Note that, in our expression for Π τ , see (11), we use the term γ α p α as in [34,35], in contrast to γp in [39]. Moreover, both sides of the second inequality (21) can be considered as different definitions for the squared sound speed of the mixture, and we prefer to use the right-hand side in this role in this paper (Appendix C helps to make the choice).

Formulas for γ andγ in relations
We first consider the sought functions of the particular form.
Proposition 2. Let 0 < C α < 1, α = 1, K, be arbitrary constants such that C α = 1. Consider the sought functions of the particular form ρ α = C α ρ with ρ > 0 (α = 1, K), u and θ > 0 and the case of d α = 0 and ψ = ψ(ρ, u, θ) for the functions ψ = τ, µ, λ, κ. For them, the above regularized system of PDEs for the gas mixture dynamics is reduced to the following regularized system of PDEs for a single-component gas dynamics for the sought functions ρ, u and θ and ( whereγ = C α R α γ α / C α R α and the above formulas (10) for Π NS and q F are in use.
Proof. For ρ α = C α ρ, under the assumptions made about C α , we clearly obtain ρ = ρ α and Thus, all the balance PDEs for the mass of components (1), after division by C α , are reduced to Equation (23). Moreover, expressions (11) for Π τ and (12) for −q τ are reduced to those given in Formulas (27) and (28), using the expression forγ in relations (21). Therefore, now the original balance PDEs for the total momentum and total energy (2) and (3) take forms (24) and (25).
This proposition establishes a particular connection between solutions to the regularized systems of PDEs for the gas mixture dynamics and single-component gas dynamics for any γ α > 1 and c Vα > 0, α = 1, K, and can be useful to check properties of the former system. In fact, it enlarges the corresponding 1D Proposition 1 in [41]. However, recall that γ = γ in Formula (27) only in the particular case γ 1 = . . . = γ K , see Proposition 1, Item 3 or Remark 1.
Applying the operation · to the mass balance equation for the mixture components (1), using the formula ae α (u − w α ) = ρ(u − w ) valid according to the first expression (8), and the property d α = 0, we obtain the important total mass balance equation Here, ρ 1 , . . . , ρ K and θ appear only implicitly since ρ = ρ α and p in w depend on them. The balance PDEs for the total momentum and total energy of the mixture (2) and (3) entail sequentially the balance PDEs for the kinetic and internal energies of the mixture where the symbol : denotes the scalar product of tensors. The derivation exploits the total mass balance equation (29) and is valid for any w , for the former equation, and exploits only the relation w = τ div(ρu)u + w (for f ≡ 0), but not explicitly Formulas (7) and (8), for the latter equation, see ( [35], Section 4.1).
The first main result of the paper concerns the total entropy balance equation. Recall that the specific entropy of the mixture is given by the formula s = ρ α ρ s α . The result corresponds to ( [35], Theorem 3.1) but concerns another definition of the regularizing velocity (9) and deals with much more general form of d α , for K 3; this form is applicable in [35,37] as well.
Theorem 1. Let d αβ = d βα 0 for any α = β. The following regularized entropy balance equation for the multicomponent mixture in the presence of diffusion fluxes holds with the entropy production P NS + P τ , where and | · | F is the Frobenius norm. Moreover, P τ is non-negative for = 0, as well as for = 1 under the condition This condition is certainly true provided that τ( Proof. According to ([35], proof of Theorem 3.1), the following preliminary equation involving the entropies of the mixture and the components holds where This equation is derived from the balance PDEs (1) and (31) and does not exploit specific expressions for w α and w in them. In the first term on the right, we have taken into account the following obvious formula see definition (14) of q d , that only slightly differs from the similar formula in [35]. Using identity (16), we can write the first and second terms on the right in Equation (32) Next, in the case of expressions (9), we can extract and collect the terms with w from the first, third and fourth terms of B τ and thus write Concerning the remainderB τ , the following formula holds 1 θB τ = P τ − ρ τθ | w| 2 , see ( [35], proof of Theorem 3.1) and the references therein, that completes the proof.
Clearly, P NS + Q θ and P τ − Q θ are the Navier-Stokes-Fourier and regularizing contributions to the entropy production. Theorem 1 remains valid for τ 0 (in particular, τ = 0, i.e., without a regularization), when one should pass to a different form for the first relaxation term:

An Auxiliary Reduction of the Balance Equations
Let f = 0 and Q 1 = . . . = Q K = 0. We introduce the vector of the sought functions z := (ρ, u, θ) and first present an important auxiliary reduction of the balance PDEs for ρ α , u and θ up to the terms O(|∇z| 2 ). Lemma 1. The following reduced PDEs hold: for the densities of the components for the velocity with χ := 1 3 µ + λ, and for the temperature Hereafter, ∆ = div ∇ is the Laplace operator andē α = R α − s α + e α , see the last formula (17).
Proof. According to ([35], Section 3.2), the balance equation for the velocity holds and the balance equation for the temperature holds as well. Their derivation does not exploit a specific form of w α . For d 1 = . . . = d K = 0, the presented reductions have recently been proved in [37] for = 0, and the similar reductions have also been accomplished in [35], where the terms with multiplier are the same as here, whereas the other terms are partially different. So, it suffices to reduce the terms with d α . In the balance PDEs for the mass of components (1), we obtain since the following chain of transformations is valid and c pα − c Vα = R α . In the balance equation for the temperature (37), we can write using the previous decomposition (38). This completes the proof.
The systems of PDEs (1)- (3) and (1), (36), (37) (taking into account formula (11)) are equivalent for classical (smooth) solutions. Below the reduced system of PDEs (33)- (35) helps to linearize the original system of PDEs and perform its parabolicity analysis. Clearly, the left-hand sides of these PDEs are independent of τ and .

Linearized Regularized System of PDEs, Its Properties and Corollaries
In the case f = 0 and Q 1 = . . . = Q K = 0, the system of PDEs (1)-(3) has constant solutions and any u 0 . We perform the linearization of the solution z to this system on z 0 and write where z := ( ρ, u, θ) with ρ := ( ρ 1 , . . . , ρ K ) is the vector of dimensionless perturbations and ρ α * > 0, u * > 0 and θ * > 0 are scaling parameters selected below. We substitute the solution in this form into the reduced system of PDEs (33)- (35) and discard the terms of the second order of smallness with respect to z and its first and second order derivatives using the formula ∇z = (ρ 1 * ∇ ρ 1 , . . . , ρ K * ∇ ρ K , u * ∇ u, θ * ∇ θ). Then, we divide the resulting PDEs by ρ α * , u * and θ * , respectively, and derive the linearized system of PDEs with constant coefficients for z: Here, moreover, the scaling factors u * and u 2 * are taken out of the convective and dissipative terms (i.e., the terms with the first and second order derivatives except for the diffusion terms), respectively, and the following notation is introduced for the components of the scaled background solution, background values of ρ, R and c V and the average value of R α γ α : In addition, d αβ0 ,ē α0 , τ 0 , µ 0 , χ 0 and κ 0 are the values of d αβ ,ē α , τ, µ, χ and κ, respectively, on the background solution z 0 , and the following PDE operator is involved For d 1 = . . . = d K = 0, the possibility of simultaneous symmetrization of the convective and dissipative terms has recently been found in [35,37] by choosing the scaling parameters with a free parameter b. We accept this choice and pass to a much simpler form of the above linearized system of PDEs Here, the following constant factors have been introduced and, for the last two terms in (40) and (42), we have taken into account the formulas Next, we study the initial-boundary value problem (IBVP) for the linearized system of PDEs (40)- (42) in the cylinder Q ∞ := Ω × (0, ∞) under the boundary and initial conditions Let L 2 (Ω) and L 2 (Ω) be, respectively, the standard Lebesgue spaces of functions and vector functions defined on Ω and denote by (·, ·) Ω = (·, ·) L 2 (Ω) , · Ω = · L 2 (Ω) , (·, ·) Ω = (·, ·) L 2 (Ω) and · Ω = · L 2 (Ω) their inner products and norms. Let H 1 (Ω) = W 1 2 (Ω) be a standard Sobolev space of vector functions defined on Ω, and H 1 0 (Ω) be the closure in the H 1 (Ω)-norm of the space of smooth vector-functions with a compact support in Ω.
In the identity, the three bilinear forms are involved where the tensors ∇ u and ∇u are considered as vectors of length n 2 , and The last bilinear form corresponds to the diffusive terms and is independent of u and u. Formally, identity (44) arises after multiplying Equations (40)-(42) by ρ α , u and θ, respectively, integrating over Ω and by parts and summing up the results. Let us study properties of the defined bilinear forms that is crucial below.
We define the weak solution z ∈ V(Q T ), for any T > 0, to the IBVP for the system of PDEs (40)- (42) in Q ∞ together with conditions (43), such that the integral identity T 0 ∂ t z(·, t), z(·, t) Ω dt + u * B Q T ( z, z) + u 2 for any z ∈ L 2 ((0, T); H 1 0 (Ω)) and any T > 0, together with the initial condition z| t=0 = z (0) ∈ L 2 (Ω) are valid. Here, the inner products in the bilinear forms B Q T , A Q T and A d Q T are taken over Q T instead of Ω as originally. Due to the well known embedding V(Q T ) ⊂ C([0, T]; L 2 (Ω)) [56], the initial condition is understood by continuity in L 2 (Ω). Formally, identity (51) arises from the previous one (44) for z = z(·, t) by integration over (0, T). Now, we are ready to state the second main result.
The proof is based on some general results from [56] together with Lemma 2 and Corollary 1 and is quite similar to that of ( [35], Theorem 3.2 and Corollary 3.1), so we omit it. In the case d = 0, the Cauchy problem can be considered similarly to ([37], Theorem 2).
Due to the equivalence (for classical solutions) between system of PDEs (1), (36) and (37) and the original quasilinear system of PDEs (1)-(3) for multicomponent gas dynamics introduced in Section 2 and some very general results from [58], only this Petrovskii parabolicity property implies the local-in-time classical (in anisotropic Hölder spaces) unique solvability of the Cauchy problem to the latter system. Its statement is identical to that of the corresponding ( [35], Theorem 3.3). Moreover, the proof of the Petrovskii parabolicity property can be made very close as well, using the remark from ( [35], proof of Lemma 3.2) that the uniform in D positive definiteness of a matrix A(z 0 , ξ), defining the property is equivalent to a result such as Corollary 1 in the case Ω = R n , which goes back to a technique using the integral Fourier transform from ( [17], Section 3) (where the Petrovskii parabolicity was studied in the single-component gas case). Notice that it was shown in [35] that the matrix A(z 0 , ξ) can be symmetrized by the same scaling as accomplished above for passing to the linearized system of PDEs (40)- (42), and thus this symmetry property is directly connected to the symmetry of the bilinear forms in Lemma 2. For brevity, here we omit details of both the statement and proof of such theorem.
In addition, the Petrovskii parabolicity property allows one to pose correctly some simple boundary conditions in IBVPs for the original system of PDEs, similarly to the single-component case. We emphasize that the presence of the diffusion fluxes is crucial for validity of this property, for without them, the original system of PDEs becomes more complicated composite hyperbolic-parabolic, not parabolic. For = 0, this has recently been checked in detail in [37]. The reason is that the quadratic form A Ω (z, z) is only nonnegative rather than positive definite in H 1 0 (Ω), and the degeneration occurs with respect to ∇ρ. For = 1, this can be performed similarly and, moreover, this is clear if D contains a point z 0 = (ρ 0 , 0, θ 0 ), since then the quadratic form A Ω (z, z) is identical for = 0 and 1 at such point z 0 , withχ 0 + τ 0 a 2 substituted forχ 0 . For systems of composite hyperbolicparabolic type, results of type ( [35], Theorem 3.3) are not valid any more, and the statement of correct boundary conditions in IBVPs for the original system of PDEs becomes more complicated and has not yet been studied.

The 1D Regularized System of PDEs for Gas Mixture Dynamics and Its Entropy Correct Spatial Discretization
Starting from this section, we pass to the particular 1D case of the above regularized system of PDEs for the gas mixture dynamics, and its constituent balance PDEs for the mass of the components and the momentum and total energy of the mixture take a simpler form with = 0, 1. The main sought functions are the component densities ρ 1 > 0, . . . , ρ K > 0 and their common velocity and temperature u and θ > 0 depending on (x, t) ∈ [−X, X] × [0, T]. We exploit the previous Formulas (4)-(6) for the main gas state variables and the total energy of the components and the mixture, where now |u| 2 = u 2 . The 1D formulas for the regularizing velocities, viscous stress and heat flux look as follows where ν := 4 3 µ + λ and α = 1, K. As above, introducing the average regularizing velocity w := ρ α ρ w α = τ ρ u∂ x (ρu) + w allows one to simplify the form of the balance PDEs for the momentum and total energy of the mixture (54) and (55): However, for further discretization, we prefer to use the original form of these PDEs since this approach allows us to derive a counterpart of Theorem 1. Let us first introduce the mesh notation. Define the uniform meshω h on [−X, X], with the nodes x i = −X + ih, 0 i N, and the step h = 2X N . Let ω h =ω h \{−X, X} be its internal part. Define also an auxuliary mesh ω * h with the nodes x i+1/2 = (i + 1/2)h, 0 i N − 1.
Let H(ω) be the space of functions defined on a mesh ω. We introduce the shifts of the argument v −,i+1/2 = v i and v +,i+1/2 = v i+1 and the averages and difference quotients and y i+1/2 = y(x i+1/2 ). First, for simplicity, let there be no body force (i.e., f = 0). Following [41,47], we apply a non-standard spatial discretization of balance PDEs for the mass of the components, the momentum and total energy of the gas mixture (53)-(55) and construct their following three-point and symmetric semi-discrete counterparts The main sought functions ρ 1 > 0, . . . , ρ K > 0, u and θ > 0 together with the functions p α , ε α and E α are defined in space on the main meshω h . In the equations, the above expressions (4) for p α , ε α and E α as well as (5) for ρ, p, ε and E, with |u| 2 = u 2 and the coefficients R and c V from Formula (6), are exploited.
We apply the following discretizations of the regularizing velocities (56): as well as of the viscous stress, and heat flux and diffusion fluxes (57)- (60): The functions w α , w, Π α , q α , τ, ν α , κ α and Q α are defined in space on the auxiliary mesh ω * h . Moreover, G α , the Gibbs potential of the component α, see Formula (15), is defined in space onω h , whereas the functions d α , e α and q d are defined in space on ω * h . Here, we apply nonstandard averages of ρ α , p α , E α and ε α of the form [41,47] [ that exploit the divided difference for the logarithmic function (67). Note that u − u + is similar to the geometric mean for u 2 (although it is negative for u − and u + of different signs). Concerning the case of τ = T ρ, u, θ , one can set, in particular, τ = T [ρ], [u], [θ] or τ = T (ρ, u, θ) in space on ω * h . In computations in Section 5 below, we apply the second formula.
This spatial discretization is close to a similar one recently constructed in ( [41], Section 5) and differs from it by expression (66) for w α (approximating formulas (9), not (7), in the 1D case) and the much more general expression (70) for d α in the case K 3. In its turn, this discretization in [41] generalises the original one from [47] in the case of the single-component gas dynamics to the considered multicomponent gas mixture dynamics PDEs.
Notice that the arising semi-discrete counterparts of Formula (17) are nontrivial: since The main result in this section is a 1D semi-discrete counterpart of the balance equation for the mixture entropy, see Theorem 1. It corresponds to ([41], Theorem 2) but concerns another definition (66) of the semi-discrete regularizing velocity and deals with a much more general form of d α , for K 3; this form is applicable in [41] as well.
The term P NS h + P τ h * in Equation (71) is the semi-discrete entropy production. The first three terms of P τ h are non-negative, and the last term is non-negative for = 0, as well as for Q. This condition is certainly true provided that Proof. The following semi-discrete balance equations for the total mass and kinetic and internal energies of the mixture hold ∂ t ρ + δ * j = 0, 0.5∂ t (ρu 2 ) + 0.5δ * (j u − u + ) + (δ * [p])u = (δ * Π )u, on ω h × [0, T], with j := j α . They are counterparts of the balance PDEs (29)-(31) and have recently been proved in ( [41], Lemma 3), and their derivations remain valid for any w α , Π and q , in particular, given by expressions (66)-(68). According to the proof of ( [41], Theorem 2), the following preliminary 1D semi-discrete balance equation for the mixture entropy holds where we have taken into account that e α plays the role of K −1 b α in the definition of q d in [41]. Recall that its derivation starts from the semi-discrete balance equations for the mass of components (63) and the internal energy of the mixture (72), and the specific form of w α does not matter in this derivation. Here, the following term and its decomposition Applying them and then identity (16), we obtain Next, using expressions (66), we can extract from A α the term [θ]A α such that and, according to ( [41], Appendix A) or [47], the following formula holds Applying the operation · to it and accomplishing the transformations we complete the proof.
As in the differential case, the entropy production remains non-negative for τ 0, where one should pass to another form for the first relaxation term in P τ h inside the curly brackets: [ρ] ([ρ][u]δu + δp) 2 . At the end of the section, following [41,47], we generalize the constructed semidiscrete method and Theorem 3 to the case of any f . Recall the general momentum and total energy balance PDEs (54) and (55) and expressions for the regularized velocities (56), and generalize the semi-discrete Equations (64) and (65) by adding, respectively, the terms to their right-hand sides, with the functions ρ * := [ρ] − τδ(ρu) and f defined in space on ω * h . We also generalize the expression for w as w = τ The new terms with f produce the following additional term on the right-hand side of the semi-discrete balance equation for the internal energy (72): To derive the semi-discrete balance equation for the mixture entropy (71), one should multiply this term by 1 θ and transform the result. The required transformation was accomplished in [41] for general w α = τ [ρ α ] [u]δ(ρ α u) + w α with any w α , and in our case where w α = w is independent of α, it leads to the formulas As a result, in the preliminary entropy balance equation (73) With the given f -dependent extensions, Theorem 3 remains valid.

Numerical Experiments
We are still dealing with the 1D system of PDEs. Let us compare three cases A, B and C of the regularizing velocities w α in the balance PDEs (53)- (55): with = 1 and α = 1, 2 (the case of binary mixtures), considered, respectively, in this paper (see Formula (56)), papers [34,35,41] (see also Formula (7) for n = 1) and [34,37,39]. Case A is the main one below, and we demonstrate some its advantages over the other two cases.
In case C, w α is independent of and α. We discretize these expressions, respectively, according to Formula (66) as well as [u]δ(ρu) + w.
Moreover, ρ 1r = ρ 2l = 0, although, in computations, we take them suitably small positive (equal 10 −10 ) instead. The parameters of the gases to the left and right of x = 0 and the final time of computations t f in can be found in Table 1. Note that there the simplest values c V1 = c V2 = 1 in Examples 1, 3 and 4 were not required originally (since in the case τ = ν = κ = 0, there exists a closed Euler-type system of PDEs for the sought functions ρ, γ, u and θ, see details in Appendix A), and we have chosen them ourselves. The initial temperature θ 0 is defined in accordance with Formulas (4) and (5): The boundary values of the sought functions at x = ±X in time are kept the same as their values given at t = 0. We also take X = 0.5 in Examples 1-3 and X = 5 in Example 4. In Examples 1-3, the initial pressure drop rapidly increases: We apply them for ε and ρ α , respectively, that leads to the formulas We introduce a non-uniform mesh in time 0 = t 0 < t 1 < . . . < t m = t f in , with the steps h tm = t m − t m−1 . We take the relaxation parameter and the artificial viscosity and heat conductivity coefficients according to Formula (19) and ν = 4 3 µ + λ: , ν = τ(a S1 p 1 + a S2 p 2 ), κ = τa Pr (γ 1 c V1 p 1 + γ 2 c V2 p 2 ).
Here, 0 < a < 1 is a parameter, h = 2X N , i τ = 0 in Examples 1 and 3 or i τ = 1 in Examples 2 and 4, and, for example, u m i = u(x i , t m ). In addition, a S1 = a S2 = 3 4 and a Pr = 1 in Examples 1-3.
To complement the above spatial discretization, we apply the simplest explicit Euler method for the temporal discretization, together with the automatic choice of the time steps h tm = t m − t m−1 according to where β is the Courant-type parameter. A linearized stability (more precisely, L 2 -dissipativity) conditions for such an explicit scheme theoretically and practically were studied in [19,49] in the single-component gas case. In every example, we adjust the parameters a and β. If the values of ρ 1 or ρ 2 less than 10 −10 arise at the upper time level, we replace them by 10 −10 . Notice that, with a code for the considered discretization of the single-component gas dynamics PDEs, as in [49], at one's disposal, it is not difficult to extend it to the case of the binary gas mixture. Proposition A2 in Appendix D (the 1D discrete counterpart of Proposition 2) was applied for initial testing of our code for mixtures.
Example 1. (the test from ( [54], p. 266)). In this first rather simple test, there is a contact discontinuity between the two gases that moves to the right with constant velocity u; the pressure p and temperature θ are also constant. The Mach number of the mixture ranges approximately from 0.14 to 0.42 and is not high. However, it is known that not all numerical methods are able to reproduce this constancy well, especially for moderate N.
In the main case A, . the results for a = 0.5, β = 0.7, N = 251 (similarly to [54]) and 1001 are given in Figure 1; note that the scales for p, u and θ are enlarged there. In this and other examples, we exhibit graphs of six functions: ρ 1 , ρ 2 , ρ, p, u and θ. The deviations from constant values near the contact discontinuity are very small in ρ 1 and p, slightly larger for u and θ. They diminish as N grows and, for N = 1001, they disappear for p, almost disappear in ρ 1 and u and become very small in θ. Hereafter, the graphs for two values of N in general almost coincide except for vicinities of the contact discontinuity and the shock wave; the differences become more visible after several magnifications. In case B, we need to take smaller a and especially β to obtain suitable results: a = 0.2 and β = 0.1. However, for N = 251, the results are still not so nice: the behaviour of ρ 1 , ρ 2 and ρ is too smooth near the contact discontinuity, the deviations from constant values near the contact discontinuity are very small for p, small for θ, but rather large in u. For N = 1001, the results become better, see Figure 2. In case C, the results for a = 0.2, β = 0.1 and N = 251, are worse than in case B: all the deviations mentioned above in it are significantly larger. Hereafter, we mainly omit the graphs in this case for brevity (except for Example 3). Table I). In the original paper, t f in was missed, and we adjusted its value to the graphs given there. The final solution contains the contact discontinuity between the two gases, with jumps in the values of ρ 1 , ρ 2 , ρ and θ, but not p and u, as well as a rarefaction wave in gas 1 to the left and a shock wave (the strong discontinuity) in gas 2 to the right of the contact discontinuity. The functions ρ 1 and p are non-increasing, whereas ρ 2 , ρ, u and θ are non-monotone, with the maximal values of ρ 2 , u and θ in front of the shock. In addition, ρ 2 is piecewise constant. The final maximal Mach number is M max := max i |u m i | c m si ≈ 1.67, so now the flow is partially supersonic (note that u = M = 0 closely to the boundaries). Note that the parameters are not scaled in this example in contrast to the rest of them, and, at first, we do not use scaling in our computations to check further our method's capabilities. We choose a = 0.5 and β = 0.4. Notice that, in this example, we can take zero artificial viscosity ν = 0 without essential changing the results.

Example 2. (a version of the Sod problem) from ([52],
In the main case A, the results for N = 501 and 2001 are given in Figure 3, and they correspond well to those from [52] excluding the very small hollow in ρ at the contact discontinuity. Notice that scaling u and θ by the natural divisors u * = θ * = 1000 (with no scaling of ρ 1 , ρ 2 and x) and consequently p by the divisor p * = 10 6 , see details in Appendix B, does not improve the results. In case B, once again, the results for the same N are rather nice, but the graph of ρ 2 is slightly more smoothed and the graph of ρ has an additional false step, both near the contact discontinuity. See Figure 4. In case C, the results for the same N are of low quality in general. Despite the fact that ρ 1 , p and θ are computed rather accurately, the graphs of ρ 2 and u have high "fingers" at the point of contact discontinuity. Additionally, there are single oscillations of relatively small amplitude in p and of high-amplitude in ρ near that point. Since the unknown functions satisfy the unified system of PDEs, it is somewhat surprising that some of them are computed accurately, while the rest are not. In this case, even in a simplified example, low quality of numerical results has recently been detected in [59] for another scheme.
We also briefly comment on two simpler Sod problems, see Examples 1 and 2 in [41], where the initial pressure drop is much less and equals 10 and 20 and the case B was studied only. If the initial pressure drop equals 10, results in cases A and B are very close and both equally correct. Even in case C, results are rather well, though with small ledges in the graphs of ρ 2 and u at the contact discontinuity.
If the initial pressure drop equals 20, the best results are in case A. In case B, they are also nice, though the graph of ρ 2 is slightly more smoothed and the graph of ρ has a slightly deeper hollow near the contact discontinuity. In case C, the results are already poor in general: though ρ 1 , p and θ are computed rather accurately, the graphs of ρ 2 and u have visible "fingers" at the point of contact discontinuity, and the graph of ρ has a single oscillation near that point. Example 3. (stiff two-gas shock-tube problem) from ( [53], Test 5.4). In [53], the values p l and p r were confused, and t f in was not specified so we adjusted its value. In this example, the initial pressure drop equals 2500 and is very high. In general, the behavior of the final solution is similar to the previous example. However, the support of the maximal value of ρ 2 is narrower, θ is non-increasing and the jumps in the values of ρ 2 , ρ and θ are high. Furthermore, M max ≈ 1.44, so the flow is partially supersonic once again. We take a = 0.25 and β = 0.4.
In the main case A, for N = 4001, the results correspond well to those in [53], see Figure 5. For smaller N = 1001, the graphs of ρ 2 and ρ are more smoothed near the contact discontinuity, as well as u and θ have very small ledges near the right end of the rarefaction wave and the contact discontinuity, respectively, though the rest of the graphs look well. In case B, the computation for the same β = 0.4 and N fails. For four times smaller β = 0.1, the results are not bad (see them for N = 4001 in ( [41], Example 3)), but now the graph of ρ 1 acquires an additional rather high, though very narrow, false step, and the graphs of ρ 2 and ρ are slightly more smoothed, both near the point of contact discontinuity.
In case C, the results are specific, see Figure 6. Now, for β = 0.2 and the same N, the graph of ρ 1 has a large "finger" at the contact discontinuity which height is about twice the value of ρ 1 to the left of that point. Surprisingly, the rest of graphs look rather accurate including even ρ (although after a magnification, the defect in its graph just to the right of the contact discontinuity becomes more noticeable), and the situation does not change for some larger times as well. This figure shows that the rather accurate computation of ρ, p, u and θ for the mixture does not guarantee the same concerning both ρ 1 and ρ 2 for the components (the latter graphs are sometimes omitted from the numerical results).    1.5698, 1.4), otherwise , thus, in gas 1 with γ 1 = 1.4, there is "the bubble" of gas 2 with γ 2 = 1.648 moving to the right.
Here, t f in = 4, and M min := min i |u m i | c m si ≈ 0.32 and M max ≈ 1.22, so the flow is transonic. In this problem, shock waves and contact discontinuities interact, that complicates computations. Let a = 0.5 and β = 0.2.
In the main case A, the results for N = 1001 and 4001 are given in Figure 7. For N = 4001, the graphs of ρ and p are very close to those only presented in [55]. The very small ledges in the graphs of ρ 2 are observed at the two points of contact discontinuities. To achieve their smallness, we have changed the values a S2 = 0.15 and a Pr = 10. In case B, the results for the same N are shown in Figure 8. On the one hand, the overall quality of the solution for N = 1001 is high, without any ledges, even for the same standard a Sα and a Pr as in previous examples. On the other hand, ρ 1 , ρ 2 , ρ and θ are too smooth near the two mentioned points, especially for N = 1001.
In case C, the graphs of ρ 2 have high "fingers" at the same two points which we could not remove by changing the parameters. Nevertheless, as in Example 3, the other graphs are correct except for a small "finger" in θ at the right of the same points.
Finally, we see that the results in the main case A are better or not worse than in the other two cases. The weakest results are in case C.

Data Availability Statement:
The datasets generated during the current study are available from the corresponding author on reasonable request.

Conflicts of Interest:
The authors declare no conflicts of interest. The funders had no role in the design of the study, in the writing of the manuscript or in the decision to publish the results.
By definition (6) of R and c V and the formula γ = R c V + 1 we immediately find ∂ t R + u · ∇R = 0, ∂ t c V + u · ∇c V = 0, Note that due to Equation (A4), we can also write the derived PDEs in the divergent form ∂ t (ρϕ) + div(ρϕu) = 0 for ϕ = C α , R, c V .
Equations (A6) and (A7) are well-known in the case of the single-component gas dynamics PDEs. In the general case, they are known too, and we have included their short derivations for completeness and the reader's convenience. The derived PDEs have an important corollary. The original system of PDEs (A1)-(A3) consists of K + n + 1 PDEs for K + n + 1 sought scalar functions ρ, u and θ. Now we see that the closed systems of PDEs exist for n + 4 sought scalar functions ρ, R, c V , u and θ, or even n + 3 sought scalar functions ρ, γ, u and ε, since p = Rρθ = (γ − 1)ρε and E = 0.5ρ|u| 2 + c V ρθ = 0.5ρ|u| 2 + ρε. Other choices of sought functions are possible as well.
In the Euler balance equation for the total energy of the mixture (A3), we first replace (E + p)u → (E + p)u + τ∂ t (E + p)u = (E + p)u + τ (∂ t E + ∂ t p)u + (E + p)∂ t u .
Notice that not only the termsΠ τ andq τ but also the term div (E + p)(u − w 1 ) are different from their respective counterparts in Equations (2) and (3). One can also add the above Navier-Stokes-Fourier terms Π NS and q F , see (10), toΠ τ andq τ and thus obtain the regularized QGD-type system. Such a system, where w 1α is replaced with w 1 , has recently been applied in [39] (without any its derivation, only by the formal analogy to the single-component gas case) and, in a sense, looks closer to the standard form of the corresponding single-component gas system of Equations (23)- (25), than the regularized system (1)-(3) from Section 2, though both the systems for mixtures are reduced to it in the case of a single component. Recall that the derivation of system (1)-(3) was accomplished in a quite different way by aggregating the regularized multi-velocity and multi-temperature gas mixture equations [30].
However, a significant theoretical drawback of the system derived in this Appendix and used in practice, with the regularizing velocity either w 1α , or w 1 instead, is the failure of attempts to prove the non-negativity of the entropy production for it, in contrast to the regularized systems considered in [34,35,37] and system (1)-(3) with the combined regularizing velocity (9). The reason of this drawback is any of the above mentioned differences in the regularized terms.
The expression M 2 = u 2 γ(γ−1)ε that has just arisen in relations (A10) is the squared Mach number in mixtures corresponding to the sound speed in mixtures defined as in Section 2. We have considered the stationary shock waves but it is well-known that non-stationary ones can be reduced to them by passing to a moving system of coordinates [60].