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

Abstract

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 L²-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.

1. 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,2,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 $\tau >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 $\tau$ 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 well-known Coleman–Noll procedure and also ensuring the entropy correctness of the QHD model. The regularizing velocity ${\mathbf{w}}_{\alpha }$ of the component $\alpha$, 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_{\alpha }$ and the total pressure p, respectively. An additional full or partial averaging of ${\mathbf{w}}_{\alpha }$ 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 hyperbolic–parabolic 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 ${\mathbf{w}}_{\alpha }$ 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 ${\mathbf{w}}_{\alpha }$ corresponds to the combined ${\mathbf{w}}_{\alpha }$, and the first sufficiently successful experiments with corresponding difference schemes are presented in [42]. An incomplete attempt to derive the combined ${\mathbf{w}}_{\alpha }$ 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 ${\mathbf{w}}_{\alpha }$ 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 single-component gas dynamics PDEs, see, in particular, reviews: Tadmor, E., Entropy stable schemes ([6], Chapter 17) and Carpenter, M.H.; Fisher, T.C.; Nielsen, E.J. et al. Entropy stable summation-by-parts formulations for compressible computational fluid dynamics ([6], Chapter 19) and references therein.

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 finite-difference counterpart of Proposition 2 is given.

2. 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:

$$\begin{array}{c}{\partial }_t{\rho }_{\alpha }+div\left({\rho }_{\alpha }(\mathbf{u}-{\mathbf{w}}_{\ell \alpha })+{\mathbf{d}}_{\alpha }\right)=0,\alpha =\overline{1,K},\end{array}$$

(1)

$$\begin{array}{c}{\partial }_t\left(\rho \mathbf{u}\right)+div\left(\rho (\mathbf{u}-{\mathbf{w}}_{\ell })\otimes \mathbf{u}\right)+\nabla p=div{\mathrm{\Pi }}_{\ell }+\left(\rho -\ell \tau div\left(\rho \mathbf{u}\right)\right)\mathbf{f},\end{array}$$

(2)

$$\begin{array}{c}{\partial }_{t}E+div\left({0.5\rho \left|\mathbf{u}\right|}^2(\mathbf{u}-{\mathbf{w}}_{\ell })+\langle {\rho }_{\alpha }{h}_{\alpha }(\mathbf{u}-{\mathbf{w}}_{\ell \alpha })\rangle \right)=div(-\mathbf{q}+{\mathrm{\Pi }}_{\ell }\mathbf{u})+\rho (\mathbf{u}-{\mathbf{w}}_{\ell })·\mathbf{f}+Q.\end{array}$$

(3)

Here, the main sought functions are the densities of the mixture components ${\rho }_1>0,\dots ,{\rho }_K>0$ ($K⩾2$ is their amount), their common velocity $\mathbf{u}=(u_1,\dots ,u_n)$ and absolute temperature $\theta >0$. These functions depend on $x=(x_1,\dots ,x_n)\in \mathsf{\Omega }$ and $t⩾0$, where $\mathsf{\Omega }$ is a domain in ${\mathbb{R}}^n$, $n=1,2,3$, and $\alpha =\overline{1,K}$ means that $\alpha =1,\dots ,K$. Vector-functions are written in bold. The operators $div=\nabla ·$ and $\nabla =({\partial }_1,\dots ,{\partial }_n)$ are taken in x, ${\partial }_t=\frac{\partial }{\partial t}$ and ${\partial }_i=\frac{\partial }{\partial 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 $\langle ·\rangle$ is the operation of summation over index $\alpha =\overline{1,K}$.

This regularized system of PDEs was derived in [34] for $K=2$ and ${\mathbf{d}}_1={\mathbf{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 $\alpha$:

$$\begin{array}{c}{p}_{\alpha }=({\gamma }_{\alpha }-1){\rho }_{\alpha }{\epsilon }_{\alpha }=R_{\alpha }{\rho }_{\alpha }\theta ,{\epsilon }_{\alpha }=c_{V\alpha }\theta ,E_{\alpha }=0.5{\rho }_{\alpha }{\left|\mathbf{u}\right|}^2+{\rho }_{\alpha }{\epsilon }_{\alpha },h_{\alpha }={\epsilon }_{\alpha }+\frac{p_{\alpha }}{{\rho }_{\alpha }}=c_{p\alpha }\theta ,\end{array}$$

(4)

with physical constants ${\gamma }_{\alpha }=\frac{R_{\alpha }}{c_{V\alpha }}+1>1$, $R_{\alpha }>0$, $c_{V\alpha }>0$ and $c_{p\alpha }=c_{V\alpha }+R_{\alpha }={\gamma }_{\alpha }{c}_{V\alpha }$, the last two of which are the specific heat capacities at constant volume and pressure, $\alpha =\overline{1,K}$. One can consider any two of four constants ${\gamma }_{\alpha }>1$, $R_{\alpha }>0$, $c_{V\alpha }>0$ and $c_{p\alpha }>0$ as the main independent ones; below, in computations in Section 5, such role is played by ${\gamma }_{\alpha }$ and $c_{V\alpha }$.

The total density and pressure, average specific internal energy and total energy of the mixture are expressed by the formulas

$$\begin{array}{c}\rho =\langle {\rho }_{\alpha }\rangle ,p=\langle p_{\alpha }\rangle =R\rho \theta ,\epsilon =〈\frac{{\rho }_{\alpha }}{\rho }{\epsilon }_{\alpha }〉=c_V\theta ,E=\langle E_{\alpha }\rangle =0.5\rho {\left|\mathbf{u}\right|}^2+\rho \epsilon ,\end{array}$$

(5)

with the average gas mixture parameters

$$\begin{array}{c}R:=〈\frac{{\rho }_{\alpha }}{\rho }{R}_{\alpha }〉,c_V:=〈\frac{{\rho }_{\alpha }}{\rho }{c}_{V\alpha }〉.\end{array}$$

(6)

The second Formula (5) is the Dalton law for mixtures. The function $\frac{{\rho }_{\alpha }}{\rho }=:C_{\alpha }$ is the mass concentration of the mixture component $\alpha$. Consequently, the important formula of the standard form for the total pressure holds as well

$$p=(\gamma -1)\rho \epsilon ,\gamma :=\frac{R}{c_V}+1=\frac{c_p}{c_V}\mathrm{with}{c}_p:=〈\frac{{\rho }_{\alpha }}{\rho }{c}_{p\alpha }〉.$$

In contrast to the single-component case, R, $c_V$ and $\gamma$ are functions, not constants, except for the particular cases $R_1=\dots =R_K$, $c_{V1}=\dots =c_{VK}$ and ${\gamma }_1=\dots {\gamma }_K$, respectively.

In the above PDEs, the following regularizing velocities for the component $\alpha$ and average ones were originally used

$$\begin{array}{c}{\mathbf{w}}_{\ell \alpha }=\ell \frac{\tau }{{\rho }_{\alpha }}div\left({\rho }_{\alpha }\mathbf{u}\right)\mathbf{u}+{\widehat{\mathbf{w}}}_{\alpha },{\widehat{\mathbf{w}}}_{\alpha }=\tau \left((\mathbf{u}·\nabla )\mathbf{u}+\frac{1}{{\rho }_{\alpha }}\nabla p_{\alpha }-\mathbf{f}\right),\end{array}$$

(7)

$$\begin{array}{c}{\mathbf{w}}_{\ell }:=〈\frac{{\rho }_{\alpha }}{\rho }{\mathbf{w}}_{\ell \alpha }〉=\ell \frac{\tau }{\rho }div\left(\rho \mathbf{u}\right)\mathbf{u}+\widehat{\mathbf{w}},\widehat{\mathbf{w}}:=〈\frac{{\rho }_{\alpha }}{\rho }{\widehat{\mathbf{w}}}_{\alpha }〉=\tau \left((\mathbf{u}·\nabla )\mathbf{u}+\frac{1}{\rho }\nabla p-\mathbf{f}\right),\end{array}$$

(8)

see [34], where $\tau =\tau (\rho ,\mathbf{u},\theta )>0$ is a regularization (relaxation) parameter which is usually a function, not constant, with $\rho :=({\rho }_1,\dots ,{\rho }_K)$. Here, $\ell =0$ or 1, and the regularization is of the so-called quasi-gasdynamic (QGD) type for $\ell =1$ or essentially simpler quasi-hydrodynamic (QHD) type for $\ell =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 ${\mathbf{w}}_{\ell \alpha }$ and ${\widehat{\mathbf{w}}}_{\alpha }$, in other words, the full and partial averaging of ${\mathbf{w}}_{\ell \alpha }$.

In this paper, we replace ${\widehat{\mathbf{w}}}_{\alpha }$ by $\widehat{\mathbf{w}}$ making ${\mathbf{w}}_{\ell \alpha }$ combined and dependent on the total pressure p instead of the partial pressure $p_{\alpha }$:

$$\begin{array}{c}{\mathbf{w}}_{\ell \alpha }=\ell \frac{\tau }{{\rho }_{\alpha }}div\left({\rho }_{\alpha }\mathbf{u}\right)\mathbf{u}+\widehat{\mathbf{w}},\widehat{\mathbf{w}}=\tau \left((\mathbf{u}·\nabla )\mathbf{u}+\frac{1}{\rho }\nabla p-\mathbf{f}\right)\end{array}$$

(9)

and analyze the effect of this replacement discussed above in Introduction. Notice that the replacement does not affect the validity of Formula (8) for ${\mathbf{w}}_{\ell }$.

The viscosity tensor and heat flux are expressed, respectively, by the formulas

$${\mathrm{\Pi }}_{\ell }={\mathrm{\Pi }}^{NS}+{\mathrm{\Pi }}_{\ell }^{\tau },\mathbf{q}={\mathbf{q}}^F+{\mathbf{q}}^d+\ell {\mathbf{q}}^{\tau }$$

and contain the standard-type terms and the regularizing ones with the superscript $\tau$. The classical Navier–Stokes viscosity tensor and the Fourier heat flux are given by the formulas

$$\begin{array}{c}{\mathrm{\Pi }}^{NS}=\mu \left(\nabla \mathbf{u}+{(\nabla \mathbf{u})}^T\right)+\left(\lambda -\frac{2}{3}\mu \right)(div\mathbf{u})\mathbb{I},-{\mathbf{q}}^F=\varkappa \nabla \theta ,\end{array}$$

(10)

where $\mu >0$, $\lambda ⩾0$ and $\varkappa >0$ are the total coefficients of dynamic and bulk viscosities and heat conductivity (which can depend on the sought functions $(\rho ,\mathbf{u},\theta )$), $\nabla \mathbf{u}={\left\{{\partial }_i{u}_j\right\}}_{i,j=1}^n$ and $\mathbb{I}$ is the unit tensor of order n.

Next, the regularizing viscosity tensor and heat flux are given by the formulas

$$\begin{array}{c}{\mathrm{\Pi }}_{\ell }^{\tau }=\rho \mathbf{u}\otimes \widehat{\mathbf{w}}+\ell \tau \left(\mathbf{u}·\nabla p+\langle {\gamma }_{\alpha }{p}_{\alpha }\rangle div\mathbf{u}-\langle {\gamma }_{\alpha }{Q}_{\alpha }\rangle +Q\right)\mathbb{I},\end{array}$$

(11)

$$\begin{array}{c}-{\mathbf{q}}^{\tau }=\tau \left\{\left(c_V\rho \nabla \theta -\theta \nabla \left(R\rho \right)\right)·\mathbf{u}-Q\right\}\mathbf{u}.\end{array}$$

(12)

The density of body force $\mathbf{f}$ and intensities of heat sources $Q_{\alpha }⩾0$ (acting on the component $\alpha$) are given functions, and $Q:=\langle Q_{\alpha }\rangle ⩾0$.

Finally, we consider the diffusion fluxes and additional respective heat flux of the form

$$\begin{array}{c}-{\mathbf{d}}_{\alpha }:={〈d_{\alpha \beta }\left(\nabla (G_{\alpha }-G_{\beta })+(e_{\alpha }-e_{\beta })\nabla \theta \right)〉}_{\beta }={〈d_{\alpha \beta }\left(\nabla G_{\alpha }+e_{\alpha }\nabla \theta -(\nabla G_{\beta }+e_{\beta }\nabla \theta )\right)〉}_{\beta },\end{array}$$

(13)

$$\begin{array}{c}{\mathbf{q}}^d=〈\left(G_{\alpha }+e_{\alpha }\theta \right){\mathbf{d}}_{\alpha }〉,\end{array}$$

(14)

$$\begin{array}{c}{G}_{\alpha }:={\epsilon }_{\alpha }+\frac{p_{\alpha }}{{\rho }_{\alpha }}-s_{\alpha }\theta =h_{\alpha }-s_{\alpha }\theta =(c_{p\alpha }-s_{\alpha })\theta ,s_{\alpha }=s_{\alpha 0}-R_{\alpha }ln\frac{{\rho }_{\alpha }}{{\rho }_{\alpha 0}}+c_{V\alpha }ln\frac{\theta }{{\theta }_0},\end{array}$$

(15)

where ${\langle ·\rangle }_{\beta }$ means the summation over index $\beta =\overline{1,K}$. The functions $G_{\alpha }$ and $s_{\alpha }$ are the usual Gibbs potential and specific entropy of the component $\alpha$, and $s_{\alpha 0}$, ${\rho }_{\alpha 0}>0$ and ${\theta }_0>0$ are constant reference values for $s_{\alpha }$, ${\rho }_{\alpha }$ and $\theta$, $\alpha =\overline{1,K}$. The functions-coefficients $d_{\alpha \beta }$ and $e_{\alpha }$ can depend on the sought functions. Their specific form is not essential below, and we only assume the symmetry property $d_{\alpha \beta }=d_{\beta \alpha }$ for any $\alpha \ne \beta$.

Let ${\langle ·\rangle }_{\alpha ,\beta }:=\langle {\langle ·\rangle }_{\beta }\rangle$ mean the summation over $\alpha ,\beta =\overline{1,K}$. Using permutations of indices $\alpha$ and $\beta$ and then this symmetry property, we obtain two identities

$$\begin{array}{c}{\langle d_{\alpha \beta }({\phi }_{\alpha }-{\phi }_{\beta })\rangle }_{\alpha ,\beta }={\langle (d_{\alpha \beta }-d_{\beta \alpha }){\phi }_{\alpha }\rangle }_{\alpha ,\beta }=0,\\ {\langle d_{\alpha \beta }({\phi }_{\alpha }-{\phi }_{\beta }){\psi }_{\alpha }\rangle }_{\alpha ,\beta }\\ =0.5{\langle d_{\alpha \beta }({\phi }_{\alpha }-{\phi }_{\beta }){\phi }_{\alpha }+d_{\beta \alpha }({\phi }_{\beta }-{\phi }_{\alpha }){\psi }_{\beta }\rangle }_{\alpha ,\beta }=0.5{〈d_{\alpha \beta }({\phi }_{\alpha }-{\phi }_{\beta })({\psi }_{\alpha }-{\psi }_{\beta })〉}_{\alpha ,\beta }\end{array}$$

(16)

for any numbers ${\phi }_1,\dots ,{\phi }_K$ and ${\psi }_1,\dots ,{\psi }_K$. The first identity implies the important physical property $\langle {\mathbf{d}}_{\alpha }\rangle =0$, and the second one will also be essential below.

We can avoid the explicit usage of $s_{\alpha }$ due to the formulas

$$\begin{array}{c}\nabla G_{\alpha }+e_{\alpha }\nabla \theta =-\theta \nabla s_{\alpha }+(c_{p\alpha }-s_{\alpha }+e_{\alpha })\nabla \theta =R_{\alpha }\frac{\theta }{{\rho }_{\alpha }}\nabla {\rho }_{\alpha }+{\overline{e}}_{\alpha }\nabla \theta \\ =\frac{1}{{\rho }_{\alpha }}\nabla p_{\alpha }+({\overline{e}}_{\alpha }-R_{\alpha })\nabla \theta \mathrm{and}{G}_{\alpha }+e_{\alpha }\theta =(c_{V\alpha }+{\overline{e}}_{\alpha })\theta ,\mathrm{with}{\overline{e}}_{\alpha }:=R_{\alpha }-s_{\alpha }+e_{\alpha }.\end{array}$$

(17)

In particular, for $e_{\alpha }=s_{\alpha }$, we obtain the simplest formulas

$$\nabla G_{\alpha }+s_{\alpha }\nabla \theta =R_{\alpha }\frac{\theta }{{\rho }_{\alpha }}\nabla {\rho }_{\alpha }+R_{\alpha }\nabla \theta =\frac{1}{{\rho }_{\alpha }}\nabla p_{\alpha },G_{\alpha }+s_{\alpha }\theta =h_{\alpha }.$$

The general multicomponent case ($K⩾2$) for $\ell =0,1$ in the presence of ${\mathbf{d}}_{\alpha }$ and ${\mathbf{q}}^d$ has recently been studied mathematically in [35] but only in the particular case $d_{\alpha \beta }\equiv d_0$ and $K{e}_{\alpha }-{\langle e_{\beta }\rangle }_{\beta }=b_{\alpha }$, that is the same for $K=2$ but much less general for $K⩾3$. The above quantities ${\mathbf{d}}_{\alpha }$ and ${\mathbf{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].

Without the regularization, i.e., for $\tau =0$, the above regularized system of PDEs is simplified and reduced to the compressible Navier–Stokes–Fourier-type system for the one-velocity and one-temperature multicomponent gas mixture dynamics for ${\mu }_{\alpha }>0$, ${\lambda }_{\alpha }>0$ and ${\varkappa }_{\alpha }>0$ or the Euler-type one for ${\mu }_{\alpha }={\lambda }_{\alpha }={\varkappa }_{\alpha }=0$, $\alpha =\overline{1,K}$, in particular, see [1,2,50] and references therein, and also Appendix A.

In [34], the above total coefficients $\mu$, $\lambda$ and $\varkappa$ are defined simply as

$$\begin{array}{c}\mu =\langle {\mu }_{\alpha }\rangle ,\lambda =\langle {\lambda }_{\alpha }\rangle ,\varkappa =\langle {\varkappa }_{\alpha }\rangle ,\end{array}$$

(18)

i.e., the sums of the corresponding coefficients of the components. These coefficients can be artificial depending on $\tau$ 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 $\tau$ and them are as follows

$$\begin{array}{c}\tau =\frac{ah}{c_s+i_{\tau }\left|\mathbf{u}\right|},\mu =\tau \langle a_{S\alpha }{p}_{\alpha }\rangle ,\lambda =\tau \langle a_{1S\alpha }{p}_{\alpha }\rangle ,\varkappa =\tau \langle a_{Pr\alpha }{\gamma }_{\alpha }{c}_{V\alpha }{p}_{\alpha }\rangle \end{array}$$

(19)

in accordance with Formula (18). Here, $0<a⩽1$ is a parameter, $a_{S\alpha }>0$ and $a_{Pr\alpha }>0$ are the Schmidt and inverse Prandtl numbers for the component $\alpha$; $a_{1S\alpha }⩾0$ is a counterpart of $a_{S\alpha }$ (in particular, $a_{1S\alpha }=0$), which can be also used as adjusting numerical parameters,

$$c_s=\sqrt{\gamma (\gamma -1)\epsilon }=\sqrt{\gamma R\theta }$$

is the sound speed of the mixture, $i_{\tau }=0$ or 1, and h is a characteristic size of the spatial mesh. In the case of $a_{S\alpha }=a_S$, $a_{1S\alpha }=a_{1S}$ and $a_{Pr\alpha }=a_{Pr}$ independent of $\alpha$, the formulas for $\mu$, $\lambda$ and $\varkappa$ are simplified:

$$\begin{array}{c}\mu =a_S\tau p,\lambda =a_{1S}\tau p,\varkappa =\tau a_{Pr}\langle {\gamma }_{\alpha }{c}_{V\alpha }{p}_{\alpha }\rangle .\end{array}$$

For the single-component gas dynamics, see such formulas, in particular, in [8,9,19].

Recall that $R_{\alpha }=\frac{R_0}{m_{\alpha }}$, where $R_0$ is the universal gas constant and $m_{\alpha }>0$ is the molecular mass of gas $\alpha$. In some cases, ${\gamma }_{\alpha }$ and $m_{\alpha }$ are taken as the two main gas constants, and the average molecular mass of the mixture m is defined by $\frac{1}{m}=〈\frac{{\rho }_{\alpha }}{\rho }\frac{1}{m_{\alpha }}〉.$ Then the other gas constants can be expressed in the form

$$c_{V\alpha }=\frac{R_0}{({\gamma }_{\alpha }-1)m_{\alpha }},R=\frac{R_0}{m},c_V=R_0〈\frac{{\rho }_{\alpha }}{\rho }\frac{1}{({\gamma }_{\alpha }-1)m_{\alpha }}〉,\gamma -1=\frac{R}{c_V}={〈\frac{{\rho }_{\alpha }}{\rho }\frac{1}{{\gamma }_{\alpha }-1}\frac{m}{m_{\alpha }}〉}^{-1}.$$

We first give some inequalities for R, $c_V$, m and $\gamma$.

 Proposition 1.

1. The two-sided bounds hold

$$\underset{\alpha =\overline{1,K}}{min}{R}_{\alpha }⩽R⩽\underset{\alpha =\overline{1,K}}{max}{R}_{\alpha },\underset{\alpha =\overline{1,K}}{min}{c}_{V\alpha }⩽c_V⩽\underset{\alpha =\overline{1,K}}{max}{c}_{V\alpha },\underset{\alpha =\overline{1,K}}{min}{m}_{\alpha }⩽m⩽\underset{\alpha =\overline{1,K}}{max}{m}_{\alpha }.$$

Moreover, all the inequalities are strict, except for the particular cases $R_1=\dots =R_K$, $c_{V1}=\dots =c_{VK}$ and $m_1=\dots =m_K$, respectively.

2. The formula and two-sided bounds also hold

$$\begin{array}{c}\underset{\alpha =\overline{1,K}}{min}{\gamma }_{\alpha }⩽\gamma =\frac{R+c_V}{c_V}=\frac{\langle {\rho }_{\alpha }{c}_{V\alpha }{\gamma }_{\alpha }\rangle }{\langle {\rho }_{\alpha }{c}_{V\alpha }\rangle }⩽\underset{\alpha =\overline{1,K}}{max}{\gamma }_{\alpha }.\end{array}$$

(20)

Moreover, both the bounds are strict excluding the particular case ${\gamma }_1=\dots ={\gamma }_K$; in that case, we have $\gamma ={\gamma }_1$ even if $c_{V\alpha }$ are not identical for all $\alpha =\overline{1,K}$.

3. The following relations hold

$$\begin{array}{c}\langle {\gamma }_{\alpha }{p}_{\alpha }\rangle =\tilde{\gamma }p⩾\gamma p\mathit{with}\tilde{\gamma }:=\frac{\langle {\rho }_{\alpha }{R}_{\alpha }{\gamma }_{\alpha }\rangle }{\langle {\rho }_{\alpha }{R}_{\alpha }\rangle }=〈\frac{{\rho }_{\alpha }}{\rho }\frac{m}{m_{\alpha }}{\gamma }_{\alpha }〉\iff 〈\frac{{\rho }_{\alpha }}{\rho }{c}_{s\alpha }^2〉⩾c_s^2=\gamma (\gamma -1)\epsilon .\end{array}$$

(21)

with $c_{s\alpha }^2={\gamma }_{\alpha }({\gamma }_{\alpha }-1){\epsilon }_{\alpha }$. Here, $\langle {\gamma }_{\alpha }{p}_{\alpha }\rangle =\gamma p$, or $\tilde{\gamma }=\gamma$, in the case ${\gamma }_1=\dots ={\gamma }_K$ only.

 Proof. 

Items 1 and 2 are elementary. Item 3 is valid since

$$\begin{array}{c}\langle {\gamma }_{\alpha }{p}_{\alpha }\rangle -\gamma p=\langle ({\gamma }_{\alpha }-1)p_{\alpha }\rangle -(\gamma -1)p=\left(〈\frac{R_{\alpha }^2}{c_{V\alpha }}{\rho }_{\alpha }〉-\frac{{\langle {\rho }_{\alpha }{R}_{\alpha }\rangle }^2}{\langle {\rho }_{\alpha }{c}_{V\alpha }\rangle }\right)\theta ⩾0\end{array}$$

(22)

owing to the Cauchy inequality

$${\langle {\rho }_{\alpha }{R}_{\alpha }\rangle }^2={〈R_{\alpha }\sqrt{\frac{{\rho }_{\alpha }}{c_{V\alpha }}}\sqrt{{\rho }_{\alpha }{c}_{V\alpha }}〉}^2⩽〈\frac{R_{\alpha }^2}{c_{V\alpha }}{\rho }_{\alpha }〉\langle c_{V\alpha }{\rho }_{\alpha }\rangle .$$

The inequality becomes an equality only for $\frac{R_{\alpha }^2}{c_{V\alpha }^2}$ independent of $\alpha$, i.e., ${\gamma }_1=\dots ={\gamma }_K$. □

 Remark 1.

Starting from Formula (22), we can accomplish the following transformations

$$\begin{array}{c}\left(\langle {\gamma }_{\alpha }{p}_{\alpha }\rangle -\gamma p\right)\frac{\langle {\rho }_{\alpha }{c}_{V\alpha }\rangle }{\theta }=〈{({\gamma }_{\alpha }-1)}^2{c}_{V\alpha }{\rho }_{\alpha }〉\langle {\rho }_{\alpha }{c}_{V\alpha }\rangle -{〈({\gamma }_{\alpha }-1){\rho }_{\alpha }{c}_{V\alpha }〉}^2\\ ={〈\left({({\gamma }_{\alpha }-1)}^2-({\gamma }_{\alpha }-1)({\gamma }_{\beta }-1)\right){\rho }_{\alpha }{\rho }_{\beta }{c}_{V\alpha }{c}_{V\beta }〉}_{\alpha ,\beta }={〈({\gamma }_{\alpha }-1)({\gamma }_{\alpha }-{\gamma }_{\beta }){\rho }_{\alpha }{\rho }_{\beta }{c}_{V\alpha }{c}_{V\beta }〉}_{\alpha ,\beta }.\end{array}$$

Permuting indexes α and β, similarly to identity (16), we derive the representation

$$\begin{array}{c}\langle {\gamma }_{\alpha }{p}_{\alpha }\rangle -\gamma p=\frac{\theta }{2\langle c_{V\alpha }{\rho }_{\alpha }\rangle }{〈{({\gamma }_{\alpha }-{\gamma }_{\beta })}^2{\rho }_{\alpha }{\rho }_{\beta }{c}_{V\alpha }{c}_{V\beta }〉}_{\alpha ,\beta }⩾0\end{array}$$

since $({\gamma }_{\alpha }-1)({\gamma }_{\alpha }-{\gamma }_{\beta })+({\gamma }_{\beta }-1)({\gamma }_{\beta }-{\gamma }_{\alpha })={({\gamma }_{\alpha }-{\gamma }_{\beta })}^2$. This formula also implies Item 3.

Formulas for $\gamma$ and $\tilde{\gamma }$ in relations (20) and (21) show that both of them are averages of ${\gamma }_1,\dots ,{\gamma }_K$, and clearly $c_{V\alpha }$ and $R_{\alpha }$ can be scaled in these formulas. For $\tilde{\gamma }$, the same bounds as in Item 2 for $\gamma$ are valid. Note that, in our expression for ${\mathrm{\Pi }}^{\tau }$, see (11), we use the term $\langle {\gamma }_{\alpha }{p}_{\alpha }\rangle$ as in [34,35], in contrast to $\gamma 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).

We first consider the sought functions of the particular form.

 Proposition 2.

Let $0<C_{\alpha }<1$, $\alpha =\overline{1,K}$, be arbitrary constants such that $\langle C_{\alpha }\rangle =1$. Consider the sought functions of the particular form ${\rho }_{\alpha }=C_{\alpha }\rho$ with $\rho >0$ ($\alpha =\overline{1,K}$), $\mathbf{u}$ and $\theta >0$ and the case of ${\mathbf{d}}_{\alpha }=0$ and $\psi =\psi (\rho ,\mathbf{u},\theta )$ for the functions $\psi =\tau ,\mu ,\lambda ,\varkappa$. 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

$$\begin{array}{c}{\partial }_t\rho +div\left(\rho (\mathbf{u}-{\mathbf{w}}_{\ell })\right)=0,\end{array}$$

(23)

$$\begin{array}{c}{\partial }_t\left(\rho \mathbf{u}\right)+div\left(\rho (\mathbf{u}-{\mathbf{w}}_{\ell })\otimes \mathbf{u}\right)+\nabla p=div{\mathrm{\Pi }}_{\ell }+\left(\rho -\ell \tau div\left(\rho \mathbf{u}\right)\right)\mathbf{f},\end{array}$$

(24)

$$\begin{array}{c}{\partial }_{t}E+div\left((E+p)(\mathbf{u}-{\mathbf{w}}_{\ell })\right)=div(-\mathbf{q}+{\mathrm{\Pi }}_{\ell }\mathbf{u})+\rho (\mathbf{u}-{\mathbf{w}}_{\ell })·\mathbf{f}+Q\end{array}$$

(25)

for the sought functions ρ, $\mathbf{u}$ and θ and $(x,t)\in \mathsf{\Omega }\times [0,T]$. Here,

$$\begin{array}{c}p=(\gamma -1)\rho \epsilon =R\rho \theta ,\epsilon =c_V\theta ,E=0.5\rho {\left|\mathbf{u}\right|}^2+\rho \epsilon ,R=\langle C_{\alpha }{R}_{\alpha }\rangle ,c_V=\langle C_{\alpha }{c}_{V\alpha }\rangle \end{array}$$

(26)

with constant $\gamma =\frac{R}{c_V}+1$, R and $c_V$, together with

$$\begin{array}{c}{\mathbf{w}}_{\ell }=\ell \frac{\tau }{\rho }div\left(\rho \mathbf{u}\right)\mathbf{u}+\widehat{\mathbf{w}},\widehat{\mathbf{w}}=\tau \left((\mathbf{u}·\nabla )\mathbf{u}+\frac{1}{\rho }\nabla p-\mathbf{f}\right),\\ {\mathrm{\Pi }}_{\ell }={\mathrm{\Pi }}^{NS}+{\mathrm{\Pi }}_{\ell }^{\tau },{\mathrm{\Pi }}_{\ell }^{\tau }=\rho \mathbf{u}\otimes \widehat{\mathbf{w}}+\ell \tau \left(\mathbf{u}·\nabla p+\tilde{\gamma }pdiv\mathbf{u}-({\gamma }_1-1)Q\right)\mathbb{I},\end{array}$$

(27)

$$\begin{array}{c}\mathbf{q}={\mathbf{q}}^F+\ell {\mathbf{q}}^{\tau },-{\mathbf{q}}^{\tau }=\tau \left\{\left(c_V\rho \nabla \theta -R\theta \nabla \rho \right)·\mathbf{u}-Q\right\}\mathbf{u}=\tau \left\{\left(\rho \nabla \epsilon -\frac{p}{\rho }\nabla \rho \right)·\mathbf{u}-Q\right\}\mathbf{u},\end{array}$$

(28)

where $\tilde{\gamma }=\langle C_{\alpha }{R}_{\alpha }{\gamma }_{\alpha }\rangle /\langle C_{\alpha }{R}_{\alpha }\rangle$ and the above formulas (10) for ${\mathrm{\Pi }}^{NS}$ and ${\mathbf{q}}^F$ are in use.

 Proof. 

For ${\rho }_{\alpha }=C_{\alpha }\rho$, under the assumptions made about $C_{\alpha }$, we clearly obtain $\rho =\langle {\rho }_{\alpha }\rangle$ and

$${\mathbf{w}}_{\ell \alpha }=\ell \frac{\tau }{\rho }div\left(\rho \mathbf{u}\right)\mathbf{u}+\widehat{\mathbf{w}}={\mathbf{w}}_{\ell },\langle {\rho }_{\alpha }{h}_{\alpha }\rangle =\rho \langle C_{\alpha }{c}_{p\alpha }\rangle \theta =\rho \epsilon +p.$$

Thus, all the balance PDEs for the mass of components (1), after division by $C_{\alpha }$, are reduced to Equation (23). Moreover, expressions (11) for ${\mathrm{\Pi }}_{\ell }^{\tau }$ and (12) for $-{\mathbf{q}}^{\tau }$ are reduced to those given in Formulas (27) and (28), using the expression for $\tilde{\gamma }$ 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 ${\gamma }_{\alpha }>1$ and $c_{V\alpha }>0$, $\alpha =\overline{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 $\tilde{\gamma }=\gamma$ in Formula (27) only in the particular case ${\gamma }_1=\dots ={\gamma }_K$, see Proposition 1, Item 3 or Remark 1.

Applying the operation $\langle ·\rangle$ to the mass balance equation for the mixture components (1), using the formula $\langle {\rho }_{\alpha }(\mathbf{u}-{\mathbf{w}}_{\ell \alpha })\rangle =\rho (\mathbf{u}-{\mathbf{w}}_{\ell })$ valid according to the first expression (8), and the property $\langle {\mathbf{d}}_{\alpha }\rangle =0$, we obtain the important total mass balance equation

$${\partial }_t\rho +div\left(\rho (\mathbf{u}-{\mathbf{w}}_{\ell })\right)=0.$$

(29)

Here, ${\rho }_1,\dots ,{\rho }_K$ and $\theta$ appear only implicitly since $\rho =\langle {\rho }_{\alpha }\rangle$ and p in ${\mathbf{w}}_{\ell }$ 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

$$\begin{array}{c}0.5{\partial }_t{\left(\rho \right|\mathbf{u}|}^2)+0.5div\left({\rho \left|\mathbf{u}\right|}^2(\mathbf{u}-{\mathbf{w}}_{\ell })\right)+(\nabla p)·\mathbf{u}=(div\mathrm{\Pi })·\mathbf{u}+\left(\rho -\ell \tau div\left(\rho \mathbf{u}\right)\right)\mathbf{f}·\mathbf{u},\end{array}$$

(30)

$$\begin{array}{c}{\partial }_t\left(\rho \epsilon \right)+div\langle {\rho }_{\alpha }{\epsilon }_{\alpha }(\mathbf{u}-{\mathbf{w}}_{\ell \alpha })\rangle +pdiv\mathbf{u}=div(-\mathbf{q}+\langle p_{\alpha }{\mathbf{w}}_{\ell \alpha }\rangle )+\mathrm{\Pi }:\nabla \mathbf{u}-\rho \widehat{\mathbf{w}}·\mathbf{f}+Q,\end{array}$$

(31)

where the symbol: denotes the scalar product of tensors. The derivation exploits the total mass balance equation (29) and is valid for any ${\mathbf{w}}_{\ell }$, for the former equation, and exploits only the relation ${\mathbf{w}}_{\ell }=\ell \tau div\left(\rho \mathbf{u}\right)\mathbf{u}+\widehat{\mathbf{w}}$ (for $\mathbf{f}\neg \equiv 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=〈\frac{{\rho }_{\alpha }}{\rho }{s}_{\alpha }〉$. 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 ${\mathbf{d}}_{\alpha }$, for $K⩾3$; this form is applicable in [35,37] as well.

 Theorem 1.

Let $d_{\alpha \beta }=d_{\beta \alpha }⩾0$ for any $\alpha \ne \beta$. The following regularized entropy balance equation for the multicomponent mixture in the presence of diffusion fluxes holds

$$\begin{array}{c}{\partial }_t\left(\rho s\right)+div\left\{\langle {\rho }_{\alpha }{s}_{\alpha }(\mathbf{u}-{\mathbf{w}}_{\ell \alpha })\rangle +\langle e_{\alpha }{\mathbf{d}}_{\alpha }\rangle +\frac{1}{\theta }({\mathbf{q}}^F+\ell {\mathbf{q}}^{\tau })\right\}={\mathcal{P}}^{NS}+{\mathcal{P}}^{\tau }\end{array}$$

with the entropy production ${\mathcal{P}}^{NS}+{\mathcal{P}}^{\tau }$, where

$$\begin{array}{c}{\mathcal{P}}^{NS}=\frac{1}{\theta }\{\frac{\mu }{2}|\nabla \mathbf{u}+\nabla {\mathbf{u}}^T{|}_F^2+\left(\lambda -\frac{2}{3}\mu \right){(div\mathbf{u})}^2+\frac{1}{\theta }\varkappa {|\nabla \theta |}^2\\ +\frac{1}{2}{〈d_{\alpha \beta }|\nabla (G_{\alpha }-G_{\beta })+(e_{\alpha }-e_{\beta })\nabla \theta {|}^2〉}_{\alpha ,\beta }\}⩾0,\\ {\mathcal{P}}^{\tau }=\frac{\rho }{\tau \theta }{|\widehat{\mathbf{w}}|}^2+\ell 〈\tau \frac{R_{\alpha }}{{\rho }_{\alpha }}{(div\left({\rho }_{\alpha }\mathbf{u}\right))}^2+\tau c_{V\alpha }{\rho }_{\alpha }{\left(\mathbf{u}·\nabla ln\theta +({\gamma }_{\alpha }-1)div\mathbf{u}-\frac{({\gamma }_{\alpha }-1)Q_{\alpha }}{2p_{\alpha }}\right)}^2〉\\ +〈\frac{Q_{\alpha }}{\theta }\left(1-\ell \frac{\tau ({\gamma }_{\alpha }-1)Q_{\alpha }}{4p_{\alpha }}\right)〉,\end{array}$$

and ${|·|}_F$ is the Frobenius norm. Moreover, ${\mathcal{P}}^{\tau }$ is non-negative for $\ell =0$, as well as for $\ell =1$ under the condition

$$\tau 〈\frac{({\gamma }_{\alpha }-1)Q_{\alpha }^2}{4p_{\alpha }}〉⩽Q.$$

This condition is certainly true provided that $\tau ({\gamma }_{\alpha }-1)Q_{\alpha }⩽4p_{\alpha }$, $\alpha =\overline{1,K}$.

 Proof. 

According to ([35], proof of Theorem 3.1), the following preliminary equation involving the entropies of the mixture and the components holds

$$\begin{array}{c}{\partial }_t\left(\rho s\right)+div\left\{\langle {\rho }_{\alpha }{s}_{\alpha }(\mathbf{u}-{\mathbf{w}}_{\ell \alpha })+e_{\alpha }{\mathbf{d}}_{\alpha }\rangle +\frac{1}{\theta }({\mathbf{q}}^F+\ell {\mathbf{q}}^{\tau })\right\}=-\frac{1}{\theta }\nabla \theta ·\langle e_{\alpha }{\mathbf{d}}_{\alpha }\rangle -\frac{1}{\theta }〈\nabla G_{\alpha }·{\mathbf{d}}_{\alpha }〉\\ +\frac{1}{{\theta }^2}\varkappa {|\nabla \theta |}^2+\frac{1}{\theta }\left\{\frac{\mu }{2}{|\nabla \mathbf{u}+\nabla {\mathbf{u}}^T|}^2+\left(\lambda -\frac{2}{3}\mu \right){(div\mathbf{u})}^2\right\}+\frac{1}{\theta }{B}^{\tau },\end{array}$$

(32)

where

$$B^{\tau }:=\langle \nabla p_{\alpha }·{\mathbf{w}}_{\ell \alpha }\rangle -\ell {\mathbf{q}}^{\tau }·\frac{1}{\theta }\nabla \theta +{\mathrm{\Pi }}^{\tau }:\nabla \mathbf{u}-\rho \widehat{\mathbf{w}}·\mathbf{f}+Q.$$

This equation is derived from the balance PDEs (1) and (31) and does not exploit specific expressions for ${\mathbf{w}}_{\ell \alpha }$ and $\widehat{\mathbf{w}}$ in them. In the first term on the right, we have taken into account the following obvious formula

$$\begin{array}{c}〈\frac{G_{\alpha }}{\theta }{\mathbf{d}}_{\alpha }〉-\frac{{\mathbf{q}}^d}{\theta }=-\langle e_{\alpha }{\mathbf{d}}_{\alpha }\rangle ,\end{array}$$

see definition (14) of ${\mathbf{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) in the form

$$\begin{array}{c}-\frac{1}{\theta }\nabla \theta ·\langle e_{\alpha }{\mathbf{d}}_{\alpha }\rangle -\frac{1}{\theta }〈\nabla G_{\alpha }·{\mathbf{d}}_{\alpha }〉=-\frac{1}{\theta }〈\left(\nabla G_{\alpha }+e_{\alpha }\nabla \theta \right)·{\mathbf{d}}_{\alpha }〉\\ =\frac{1}{2\theta }{〈d_{\alpha \beta }{\left(\nabla G_{\alpha }+e_{\alpha }\nabla \theta -(\nabla G_{\beta }+e_{\beta }\nabla \theta )\right)}^2〉}_{\alpha ,\beta }.\end{array}$$

Next, in the case of expressions (9), we can extract and collect the terms with $\widehat{\mathbf{w}}$ from the first, third and fourth terms of $B^{\tau }$ and thus write

$$\begin{array}{c}{B}^{\tau }=\langle \nabla p_{\alpha }·\widehat{\mathbf{w}}\rangle +(\mathbf{u}·\nabla )\mathbf{u}·\rho \widehat{\mathbf{w}}-\rho \widehat{\mathbf{w}}·\mathbf{f}+{\tilde{B}}^{\tau }\\ =\left(\rho (\mathbf{u}·\nabla )\mathbf{u}+\nabla p-\rho \mathbf{f}\right)·\widehat{\mathbf{w}}+{\tilde{B}}^{\tau }=\frac{\rho }{\tau }{\left|\widehat{\mathbf{w}}\right|}^2+{\tilde{B}}^{\tau }.\end{array}$$

Concerning the remainder ${\tilde{B}}^{\tau }$, the following formula holds $\frac{1}{\theta }{\tilde{B}}^{\tau }={\mathcal{P}}^{\tau }-\frac{\rho }{\tau \theta }{\left|\widehat{\mathbf{w}}\right|}^2,$ see ([35], proof of Theorem 3.1) and the references therein, that completes the proof. □

Clearly, ${\mathcal{P}}^{NS}+\frac{Q}{\theta }$ and ${\mathcal{P}}^{\tau }-\frac{Q}{\theta }$ are the Navier–Stokes–Fourier and regularizing contributions to the entropy production. Theorem 1 remains valid for $\tau ⩾0$ (in particular, $\tau =0$, i.e., without a regularization), when one should pass to a different form for the first relaxation term:

$$\frac{\rho }{\tau \theta }|\widehat{\mathbf{w}}{|}^2=\frac{\tau }{\rho \theta }{|\rho (\mathbf{u}·\nabla )\mathbf{u}+\nabla p-\rho \mathbf{f}|}^2.$$

3. Linearized Regularized System of PDEs for Gas Mixture Dynamics, Its Properties and Corollaries

3.1. An Auxiliary Reduction of the Balance Equations

Let $\mathbf{f}=0$ and $Q_1=\dots =Q_K=0$. We introduce the vector of the sought functions $\mathbf{z}:=(\rho ,\mathbf{u},\theta )$ and first present an important auxiliary reduction of the balance PDEs for ${\rho }_{\alpha }$, $\mathbf{u}$ and $\theta$ up to the terms $O\left(\right|\nabla \mathbf{z}{|}^2)$.

 Lemma 1.

The following reduced PDEs hold: for the densities of the components

$$\begin{array}{c}{\partial }_t{\rho }_{\alpha }+\nabla {\rho }_{\alpha }·\mathbf{u}+{\rho }_{\alpha }div\mathbf{u}\\ =\tau \left[\frac{{\rho }_{\alpha }\theta }{\rho }{\langle R_{\beta }\Delta {\rho }_{\beta }\rangle }_{\beta }+\ell [\mathbf{u}·(\mathbf{u}·\nabla )\nabla ]{\rho }_{\alpha }+(\ell +1){\rho }_{\alpha }(\mathbf{u}·\nabla )div\mathbf{u}+R{\rho }_{\alpha }\Delta \theta \right]\\ +\theta {〈d_{\alpha \beta }\left(\frac{R_{\alpha }}{{\rho }_{\alpha }}\Delta {\rho }_{\alpha }-\frac{R_{\beta }}{{\rho }_{\beta }}\Delta {\rho }_{\beta }\right)〉}_{\beta }+{〈d_{\alpha \beta }({\overline{e}}_{\alpha }-{\overline{e}}_{\beta })〉}_{\beta }\Delta \theta +{O\left(\right|\nabla \mathbf{z}|}^2),\alpha =\overline{1,K},\end{array}$$

(33)

for the velocity

$$\begin{array}{c}{\partial }_t\mathbf{u}+\frac{\theta }{\rho }\langle R_{\alpha }\nabla {\rho }_{\alpha }\rangle +(\mathbf{u}·\nabla )\mathbf{u}+R\nabla \theta \\ =(\ell +1)\tau \frac{\theta }{\rho }(\mathbf{u}·\nabla )\langle R_{\alpha }\nabla {\rho }_{\alpha }\rangle +\frac{\mu }{\rho }\Delta \mathbf{u}+\frac{\chi }{\rho }\nabla div\mathbf{u}+\ell \tau \frac{\langle {\gamma }_{\alpha }{p}_{\alpha }\rangle }{\rho }\nabla div\mathbf{u}\\ +\tau [\mathbf{u}·(\mathbf{u}·\nabla )\nabla ]\mathbf{u}+(\ell +1)\tau R(\mathbf{u}·\nabla )\nabla \theta +O\left(\right|\nabla \mathbf{z}{|}^2),\end{array}$$

(34)

with $\chi :=\frac{1}{3}\mu +\lambda$, and for the temperature

$$\begin{array}{c}{\partial }_t\theta +\frac{R}{c_V}\theta div\mathbf{u}+\mathbf{u}·\nabla \theta \\ =\tau \frac{R{\theta }^2}{c_V\rho }\langle R_{\alpha }\Delta {\rho }_{\alpha }\rangle +(\ell +1)\tau \frac{R\theta }{c_V}(\mathbf{u}·\nabla )div\mathbf{u}+\ell \tau [\mathbf{u}·(\mathbf{u}·\nabla )\nabla ]\theta +\left(\frac{\varkappa }{c_V\rho }+\tau \frac{R^2\theta }{c_V\rho }\right)\Delta \theta \\ +\frac{\theta }{c_V\rho }{〈d_{\alpha \beta }\left(\frac{R_{\alpha }}{{\rho }_{\alpha }}\Delta {\rho }_{\alpha }-\frac{R_{\beta }}{{\rho }_{\beta }}\Delta {\rho }_{\beta }\right){\overline{e}}_{\alpha }〉}_{\alpha ,\beta }+\frac{1}{c_V\rho }{〈d_{\alpha \beta }({\overline{e}}_{\alpha }-{\overline{e}}_{\beta }){\overline{e}}_{\alpha }〉}_{\alpha ,\beta }\Delta \theta +{O\left(\right|\nabla \mathbf{z}|}^2).\end{array}$$

(35)

Hereafter, $\Delta =div\nabla$ is the Laplace operator and ${\overline{e}}_{\alpha }=R_{\alpha }-s_{\alpha }+e_{\alpha }$, see the last formula (17).

 Proof. 

According to ([35], Section 3.2), the balance equation for the velocity holds

$$\begin{array}{c}{\partial }_t\mathbf{u}+((\mathbf{u}-{\mathbf{w}}_{\ell })·\nabla )\mathbf{u}+\frac{1}{\rho }\nabla p=\frac{1}{\rho }\{div{\mathrm{\Pi }}^{NS}+(\mathbf{u}·\nabla )\left(\rho \widehat{\mathbf{w}}\right)+(div\mathbf{u})\left(\rho \widehat{\mathbf{w}}+\ell \nabla \langle \tau {\gamma }_{\alpha }{p}_{\alpha }\rangle \right)\\ +\ell \tau \langle {\gamma }_{\alpha }{p}_{\alpha }\rangle \nabla div\mathbf{u}+\ell \nabla \left(\tau \mathbf{u}·\nabla p-\tau (\langle {\gamma }_{\alpha }{Q}_{\alpha }\rangle -Q)\right)\}+\left(1-\ell \frac{\tau }{\rho }div\left(\rho \mathbf{u}\right)\right)\mathbf{f},\end{array}$$

(36)

and the balance equation for the temperature holds

$$\begin{array}{c}{\partial }_t\theta +\left(\mathbf{u}-\frac{\langle c_{V\alpha }{\rho }_{\alpha }{\mathbf{w}}_{\ell \alpha }\rangle }{c_V\rho }\right)·\nabla \theta +\frac{R}{c_V}\theta div\mathbf{u}\\ =\frac{1}{c_V\rho }\left\{\langle c_{V\alpha }div{\mathbf{d}}_{\alpha }\rangle \theta +div(-\mathbf{q}+\langle p_{\alpha }{\mathbf{w}}_{\ell \alpha }\rangle )+\mathrm{\Pi }:\nabla \mathbf{u}-\rho \widehat{\mathbf{w}}·\mathbf{f}+Q\right\}\end{array}$$

(37)

as well. Their derivation does not exploit a specific form of ${\mathbf{w}}_{\ell \alpha }$. For ${\mathbf{d}}_1=\dots ={\mathbf{d}}_K=0$, the presented reductions have recently been proved in [37] for $\ell =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 ${\mathbf{d}}_{\alpha }$. In the balance PDEs for the mass of components (1), we obtain

$$\begin{array}{c}-div{\mathbf{d}}_{\alpha }={〈d_{\alpha \beta }\left(\Delta (G_{\alpha }-G_{\beta })+(e_{\alpha }-e_{\beta })\Delta \theta \right)〉}_{\beta }+{O\left(\right|\nabla \mathbf{z}|}^2)\\ =\theta {〈d_{\alpha \beta }\left(\frac{R_{\alpha }}{{\rho }_{\alpha }}\Delta {\rho }_{\alpha }-\frac{R_{\beta }}{{\rho }_{\beta }}\Delta {\rho }_{\beta }\right)〉}_{\beta }+{〈d_{\alpha \beta }({\overline{e}}_{\alpha }-{\overline{e}}_{\beta })〉}_{\beta }\Delta \theta +{O\left(\right|\nabla \mathbf{z}|}^2),\alpha =\overline{1,K},\end{array}$$

(38)

since the following chain of transformations

$$\begin{array}{c}\Delta G_{\alpha }=-\theta \Delta s_{\alpha }+(c_{p\alpha }-s_{\alpha })\Delta \theta +{O\left(\right|\nabla \mathbf{z}|}^2)=\theta \left(\frac{R_{\alpha }}{{\rho }_{\alpha }}\Delta {\rho }_{\alpha }-\frac{c_{V\alpha }}{\theta }\Delta \theta \right)+(c_{p\alpha }-s_{\alpha })\Delta \theta +{O\left(\right|\nabla \mathbf{z}|}^2)\\ =\theta \frac{R_{\alpha }}{{\rho }_{\alpha }}\Delta {\rho }_{\alpha }+(c_{p\alpha }-c_{V\alpha }-s_{\alpha })\Delta \theta +{O\left(\right|\nabla \mathbf{z}|}^2)\end{array}$$

is valid and $c_{p\alpha }-c_{V\alpha }=R_{\alpha }$. In the balance equation for the temperature (37), we can write

$$\begin{array}{c}\langle c_{V\alpha }div{\mathbf{d}}_{\alpha }\rangle \theta +div(-{\mathbf{q}}^d)=\langle c_{V\alpha }div{\mathbf{d}}_{\alpha }\rangle \theta -(G_{\alpha }+e_{\alpha }\theta )div{\mathbf{d}}_{\alpha }+{O\left(\right|\nabla \mathbf{z}|}^2)\\ =-〈(div{\mathbf{d}}_{\alpha })(c_{p\alpha }-c_{V\alpha }-s_{\alpha }+e_{\alpha })\theta 〉+{O\left(\right|\nabla \mathbf{z}|}^2)\\ ={\theta }^2{〈d_{\alpha \beta }\left(\frac{R_{\alpha }}{{\rho }_{\alpha }}\Delta {\rho }_{\alpha }-\frac{R_{\beta }}{{\rho }_{\beta }}\Delta {\rho }_{\beta }\right){\overline{e}}_{\alpha }〉}_{\alpha ,\beta }+\theta {〈d_{\alpha \beta }({\overline{e}}_{\alpha }-{\overline{e}}_{\beta }){\overline{e}}_{\alpha }〉}_{\alpha ,\beta }\Delta \theta +{O\left(\right|\nabla \mathbf{z}|}^2)\end{array}$$

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 $\tau$ and ℓ.

3.2. Linearized Regularized System of PDEs, Its Properties and Corollaries

In the case $\mathbf{f}=0$ and $Q_1=\dots =Q_K=0$, the system of PDEs (1)–(3) has constant solutions

$$(\rho ,\mathbf{u},\theta )(x,t)\equiv {\mathbf{z}}_0=({\rho }_0,{\mathbf{u}}_0,{\theta }_0),\mathrm{with}\mathrm{any}{\rho }_0:=({\rho }_{10},\dots ,{\rho }_{K0}),{\rho }_{10}>0,\dots ,{\rho }_{K0}>0,{\theta }_0>0$$

and any ${\mathbf{u}}_0$. We perform the linearization of the solution $\mathbf{z}$ to this system on ${\mathbf{z}}_0$ and write

$$\begin{array}{c}{\rho }_{\alpha }={\rho }_{\alpha 0}+{\rho }_{\alpha *}{\tilde{\rho }}_{\alpha }(\alpha =\overline{1,K}),\mathbf{u}={\mathbf{u}}_0+u_{*}\tilde{\mathbf{u}},\theta ={\theta }_0+{\theta }_{*}\tilde{\theta },\end{array}$$

where $\tilde{\mathbf{z}}:=(\tilde{\rho },\tilde{\mathbf{u}},\tilde{\theta })$ with $\tilde{\rho }:=({\tilde{\rho }}_1,\dots ,{\tilde{\rho }}_K)$ is the vector of dimensionless perturbations and ${\rho }_{\alpha *}>0$, $u_{*}>0$ and ${\theta }_{*}>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 $\tilde{\mathbf{z}}$ and its first and second order derivatives using the formula $\nabla \mathbf{z}=({\rho }_{1*}\nabla {\tilde{\rho }}_1,\dots ,{\rho }_{K*}\nabla {\tilde{\rho }}_K,u_{*}\nabla \tilde{\mathbf{u}},{\theta }_{*}\nabla \tilde{\theta }).$ Then, we divide the resulting PDEs by ${\rho }_{\alpha *}$, $u_{*}$ and ${\theta }_{*}$, respectively, and derive the linearized system of PDEs with constant coefficients for $\tilde{\mathbf{z}}$:

$$\begin{array}{c}{\partial }_t{\tilde{\rho }}_{\alpha }+u_{*}\left({\widehat{\mathbf{u}}}_0·\nabla {\tilde{\rho }}_{\alpha }+{\widehat{\rho }}_{\alpha 0}div\tilde{\mathbf{u}}\right)\\ ={\tau }_0{u}_{*}^2\left\{\frac{{\widehat{\rho }}_{\alpha 0}{\theta }_0}{{\rho }_0{u}_{*}^2}{\langle R_{\beta }{\rho }_{\beta *}\Delta {\tilde{\rho }}_{\beta }\rangle }_{\beta }+\ell {({\widehat{\mathbf{u}}}_0·\nabla )}^2{\tilde{\rho }}_{\alpha }+(\ell +1){\widehat{\rho }}_{\alpha 0}({\widehat{\mathbf{u}}}_0·\nabla )div\tilde{\mathbf{u}}+\frac{R_0{\widehat{\rho }}_{\alpha 0}{\theta }_{*}}{u_{*}^2}\Delta \tilde{\theta }\right\}\\ +{\theta }_0{〈d_{\alpha \beta 0}\left(\frac{R_{\alpha }}{{\rho }_{\alpha 0}}\Delta {\tilde{\rho }}_{\alpha }-\frac{R_{\beta }{\rho }_{\beta *}}{{\rho }_{\beta 0}{\rho }_{\alpha *}}\Delta {\tilde{\rho }}_{\beta }\right)〉}_{\beta }+\frac{{\theta }_{*}}{{\rho }_{\alpha *}}{〈d_{\alpha \beta 0}({\overline{e}}_{\alpha 0}-{\overline{e}}_{\beta 0})〉}_{\beta }\Delta \tilde{\theta },\alpha =\overline{1,K},\\ {\partial }_t\tilde{\mathbf{u}}+u_{*}\left(\frac{{\theta }_0}{{\rho }_0{u}_{*}^2}{\langle R_{\alpha }{\rho }_{\alpha *}\nabla {\tilde{\rho }}_{\alpha }\rangle }_{\beta }+({\widehat{\mathbf{u}}}_0·\nabla )\tilde{\mathbf{u}}+\frac{R_0{\theta }_{*}}{u_{*}^2}\nabla \tilde{\theta }\right)=u_{*}^2\{(\ell +1){\tau }_0\frac{{\theta }_0}{{\rho }_0{u}_{*}^2}({\widehat{\mathbf{u}}}_0·\nabla )\langle R_{\alpha }{\rho }_{\alpha *}\nabla {\tilde{\rho }}_{\alpha }\rangle \\ +\frac{{\mu }_0}{{\rho }_0{u}_{*}^2}\Delta \tilde{\mathbf{u}}+\left(\frac{{\chi }_0}{{\rho }_0{u}_{*}^2}+\ell {\tau }_0\frac{{\left(R\gamma \right)}_0{\theta }_0}{u_{*}^2}\right)\nabla div\tilde{\mathbf{u}}+{\tau }_0{({\widehat{\mathbf{u}}}_0·\nabla )}^2\tilde{\mathbf{u}}+(\ell +1){\tau }_0\frac{R_0{\theta }_{*}}{u_{*}^2}({\widehat{\mathbf{u}}}_0·\nabla )\nabla \tilde{\theta }\},\\ {\partial }_t\tilde{\theta }+u_{*}\left(\frac{R_0{\widehat{\theta }}_0}{c_{V0}}div\tilde{\mathbf{u}}+{\widehat{\mathbf{u}}}_0·\nabla \tilde{\theta }\right)=u_{*}^2\{{\tau }_0\frac{R_0{\widehat{\theta }}_0^2{\theta }_{*}}{c_{V0}{\rho }_0{u}_{*}^2}\langle R_{\alpha }{\rho }_{\alpha *}\Delta {\tilde{\rho }}_{\alpha }\rangle \\ +(\ell +1){\tau }_0\frac{R_0{\widehat{\theta }}_0}{c_{V0}}({\widehat{\mathbf{u}}}_0·\nabla )div\tilde{\mathbf{u}}+\ell {\tau }_0{({\widehat{\mathbf{u}}}_0·\nabla )}^2\tilde{\theta }+\left(\frac{{\varkappa }_0}{c_{V0}{\rho }_0{u}_{*}^2}+{\tau }_0\frac{R_0^2{\theta }_0}{c_{V0}{u}_{*}^2}\right)\Delta \tilde{\theta }\}\\ +\frac{{\widehat{\theta }}_0{\theta }_0{\rho }_{\alpha *}}{c_{V0}{\rho }_0}{〈d_{\alpha \beta 0}\left(\frac{R_{\alpha }}{{\rho }_{\alpha 0}}\Delta {\tilde{\rho }}_{\alpha }-\frac{R_{\beta }{\rho }_{\beta *}}{{\rho }_{\beta 0}{\rho }_{\alpha *}}\Delta {\tilde{\rho }}_{\beta }\right){\overline{e}}_{\alpha 0}〉}_{\alpha ,\beta }+\frac{{\theta }_0}{c_{V0}{\rho }_0}{〈d_{\alpha \beta 0}({\overline{e}}_{\alpha 0}-{\overline{e}}_{\beta 0}){\overline{e}}_{\alpha 0}〉}_{\alpha ,\beta }\Delta \tilde{\theta }\}.\end{array}$$

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 $\rho$, R and $c_V$ and the average value of $R_{\alpha }{\gamma }_{\alpha }$:

$$\begin{array}{c}{\widehat{\rho }}_{\alpha 0}:=\frac{{\rho }_{\alpha 0}}{{\rho }_{\alpha *}},{\widehat{\mathbf{u}}}_0=({\widehat{u}}_{10},\dots ,{\widehat{u}}_{n0}):=\frac{{\mathbf{u}}_0}{u_{*}},{\widehat{\theta }}_0:=\frac{{\theta }_0}{{\theta }_{*}},\\ {\rho }_0:=\langle {\rho }_{\alpha 0}\rangle ,R_0:=〈\frac{{\rho }_{\alpha 0}}{{\rho }_0}{R}_{\alpha }〉,c_{V0}:=〈\frac{{\rho }_{\alpha 0}}{{\rho }_0}{c}_{V\alpha }〉,{\left(R\gamma \right)}_0=〈\frac{{\rho }_{\alpha 0}}{{\rho }_0}{R}_{\alpha }{\gamma }_{\alpha }〉.\end{array}$$

In addition, $d_{\alpha \beta 0}$, ${\overline{e}}_{\alpha 0}$, ${\tau }_0$, ${\mu }_0$, ${\chi }_0$ and ${\varkappa }_0$ are the values of $d_{\alpha \beta }$, ${\overline{e}}_{\alpha }$, $\tau$, $\mu$, $\chi$ and $\varkappa$, respectively, on the background solution ${\mathbf{z}}_0$, and the following PDE operator is involved

$${({\widehat{\mathbf{u}}}_0·\nabla )}^2:=({\widehat{\mathbf{u}}}_0·\nabla )({\widehat{\mathbf{u}}}_0·\nabla )=\sum _{i,j=1}^n{\widehat{u}}_{0i}{\widehat{u}}_{0j}{\partial }_i{\partial }_j.$$

For ${\mathbf{d}}_1=\dots ={\mathbf{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

$$\begin{array}{c}{\rho }_{\alpha *}=b\sqrt{\frac{{\rho }_{\alpha 0}{c}_{V0}{\rho }_0}{R_{\alpha }}},\alpha =\overline{1,K},u_{*}=b\sqrt{c_{V0}{\theta }_0},{\theta }_{*}=b{\theta }_0\forall b>0,\end{array}$$

(39)

with a free parameter b. We accept this choice and pass to a much simpler form of the above linearized system of PDEs

$$\begin{array}{c}{\partial }_t{\tilde{\rho }}_{\alpha }+u_{*}\left({\widehat{\mathbf{u}}}_0·\nabla {\tilde{\rho }}_{\alpha }+{\widehat{\rho }}_{\alpha 0}div\tilde{\mathbf{u}}\right)\\ ={\tau }_0{u}_{*}^2\left\{{\widehat{\rho }}_{\alpha 0}{\langle {\widehat{\rho }}_{\beta 0}\Delta {\tilde{\rho }}_{\beta }\rangle }_{\beta }+\ell {({\widehat{\mathbf{u}}}_0·\nabla )}^2{\tilde{\rho }}_{\alpha }+(\ell +1){\widehat{\rho }}_{\alpha 0}({\widehat{\mathbf{u}}}_0·\nabla )div\tilde{\mathbf{u}}+a_0{\widehat{\rho }}_{\alpha 0}\Delta \tilde{\theta }\right\}\\ +b_{\alpha }{〈d_{\alpha \beta 0}(b_{\alpha }\Delta {\tilde{\rho }}_{\alpha }-b_{\beta }\Delta {\tilde{\rho }}_{\beta })〉}_{\beta }+b_{\alpha }{b}_{\left(\theta \right)}{〈d_{\alpha \beta 0}({\overline{e}}_{\alpha 0}-{\overline{e}}_{\beta 0})〉}_{\beta }\Delta \tilde{\theta },\alpha =\overline{1,K},\end{array}$$

(40)

$$\begin{array}{c}{\partial }_t\tilde{\mathbf{u}}+u_{*}\left(\langle {\widehat{\rho }}_{\alpha 0}\nabla {\tilde{\rho }}_{\alpha }\rangle +({\widehat{\mathbf{u}}}_0·\nabla )\tilde{\mathbf{u}}+a_0\nabla \tilde{\theta }\right)=u_{*}^2\{(\ell +1){\tau }_0({\widehat{\mathbf{u}}}_0·\nabla )\langle {\widehat{\rho }}_{\alpha 0}\nabla {\tilde{\rho }}_{\alpha }\rangle \\ +{\overline{\mu }}_0\Delta \tilde{\mathbf{u}}+\left({\overline{\chi }}_0+\ell {\tau }_0{\widehat{\theta }}_0{\left(a\gamma \right)}_0\right)\nabla div\tilde{\mathbf{u}}+{\tau }_0{({\widehat{\mathbf{u}}}_0·\nabla )}^2\tilde{\mathbf{u}}+(\ell +1){\tau }_0{a}_0({\widehat{\mathbf{u}}}_0·\nabla )\nabla \tilde{\theta }\},\end{array}$$

(41)

$$\begin{array}{c}{\partial }_t\tilde{\theta }+u_{*}\left(a_{0}div\tilde{\mathbf{u}}+{\widehat{\mathbf{u}}}_0·\nabla \tilde{\theta }\right)\\ =u_{*}^2\left\{{\tau }_0{a}_0\langle {\widehat{\rho }}_{\alpha 0}\Delta {\tilde{\rho }}_{\alpha }\rangle +(\ell +1){\tau }_0{a}_0({\widehat{\mathbf{u}}}_0·\nabla )div\tilde{\mathbf{u}}+\ell {\tau }_0{({\widehat{\mathbf{u}}}_0·\nabla )}^2\tilde{\theta }+\left({\overline{\varkappa }}_0+{\tau }_0{a}_0^2\right)\Delta \tilde{\theta }\right\},\\ +b_{\left(\theta \right)}{〈d_{\alpha \beta 0}(b_{\alpha }\Delta {\tilde{\rho }}_{\alpha }-b_{\beta }\Delta {\tilde{\rho }}_{\beta }){\overline{e}}_{\alpha 0}〉}_{\alpha ,\beta }+b_{\left(\theta \right)}^2{〈d_{\alpha \beta 0}({\overline{e}}_{\alpha 0}-{\overline{e}}_{\beta 0}){\overline{e}}_{\alpha 0}〉}_{\alpha ,\beta }\Delta \tilde{\theta }\}.\end{array}$$

(42)

Here, the following constant factors have been introduced

$$\begin{array}{c}{a}_{\alpha }:=\frac{R_{\alpha }{\theta }_{*}}{u_{*}^2},a_0:=〈\frac{{\rho }_{\alpha 0}}{{\rho }_0}{a}_{\alpha }〉=\frac{R_0{\theta }_{*}}{u_{*}^2},{\left(a\gamma \right)}_0=〈\frac{{\rho }_{\alpha 0}}{{\rho }_0}{a}_{\alpha }{\gamma }_{\alpha }〉,\\ b_{\alpha }:={\left(\frac{R_{\alpha }{\theta }_0}{{\rho }_{\alpha 0}}\right)}^{1/2},b_{\left(\theta \right)}:={\left(\frac{{\theta }_0}{c_{V0}{\rho }_0}\right)}^{1/2},{\overline{\mu }}_0:=\frac{{\mu }_0}{{\rho }_0{u}_{*}^2},{\overline{\chi }}_0:=\frac{{\chi }_0}{{\rho }_0{u}_{*}^2},{\overline{\varkappa }}_0:=\frac{{\varkappa }_0}{c_{V0}{\rho }_0{u}_{*}^2},\end{array}$$

and, for the last two terms in (40) and (42), we have taken into account the formulas

$${\theta }_0\frac{R_{\beta }{\rho }_{\beta *}}{{\rho }_{\beta 0}{\rho }_{\alpha *}}=b_{\alpha }{b}_{\beta },\frac{{\theta }_{*}}{{\rho }_{\alpha *}}=b_{\alpha }{b}_{\left(\theta \right)},\frac{{\widehat{\theta }}_0{\rho }_{\alpha *}}{c_{V0}{\rho }_0}=\frac{b_{\left(\theta \right)}}{b_{\alpha }},\alpha ,\beta =\overline{1,K}.$$

Next, we study the initial-boundary value problem (IBVP) for the linearized system of PDEs (40)–(42) in the cylinder $Q_{\infty }:=\mathsf{\Omega }\times (0,\infty )$ under the boundary and initial conditions

$$\begin{array}{c}\tilde{\mathbf{z}}{{|}_{\partial \mathsf{\Omega }\times (0,\infty )}=0,\tilde{\mathbf{z}}|}_{t=0}={\tilde{\mathbf{z}}}^{\left(0\right)}\left(x\right).\end{array}$$

(43)

Let $L^2\left(\mathsf{\Omega }\right)$ and ${\mathbf{L}}^2\left(\mathsf{\Omega }\right)$ be, respectively, the standard Lebesgue spaces of functions and vector functions defined on $\mathsf{\Omega }$ and denote by ${(·,·)}_{\mathsf{\Omega }}={(·,·)}_{L^2\left(\mathsf{\Omega }\right)}$, ${\parallel ·\parallel }_{\mathsf{\Omega }}={\parallel ·\parallel }_{L^2\left(\mathsf{\Omega }\right)}$, ${(·,·)}_{\mathsf{\Omega }}={(·,·)}_{{\mathbf{L}}^2\left(\mathsf{\Omega }\right)}$ and ${\parallel ·\parallel }_{\mathsf{\Omega }}={\parallel ·\parallel }_{{\mathbf{L}}^2\left(\mathsf{\Omega }\right)}$ their inner products and norms. Let ${\mathbf{H}}^1\left(\mathsf{\Omega }\right)={\mathbf{W}}_2^1\left(\mathsf{\Omega }\right)$ be a standard Sobolev space of vector functions defined on $\mathsf{\Omega }$, and ${\mathbf{H}}_0^1\left(\mathsf{\Omega }\right)$ be the closure in the ${\mathbf{H}}^1\left(\mathsf{\Omega }\right)$-norm of the space of smooth vector-functions with a compact support in $\mathsf{\Omega }$.

For ${\partial }_t\tilde{\mathbf{z}}(·,t),\nabla \tilde{\mathbf{z}}(·,t)\in {\mathbf{L}}^2\left(\mathsf{\Omega }\right)$, PDEs (40)–(42) correspond to the integral identity

$$\begin{array}{c}{\left({\partial }_t\tilde{\mathbf{z}}(·,t),\mathbf{z}\right)}_{\mathsf{\Omega }}+u_{*}{\mathcal{B}}_{\mathsf{\Omega }}(\tilde{\mathbf{z}}(·,t),\mathbf{z})+u_{*}^2{\mathcal{A}}_{\mathsf{\Omega }}(\tilde{\mathbf{z}}(·,t),\mathbf{z})+{\mathcal{A}}_{\mathsf{\Omega }}^d(\tilde{\mathbf{z}}(·,t),\mathbf{z})=0\forall \mathbf{z}\in {\mathbf{H}}_0^1\left(\mathsf{\Omega }\right),\end{array}$$

(44)

for any vector-function $\mathbf{z}=(\rho ,\mathbf{u},\theta )\left(x\right)\in {\mathbf{H}}_0^1\left(\mathsf{\Omega }\right)$ (which here is not the solution to the original quasilinear system of PDEs) and almost all $t>0$.

In the identity, the three bilinear forms are involved

$$\begin{array}{c}{\mathcal{B}}_{\mathsf{\Omega }}(\tilde{\mathbf{z}},\mathbf{z}):=〈{\left({\widehat{\mathbf{u}}}_0·\nabla {\tilde{\rho }}_{\alpha }+{\widehat{\rho }}_{\alpha 0}div\tilde{\mathbf{u}},{\rho }_{\alpha }\right)}_{\mathsf{\Omega }}〉\\ +{\left(\langle {\widehat{\rho }}_{\alpha 0}\nabla {\tilde{\rho }}_{\alpha }\rangle +({\widehat{\mathbf{u}}}_0·\nabla )\tilde{\mathbf{u}}+a_0\nabla \tilde{\theta },\mathbf{u}\right)}_{\mathsf{\Omega }}+{\left(a_{0}div\tilde{\mathbf{u}}+{\widehat{\mathbf{u}}}_0·\nabla \tilde{\theta },\theta \right)}_{\mathsf{\Omega }},\\ {\mathcal{A}}_{\mathsf{\Omega }}(\tilde{\mathbf{z}},\mathbf{z}):={\overline{\mu }}_0{\left(\nabla \tilde{\mathbf{u}},\nabla \mathbf{u}\right)}_{\mathsf{\Omega }}+{\overline{\chi }}_0{\left(div\tilde{\mathbf{u}},div\mathbf{u}\right)}_{\mathsf{\Omega }}+{\overline{\varkappa }}_0{\left(\nabla \tilde{\theta },\nabla \theta \right)}_{\mathsf{\Omega }}\\ +{\tau }_0\{{\left(\langle {\widehat{\rho }}_{\alpha 0}\nabla {\tilde{\rho }}_{\alpha }\rangle ,\langle {\widehat{\rho }}_{\alpha 0}\nabla {\rho }_{\alpha }\rangle \right)}_{\mathsf{\Omega }}+\ell 〈{\left({\widehat{\mathbf{u}}}_0·\nabla {\tilde{\rho }}_{\alpha },({\widehat{\mathbf{u}}}_0·\nabla ){\rho }_{\alpha }\right)}_{\mathsf{\Omega }}〉+(\ell +1){\left(({\widehat{\mathbf{u}}}_0·\nabla )\tilde{\mathbf{u}},\langle {\widehat{\rho }}_{\alpha 0}\nabla {\rho }_{\alpha }\rangle \right)}_{\mathsf{\Omega }}\\ +{\left(a_0\nabla \tilde{\theta },\langle {\widehat{\rho }}_{\alpha 0}\nabla {\rho }_{\alpha }\rangle \right)}_{\mathsf{\Omega }}+(\ell +1){\left(\langle {\widehat{\rho }}_{\alpha 0}\nabla {\tilde{\rho }}_{\alpha }\rangle ,({\widehat{\mathbf{u}}}_0·\nabla )\mathbf{u}\right)}_{\mathsf{\Omega }}+\ell {\left({\widehat{\theta }}_0{\left(a\gamma \right)}_{0}div\tilde{\mathbf{u}},div\mathbf{u}\right)}_{\mathsf{\Omega }}\\ +{\left(({\widehat{\mathbf{u}}}_0·\nabla )\tilde{\mathbf{u}},({\widehat{\mathbf{u}}}_0·\nabla )\mathbf{u}\right)}_{\mathsf{\Omega }}+(\ell +1){\left(a_0\nabla \tilde{\theta },({\widehat{\mathbf{u}}}_0·\nabla )\mathbf{u}\right)}_{\mathsf{\Omega }}\\ +{\left(〈{\widehat{\rho }}_{\alpha 0}\nabla {\tilde{\rho }}_{\alpha }〉,a_0\nabla \theta \right)}_{\mathsf{\Omega }}+(\ell +1){\left(({\widehat{\mathbf{u}}}_0·\nabla )\tilde{\mathbf{u}},a_0\nabla \theta \right)}_{\mathsf{\Omega }}+\ell {\left({\widehat{\mathbf{u}}}_0·\nabla \tilde{\theta },{\widehat{\mathbf{u}}}_0·\nabla \theta \right)}_{\mathsf{\Omega }}+{\left(a_0\nabla \tilde{\theta },a_0\nabla \theta \right)}_{\mathsf{\Omega }}\},\end{array}$$

where the tensors $\nabla \tilde{\mathbf{u}}$ and $\nabla \mathbf{u}$ are considered as vectors of length $n^2$, and

$$\begin{array}{c}{\mathcal{A}}_{\mathsf{\Omega }}^d(\tilde{\mathbf{z}},\mathbf{z}):={〈{\left(d_{\alpha \beta 0}(b_{\alpha }\nabla {\tilde{\rho }}_{\alpha }-b_{\beta }\nabla {\tilde{\rho }}_{\beta }),b_{\alpha }\nabla {\rho }_{\alpha }\right)}_{\mathsf{\Omega }}〉}_{\alpha ,\beta }+{〈{\left(d_{\alpha \beta 0}{b}_{\left(\theta \right)}({\overline{e}}_{\alpha 0}-{\overline{e}}_{\beta 0})\nabla \tilde{\theta },b_{\alpha }\nabla {\rho }_{\alpha }\right)}_{\mathsf{\Omega }}〉}_{\alpha ,\beta }\\ +{〈{\left(d_{\alpha \beta 0}(b_{\alpha }\nabla {\tilde{\rho }}_{\alpha }-b_{\beta }\nabla {\tilde{\rho }}_{\beta }),b_{\left(\theta \right)}{\overline{e}}_{\alpha 0}\nabla \theta \right)}_{\mathsf{\Omega }}〉}_{\alpha ,\beta }+{〈{\left(d_{\alpha \beta 0}{b}_{\left(\theta \right)}({\overline{e}}_{\alpha 0}-{\overline{e}}_{\beta 0})\nabla \tilde{\theta },b_{\left(\theta \right)}{\overline{e}}_{\alpha 0}\nabla \theta \right)}_{\mathsf{\Omega }}〉}_{\alpha ,\beta }.\end{array}$$

The last bilinear form corresponds to the diffusive terms and is independent of $\tilde{\mathbf{u}}$ and $\mathbf{u}$. Formally, identity (44) arises after multiplying Equations (40)–(42) by ${\rho }_{\alpha }$, $\mathbf{u}$ and $\theta$, respectively, integrating over $\mathsf{\Omega }$ and by parts and summing up the results.

Let us study properties of the defined bilinear forms that is crucial below.

 Lemma 2.

Let $d_{\alpha \beta 0}=d_{\beta \alpha 0}$ for any $\alpha \ne \beta$. The following skew symmetry and symmetry properties hold

$$\begin{array}{c}{\mathcal{B}}_{\mathsf{\Omega }}(\tilde{\mathbf{z}},\mathbf{z})=-{\mathcal{B}}_{\mathsf{\Omega }}(\mathbf{z},\tilde{\mathbf{z}})\forall \tilde{\mathbf{z}}\in {\mathbf{H}}^1\left(\mathsf{\Omega }\right),\mathbf{z}\in {\mathbf{H}}_0^1\left(\mathsf{\Omega }\right),\end{array}$$

(45)

$$\begin{array}{c}{\mathcal{A}}_{\mathsf{\Omega }}(\tilde{\mathbf{z}},\mathbf{z})={\mathcal{A}}_{\mathsf{\Omega }}(\mathbf{z},\tilde{\mathbf{z}}),{\mathcal{A}}_{\mathsf{\Omega }}^d(\tilde{\mathbf{z}},\mathbf{z})={\mathcal{A}}_{\mathsf{\Omega }}^d(\mathbf{z},\tilde{\mathbf{z}})\forall \tilde{\mathbf{z}},\mathbf{z}\in {\mathbf{H}}^1\left(\mathsf{\Omega }\right).\end{array}$$

(46)

Moreover, let $\mathbf{z}=(\rho ,\mathbf{u},\theta )\in {\mathbf{H}}^1\left(\mathsf{\Omega }\right)$ with $\mathbf{u}\in {\mathbf{H}}_0^1\left(\mathsf{\Omega }\right)$. The following representations for the quadratic forms hold

$$\begin{array}{c}{\mathcal{A}}_{\mathsf{\Omega }}(\mathbf{z},\mathbf{z})={\overline{\mu }}_0{\parallel \nabla \mathbf{u}\parallel }_{\mathsf{\Omega }}^2+{\overline{\chi }}_0{\parallel div\mathbf{u}\parallel }_{\mathsf{\Omega }}^2+{\overline{\varkappa }}_0{\parallel \nabla \theta \parallel }_{\mathsf{\Omega }}^2+{\tau }_0{\parallel \langle {\widehat{\rho }}_{\alpha 0}\nabla {\rho }_{\alpha }\rangle +({\widehat{\mathbf{u}}}_0·\nabla )\mathbf{u}+a_0\nabla \theta \parallel }_{\mathsf{\Omega }}^2\\ +\ell {\tau }_0\left(〈\parallel {\widehat{\mathbf{u}}}_0·\nabla {\rho }_{\alpha }+{\widehat{\rho }}_{\alpha 0}{div\mathbf{u}\parallel }_{\mathsf{\Omega }}^2〉+\parallel a_{0}div\mathbf{u}+{\widehat{\mathbf{u}}}_0·\nabla \theta {\parallel }_{\mathsf{\Omega }}^2+g_0{\parallel div\mathbf{u}\parallel }_{\mathsf{\Omega }}^2\right),\end{array}$$

(47)

$$\begin{array}{c}{\mathcal{A}}_{\mathsf{\Omega }}^d(\mathbf{z},\mathbf{z})=\frac{1}{2}{〈d_{\alpha \beta 0}{\parallel b_{\alpha }\nabla {\rho }_{\alpha }-b_{\beta }\nabla {\rho }_{\beta }+b_{\left(\theta \right)}({\overline{e}}_{\alpha 0}-{\overline{e}}_{\beta 0})\nabla \theta \parallel }_{\mathsf{\Omega }}^2〉}_{\alpha ,\beta },\end{array}$$

(48)

with $g_0:=\frac{1}{u_{*}^2}(\langle {\gamma }_{\alpha }{p}_{\alpha 0}\rangle -{\gamma }_0{p}_0)⩾0$, see relations (22), where $p_{\alpha 0}$, ${\gamma }_0$ and $p_0$ are the values of $p_{\alpha }$, γ and p on the background solution ${\mathbf{z}}_0$.

 Proof. 

To prove the skew symmetry property (45), we integrate by parts term by term in the definition of ${\mathcal{B}}_{\mathsf{\Omega }}(\tilde{\mathbf{z}},\mathbf{z})$, rearrange the summands and obtain the equalities

$$\begin{array}{c}{\mathcal{B}}_{\mathsf{\Omega }}(\tilde{\mathbf{z}},\mathbf{z})=-\langle {\left({\tilde{\rho }}_{\alpha },{\widehat{\mathbf{u}}}_0·\nabla {\rho }_{\alpha }\right)}_{\mathsf{\Omega }}\rangle -\left(\tilde{\mathbf{u}},\langle {\widehat{\rho }}_{\alpha 0}\nabla {\rho }_{\alpha }\right.〉{)}_{\mathsf{\Omega }}-〈{\left({\tilde{\rho }}_{\alpha },{\widehat{\rho }}_{\alpha 0}div\mathbf{u}\right)}_{\mathsf{\Omega }}〉-{\left(\tilde{\mathbf{u}},({\widehat{\mathbf{u}}}_0·\nabla )\mathbf{u}\right)}_{\mathsf{\Omega }}\\ -{\left(\tilde{\theta },a_{0}div\mathbf{u}\right)}_{\mathsf{\Omega }}-(\tilde{\mathbf{u}},a_0{\nabla \theta )}_{\mathsf{\Omega }}-{\left(\tilde{\theta },{\widehat{\mathbf{u}}}_0·\nabla \theta \right)}_{\mathsf{\Omega }}=-{\mathcal{B}}_{\mathsf{\Omega }}(\mathbf{z},\tilde{\mathbf{z}}).\end{array}$$

The first symmetry property (46) is obvious. Due to identity (16) applied to each term, the following formulas hold

$$\begin{array}{c}2{\mathcal{A}}_{\mathsf{\Omega }}^d(\tilde{\mathbf{z}},\mathbf{z})={〈{\left(d_{\alpha \beta 0}(b_{\alpha }\nabla {\tilde{\rho }}_{\alpha }-b_{\beta }\nabla {\tilde{\rho }}_{\beta }),b_{\alpha }\nabla {\rho }_{\alpha }-b_{\beta }\nabla {\rho }_{\beta }\right)}_{\mathsf{\Omega }}〉}_{\alpha ,\beta }\\ +{〈{\left(d_{\alpha \beta 0}{b}_{\left(\theta \right)}({\overline{e}}_{\alpha 0}-{\overline{e}}_{\beta 0})\nabla \tilde{\theta },b_{\alpha }\nabla {\rho }_{\alpha }-b_{\beta }\nabla {\rho }_{\beta }\right)}_{\mathsf{\Omega }}〉}_{\alpha ,\beta }\\ +{〈{\left(d_{\alpha \beta 0}(b_{\alpha }\nabla {\tilde{\rho }}_{\alpha }-b_{\beta }\nabla {\tilde{\rho }}_{\beta }),b_{\left(\theta \right)}({\overline{e}}_{\alpha 0}-{\overline{e}}_{\beta 0})\nabla \theta \right)}_{\mathsf{\Omega }}〉}_{\alpha ,\beta }\\ +{〈{\left(d_{\alpha \beta 0}{b}_{\left(\theta \right)}({\overline{e}}_{\alpha 0}-{\overline{e}}_{\beta 0})\nabla \tilde{\theta },b_{\left(\theta \right)}({\overline{e}}_{\alpha 0}-{\overline{e}}_{\beta 0})\nabla \theta \right)}_{\mathsf{\Omega }}〉}_{\alpha ,\beta }\\ ={〈{\left(d_{\alpha \beta 0}\left(b_{\alpha }\nabla {\tilde{\rho }}_{\alpha }-b_{\beta }\nabla {\tilde{\rho }}_{\beta }+b_{\left(\theta \right)}({\overline{e}}_{\alpha 0}-{\overline{e}}_{\beta 0})\nabla \tilde{\theta }\right),b_{\alpha }\nabla {\rho }_{\alpha }-b_{\beta }\nabla {\rho }_{\beta }+b_{\left(\theta \right)}({\overline{e}}_{\alpha 0}-{\overline{e}}_{\beta 0})\nabla \theta \right)}_{\mathsf{\Omega }}〉}_{\alpha ,\beta }.\end{array}$$

They imply the symmetry property (46) for ${\mathcal{A}}_{\mathsf{\Omega }}^d(\tilde{\mathbf{z}},\mathbf{z})$ and representation (48) for ${\mathcal{A}}_{\mathsf{\Omega }}^d(\mathbf{z},\mathbf{z})$.

Property (47) for $\ell =0$ has recently been checked in [37] (for $\mathsf{\Omega }={\mathbb{R}}^n$ that is not essential). The rest of the terms in ${\mathcal{A}}_{\mathsf{\Omega }}(\mathbf{z},\mathbf{z})$ are as follows

$$\begin{array}{c}\ell {\tau }_0\{〈\parallel {\widehat{\mathbf{u}}}_0·\nabla {\rho }_{\alpha }{\parallel }_{\mathsf{\Omega }}^2〉+2{\left(\langle {\widehat{\rho }}_{\alpha 0}\nabla {\rho }_{\alpha }\rangle ,({\widehat{\mathbf{u}}}_0·\nabla )\mathbf{u}\right)}_{\mathsf{\Omega }}+2{\left(({\widehat{\mathbf{u}}}_0·\nabla )\mathbf{u},a_0\nabla \theta \right)}_{\mathsf{\Omega }}\\ +{\widehat{\theta }}_0{\left(a\gamma \right)}_0{\parallel div\mathbf{u}\parallel }_{\mathsf{\Omega }}^2+{\parallel {\widehat{\mathbf{u}}}_0·\nabla \theta \parallel }_{\mathsf{\Omega }}^2\}=:{\mathcal{A}}_{\ell \mathsf{\Omega }}(\mathbf{z},\mathbf{z}).\end{array}$$

Next, we recall the following algebraic formula and integral identity

$$\begin{array}{c}〈{\left({\widehat{\mathbf{u}}}_0·\nabla {\rho }_{\alpha }+{\widehat{\rho }}_{\alpha 0}div\mathbf{u}\right)}^2〉=〈{\left({\widehat{\mathbf{u}}}_0·\nabla {\rho }_{\alpha }\right)}^2〉\\ +\left(〈{\widehat{\rho }}_{\alpha 0}^2〉+a_0^2\right){(div\mathbf{u})}^2+{\left({\widehat{\mathbf{u}}}_0·\nabla \theta \right)}^2+2\left({\widehat{\mathbf{u}}}_0·\langle {\widehat{\rho }}_{\alpha 0}\nabla {\rho }_{\alpha }\rangle \right)div\mathbf{u}+2a_0({\widehat{\mathbf{u}}}_0·\nabla \theta )div\mathbf{u},\\ {({\widehat{\mathbf{u}}}_0·\nabla \phi ,div\mathbf{u})}_{\mathsf{\Omega }}={(\nabla \phi ,({\widehat{\mathbf{u}}}_0·\nabla )\mathbf{u})}_{\mathsf{\Omega }}\forall \phi \in H^1\left(\mathsf{\Omega }\right),\mathbf{u}\in {\mathbf{H}}_0^1\left(\mathsf{\Omega }\right),\end{array}$$

(49)

see ([35], formulas (4.14) and (4.18)). Integrating Formula (49) over $\mathsf{\Omega }$ and applying the last identity, we find

$${\mathcal{A}}_{\ell \mathsf{\Omega }}(\mathbf{z},\mathbf{z})=\ell {\tau }_0\left(〈\parallel {\widehat{\mathbf{u}}}_0·\nabla {\rho }_{\alpha }+{\widehat{\rho }}_{\alpha 0}{div\mathbf{u}\parallel }_{\mathsf{\Omega }}^2〉+\parallel a_{0}div\mathbf{u}+{\widehat{\mathbf{u}}}_0·\nabla \theta {\parallel }_{\mathsf{\Omega }}^2+g_0{\parallel div\mathbf{u}\parallel }_{\mathsf{\Omega }}^2\right)$$

with $g_0:={\widehat{\theta }}_0{\left(a\gamma \right)}_0-(〈{\widehat{\rho }}_{\alpha 0}^2〉+a_0^2)$. According to ([35], proof of Lemma 3.1) and relations (22), we obtain $g_0=\frac{1}{u_{*}^2}(\langle {\gamma }_{\alpha }{p}_{\alpha 0}\rangle -{\gamma }_0{p}_0)⩾0$ that completes the proof of representation (47). □

 Corollary 1.

Let $d_{\alpha \beta 0}=d_{\beta \alpha 0}⩾d>0$ for any $\alpha \ne \beta$. The following positive definiteness inequality holds

$$\begin{array}{c}{u}_{*}^2{\mathcal{A}}_{\mathsf{\Omega }}(\mathbf{z},\mathbf{z})+{\mathcal{A}}_{\mathsf{\Omega }}^d(\mathbf{z},\mathbf{z})⩾max\left\{{\delta }_1〈\parallel \nabla {\rho }_{\alpha }{\parallel }_{\mathsf{\Omega }}^2〉,{\delta }_0\left({\parallel \nabla \mathbf{u}\parallel }_{\mathsf{\Omega }}^2+{\parallel \nabla \theta \parallel }_{\mathsf{\Omega }}^2\right)\right\},\end{array}$$

(50)

for any $\mathbf{z}=(\rho ,\mathbf{u},\theta )\in {\mathbf{H}}^1\left(\mathsf{\Omega }\right)$ with $\mathbf{u}\in {\mathbf{H}}_0^1\left(\mathsf{\Omega }\right)$, with ${\delta }_0:=u_{*}^{2}min\{{\overline{\mu }}_0,{\overline{\varkappa }}_0\}>0$ and some ${\delta }_1>0$.

 Proof. 

Clearly, representations (47) and (48) imply the lower bound

$$\begin{array}{c}{u}_{*}^2{\mathcal{A}}_{\mathsf{\Omega }}(\mathbf{z},\mathbf{z})+{\mathcal{A}}_{\mathsf{\Omega }}^d(\mathbf{z},\mathbf{z})⩾u_{*}^2\left({\overline{\mu }}_0{\parallel \nabla \mathbf{u}\parallel }_{\mathsf{\Omega }}^2+{\overline{\chi }}_0{\parallel div\mathbf{u}\parallel }_{\mathsf{\Omega }}^2+{\overline{\varkappa }}_0{\parallel \nabla \theta \parallel }_{\mathsf{\Omega }}^2\right)\\ +u_{*}^2{\tau }_0{\parallel \langle {\widehat{\rho }}_{\alpha 0}\nabla {\rho }_{\alpha }\rangle +({\widehat{\mathbf{u}}}_0·\nabla )\mathbf{u}+a_0\nabla \theta \parallel }_{\mathsf{\Omega }}^2+\frac{1}{2}d{〈\parallel b_{\alpha }\nabla {\rho }_{\alpha }-b_{\beta }\nabla {\rho }_{\beta }+b_{\left(\theta \right)}({\overline{e}}_{\alpha 0}-{\overline{e}}_{\beta 0}){\nabla \theta \parallel }_{\mathsf{\Omega }}^2〉}_{\alpha ,\beta }.\end{array}$$

We further apply simple bounds for the terms containing $\nabla {\rho }_{\alpha }$:

$$\begin{array}{c}\parallel \langle {\widehat{\rho }}_{\alpha 0}\nabla {\rho }_{\alpha }\rangle {\parallel }_{\mathsf{\Omega }}^2⩽2\left(\parallel \langle {\widehat{\rho }}_{\alpha 0}\nabla {\rho }_{\alpha }\rangle +({\widehat{\mathbf{u}}}_0·\nabla )\mathbf{u}+a_0{\nabla \theta \parallel }_{\mathsf{\Omega }}^2+2|{\widehat{\mathbf{u}}}_0{|}^2{\parallel \nabla \mathbf{u}\parallel }_{\mathsf{\Omega }}^2+a_0^2{\parallel \nabla \theta \parallel }_{\mathsf{\Omega }}^2\right),\\ \parallel b_{\alpha }\nabla {\rho }_{\alpha }-b_{\beta }\nabla {\rho }_{\beta }{\parallel }_{\mathsf{\Omega }}^2⩽2\parallel b_{\alpha }\nabla {\rho }_{\alpha }-b_{\beta }\nabla {\rho }_{\beta }+b_{\left(\theta \right)}({\overline{e}}_{\alpha 0}-{\overline{e}}_{\beta 0}){\nabla \theta \parallel }_{\mathsf{\Omega }}^2+2b_{\left(\theta \right)}^2{({\overline{e}}_{\alpha 0}-{\overline{e}}_{\beta 0})}^2{\parallel \nabla \theta \parallel }_{\mathsf{\Omega }}^2\end{array}$$

and

$$\begin{array}{c}{b}_{\alpha }^2{\parallel \nabla {\rho }_{\alpha }\parallel }_{\mathsf{\Omega }}^2⩽2\left(\parallel b_{\alpha }\nabla {\rho }_{\alpha }-\frac{1}{{\langle {\check{\rho }}_{\beta 0}\rangle }_{\beta }}{\langle {\widehat{\rho }}_{\beta 0}\nabla {\rho }_{\beta }\rangle }_{\beta }{\parallel }_{\mathsf{\Omega }}^2+\parallel \frac{1}{{\langle {\check{\rho }}_{\beta 0}\rangle }_{\beta }}{\langle {\widehat{\rho }}_{\beta 0}\nabla {\rho }_{\beta }\rangle }_{\beta }{\parallel }_{\mathsf{\Omega }}^2\right)\\ =\frac{2}{{\langle {\check{\rho }}_{\beta 0}\rangle }_{\beta }^2}\left(\parallel {\langle {\check{\rho }}_{\beta 0}(b_{\alpha }\nabla {\rho }_{\alpha }-b_{\beta }\nabla {\rho }_{\beta })\rangle }_{\beta }{\parallel }_{\mathsf{\Omega }}^2+{\parallel {\langle {\widehat{\rho }}_{\beta 0}\nabla {\rho }_{\beta }\rangle }_{\beta }\parallel }_{\mathsf{\Omega }}^2\right)\\ ⩽\frac{2}{{\langle {\check{\rho }}_{\beta 0}\rangle }_{\beta }^2}\left({\langle {\check{\rho }}_{\beta 0}^2\rangle }_{\beta }\langle \parallel {(b_{\alpha }\nabla {\rho }_{\alpha }-b_{\beta }\nabla {\rho }_{\beta }{\parallel }_{\mathsf{\Omega }}^2\rangle }_{\beta }+{\parallel {\langle {\widehat{\rho }}_{\beta 0}\nabla {\rho }_{\beta }\rangle }_{\beta }\parallel }_{\mathsf{\Omega }}^2\right),\end{array}$$

where ${\check{\rho }}_{\beta 0}:={\widehat{\rho }}_{\beta 0}/b_{\beta }$. According to these bounds, we derive

$$\begin{array}{c}〈\parallel \nabla {\rho }_{\alpha }{\parallel }_{\mathsf{\Omega }}^2〉⩽\frac{2}{{\langle {\check{\rho }}_{\beta 0}\rangle }_{\beta }^2{min}_{\alpha =\overline{1,K}}{b}_{\alpha }^2}\left({\langle {\check{\rho }}_{\beta 0}^2\rangle }_{\beta }\langle \parallel {(b_{\alpha }\nabla {\rho }_{\alpha }-b_{\beta }\nabla {\rho }_{\beta }{\parallel }_{\mathsf{\Omega }}^2\rangle }_{\alpha ,\beta }+K{\parallel \langle {\widehat{\rho }}_{\alpha 0}\nabla {\rho }_{\alpha }\rangle \parallel }_{\mathsf{\Omega }}^2\right)\\ ⩽{\tilde{\delta }}_1{(\parallel \nabla \mathbf{u}\parallel }_{\mathsf{\Omega }}^2+{\parallel \nabla \theta \parallel }_{\mathsf{\Omega }}^2+{\parallel \langle {\widehat{\rho }}_{\alpha 0}\nabla {\rho }_{\alpha }\rangle +({\widehat{\mathbf{u}}}_0·\nabla )\mathbf{u}+a_0\nabla \theta \parallel }_{\mathsf{\Omega }}^2\\ +{〈\parallel b_{\alpha }\nabla {\rho }_{\alpha }-b_{\beta }\nabla {\rho }_{\beta }+b_{\left(\theta \right)}({\overline{e}}_{\alpha 0}-{\overline{e}}_{\beta 0}){\nabla \theta \parallel }_{\mathsf{\Omega }}^2〉}_{\alpha ,\beta }),\end{array}$$

with ${\tilde{\delta }}_1$ depending on ${\widehat{\rho }}_{\alpha 0}$, $|{\widehat{\mathbf{u}}}_0|$, ${\overline{e}}_{\alpha 0}$ ($\alpha =\overline{1,K}$), $a_0$, $b_{\alpha }$ and $b_{\left(\theta \right)}$. This estimate and the above lower bound imply the positive definiteness inequality (50). □

Let $\mathsf{\Omega }$ be a bounded domain in ${\mathbb{R}}^n$. Define the dual space ${\mathbf{H}}^{-1}\left(\mathsf{\Omega }\right)={\left({\mathbf{H}}_0^1\left(\mathsf{\Omega }\right)\right)}^{*}$ and the duality relation ${\langle ·,·\rangle }_{\mathsf{\Omega }}$ on ${\mathbf{H}}^{-1}\left(\mathsf{\Omega }\right)\times {\mathbf{H}}_0^1\left(\mathsf{\Omega }\right)$. Denote by $\mathbf{V}\left(Q_T\right)$, with $Q_T=\mathsf{\Omega }\times (0,T)$, the space of vector functions $\tilde{\mathbf{z}}\in L^2((0,T);{\mathbf{H}}_0^1\left(\mathsf{\Omega }\right))$ possessing a distributional derivative ${\partial }_t\tilde{\mathbf{z}}\in L^2((0,T);{\mathbf{H}}^{-1}\left(\mathsf{\Omega }\right))$, see these notions, for example, in [56].

We define the weak solution $\tilde{\mathbf{z}}\in \mathbf{V}\left(Q_T\right)$, for any $T>0$, to the IBVP for the system of PDEs (40)–(42) in $Q_{\infty }$ together with conditions (43), such that the integral identity

$$\begin{array}{c}{\int }_0^T{〈{\partial }_t\tilde{\mathbf{z}}(·,t),\mathbf{z}(·,t)〉}_{\mathsf{\Omega }}dt+u_{*}{\mathcal{B}}_{Q_T}(\tilde{\mathbf{z}},\mathbf{z})+u_{*}^2{\mathcal{A}}_{Q_T}(\tilde{\mathbf{z}},\mathbf{z})+{\mathcal{A}}_{Q_T}^d(\tilde{\mathbf{z}},\mathbf{z})=0,\end{array}$$

(51)

for any $\mathbf{z}\in {\mathbf{L}}^2((0,T);{\mathbf{H}}_0^1\left(\mathsf{\Omega }\right))$ and any $T>0$, together with the initial condition $\tilde{\mathbf{z}}{|}_{t=0}={\tilde{\mathbf{z}}}^{\left(0\right)}\in {\mathbf{L}}^2\left(\mathsf{\Omega }\right)$ are valid. Here, the inner products in the bilinear forms ${\mathcal{B}}_{Q_T}$, ${\mathcal{A}}_{Q_T}$ and ${\mathcal{A}}_{Q_T}^d$ are taken over $Q_T$ instead of $\mathsf{\Omega }$ as originally. Due to the well known embedding $\mathbf{V}\left(Q_T\right)\subset C([0,T];{\mathbf{L}}^2\left(\mathsf{\Omega }\right))$ [56], the initial condition is understood by continuity in ${\mathbf{L}}^2\left(\mathsf{\Omega }\right)$. Formally, identity (51) arises from the previous one (44) for $\mathbf{z}=\mathbf{z}(·,t)$ by integration over $(0,T)$.

Now, we are ready to state the second main result.

 Theorem 2.

Let $d_{\alpha \beta 0}=d_{\beta \alpha 0}⩾d>0$ for any $\alpha \ne \beta$. The defined weak solution $\tilde{\mathbf{z}}\in \mathbf{V}\left(Q_T\right)$, for any $T>0$, to the IBVP (40)–(43) for the linearized system of PDEs exists and is unique. It satisfies the energy equality and bound

$$\begin{array}{c}\frac{1}{2}\parallel \tilde{\mathbf{z}}{(·,T)\parallel }_{{\mathbf{L}}^2\left(\mathsf{\Omega }\right)}^2+u_{*}^2{\mathcal{A}}_{Q_T}(\tilde{\mathbf{z}},\tilde{\mathbf{z}})+{\mathcal{A}}_{Q_T}^d(\tilde{\mathbf{z}},\tilde{\mathbf{z}})=\frac{1}{2}{\parallel {\tilde{\mathbf{z}}}^{\left(0\right)}\parallel }_{{\mathbf{L}}^2\left(\mathsf{\Omega }\right)}^2\forall T>0,\\ max\{\underset{t⩾0}{max}{\parallel \tilde{\mathbf{z}}(·,t)\parallel }_{{\mathbf{L}}^2\left(\mathsf{\Omega }\right)},\sqrt{2{\delta }_1}{〈\parallel \nabla {\tilde{\rho }}_{\alpha }{\parallel }_{{\mathbf{L}}^2\left(Q\right)}^2〉}^{1/2},\\ \sqrt{2}{u}_{*}{\left({\overline{\mu }}_0{\parallel \nabla \mathbf{u}\parallel }_{{\mathbf{L}}^2\left(Q\right)}^2+{\overline{\chi }}_0{\parallel div\mathbf{u}\parallel }_{L^2\left(Q\right)}^2+{\overline{\varkappa }}_0{\parallel \nabla \theta \parallel }_{{\mathbf{L}}^2\left(Q\right)}^2\right)}^{1/2}\}⩽{\parallel {\tilde{\mathbf{z}}}^{\left(0\right)}\parallel }_{{\mathbf{L}}^2\left(\mathsf{\Omega }\right)}.\end{array}$$

In addition, the derivative ${\partial }_t\left(\parallel \tilde{\mathbf{z}}{(·,t)\parallel }_{{\mathbf{L}}^2\left(\mathsf{\Omega }\right)}^2\right)\in L^1(0,\infty )$ exists, and another form of the energy equality holds

$$\begin{array}{c}\frac{1}{2}{\partial }_t(\parallel \tilde{\mathbf{z}}(·,t){\parallel }_{{\mathbf{L}}^2\left(\mathsf{\Omega }\right)}^2)+u_{*}^2{\mathcal{A}}_{\mathsf{\Omega }}(\tilde{\mathbf{z}}(·,t),\tilde{\mathbf{z}}(·,t))+{\mathcal{A}}_{\mathsf{\Omega }}^d(\tilde{\mathbf{z}}(·,t),\tilde{\mathbf{z}}(·,t))=0\end{array}$$

(52)

for almost all $t>0$ and, consequently, the following strong ${\mathbf{L}}^2\left(\mathsf{\Omega }\right)$-dissipativity property holds

$${\partial }_t\left(\parallel \tilde{\mathbf{z}}{(·,t)\parallel }_{{\mathbf{L}}^2\left(\mathsf{\Omega }\right)}^2\right)⩽0\mathrm{for}\mathrm{almost}\mathrm{all}t>0.$$

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).

Notice that the exponential decay $\parallel \tilde{\mathbf{z}}{(·,t)\parallel }_{{\mathbf{L}}^2\left(\mathsf{\Omega }\right)}⩽e^{-{\delta }_{2}t}{\parallel {\tilde{\mathbf{z}}}^{\left(0\right)}\parallel }_{{\mathbf{L}}^2\left(\mathsf{\Omega }\right)}$ for $t⩾0$, with ${\delta }_2=\frac{1}{2}({\delta }_0+{\delta }_1)>0$, follows from the energy equality (52) and Corollary 1. Methods developed in theory of linear parabolic PDEs (see, for example, references [56,57,58]) allow one to derive various regularity properties for $\tilde{\mathbf{z}}$ which we do not consider here.

Lemmas 1 and 2 and Corollary 1 lead also to other important corollaries which we briefly describe now. First, it can be checked that the system of PDEs (1), (36) and (37) is uniformly parabolic in the Petrovskii sense [58] in any bounded subdomain $\mathcal{D}\subset {\mathcal{D}}_{+}:={(0,\infty )}^K\times {\mathbb{R}}^n\times (0,\infty )$ of values of its solutions $\mathbf{z}=(\rho ,\mathbf{u},\theta )$ under additional assumptions on the diffusion fluxes $d_{\alpha \beta }=d_{\beta \alpha }>0$ and $d_{\alpha \beta },e_{\alpha }\in C^2\left({\mathcal{D}}_{+}\right)$ for all $\alpha \ne \beta$.

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 $\mathcal{D}$ positive definiteness of a matrix $A({\mathbf{z}}_0,\xi )$, defining the property is equivalent to a result such as Corollary 1 in the case $\mathsf{\Omega }={\mathbb{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({\mathbf{z}}_0,\xi )$ 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 $\ell =0$, this has recently been checked in detail in [37]. The reason is that the quadratic form ${\mathcal{A}}_{\mathsf{\Omega }}(\mathbf{z},\mathbf{z})$ is only non-negative rather than positive definite in ${\mathbf{H}}_0^1\left(\mathsf{\Omega }\right)$, and the degeneration occurs with respect to $\nabla \rho$. For $\ell =1$, this can be performed similarly and, moreover, this is clear if $\mathcal{D}$ contains a point ${\mathbf{z}}_0=({\rho }_0,0,{\theta }_0)$, since then the quadratic form ${\mathcal{A}}_{\mathsf{\Omega }}(\mathbf{z},\mathbf{z})$ is identical for $\ell =0$ and 1 at such point ${\mathbf{z}}_0$, with ${\overline{\chi }}_0+{\tau }_0{a}^2$ substituted for ${\overline{\chi }}_0$. For systems of composite hyperbolic–parabolic 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.

4. 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

$$\begin{array}{c}{\partial }_t{\rho }_{\alpha }+{\partial }_x\left({\rho }_{\alpha }(u-w_{\ell \alpha })+d_{\alpha }\right)=0,\alpha =\overline{1,K},\end{array}$$

(53)

$$\begin{array}{c}{\partial }_t\left(\rho u\right)+{\partial }_x\left(\langle {\rho }_{\alpha }(u-w_{\ell \alpha })\rangle u+p\right)={\partial }_x{\mathrm{\Pi }}_{\ell }+\left(\rho -\ell \tau {\partial }_x\left(\rho u\right)\right)f,\end{array}$$

(54)

$$\begin{array}{c}{\partial }_{t}E+{\partial }_x\langle (E_{\alpha }+p_{\alpha })(u-w_{\ell \alpha })\rangle ={\partial }_x(-q_{\ell }+{\mathrm{\Pi }}_{\ell }u)+\langle {\rho }_{\alpha }(u-w_{\ell \alpha })\rangle f+Q,\end{array}$$

(55)

with $\ell =0,1$. The main sought functions are the component densities ${\rho }_1>0,\dots ,{\rho }_K>0$ and their common velocity and temperature u and $\theta >0$ depending on $(x,t)\in [-X,X]\times [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 ${\left|\mathbf{u}\right|}^2=u^2$.

The 1D formulas for the regularizing velocities, viscous stress and heat flux look as follows

$$\begin{array}{c}{w}_{\ell \alpha }=\ell \frac{\tau }{{\rho }_{\alpha }}u{\partial }_x\left({\rho }_{\alpha }u\right)+\widehat{w},\widehat{w}=\frac{\tau }{\rho }(\rho u{\partial }_{x}u+{\partial }_{x}p-\rho f),\end{array}$$

(56)

$$\begin{array}{c}{\mathrm{\Pi }}_{\ell }=\nu {\partial }_{x}u+{\mathrm{\Pi }}_{\ell }^{\tau },{\mathrm{\Pi }}_{\ell }^{\tau }=\rho u\widehat{w}+\ell \tau \{u{\partial }_{x}p+\langle {\gamma }_{\alpha }{p}_{\alpha }\rangle {\partial }_{x}u-\langle ({\gamma }_{\alpha }-1)Q_{\alpha }\rangle \},\end{array}$$

(57)

$$\begin{array}{c}{q}_{\ell }=-\varkappa {\partial }_x\theta +\ell q^{\tau }+q^d,\end{array}$$

(58)

$$\begin{array}{c}-q^{\tau }=\tau \left\{\left(\langle c_{V\alpha }{\rho }_{\alpha }\rangle {\partial }_x\theta -\theta {\partial }_x\langle R_{\alpha }{\rho }_{\alpha }\rangle \right)u^2-Qu\right\}=\tau \left\{\left(c_V\rho {\partial }_x\theta -\theta {\partial }_x\left(R\rho \right)\right)u^2-Qu\right\},\end{array}$$

(59)

$$\begin{array}{c}-d_{\alpha }={〈d_{\alpha \beta }\left({\partial }_x(G_{\alpha }-G_{\beta })+(e_{\alpha }-e_{\beta }){\partial }_x\theta \right)〉}_{\beta },q^d=〈\left(G_{\alpha }+e_{\alpha }\theta \right)d_{\alpha }〉,\end{array}$$

(60)

where $\nu :=\frac{4}{3}\mu +\lambda$ and $\alpha =\overline{1,K}$.

As above, introducing the average regularizing velocity $w_{\ell }:=〈\frac{{\rho }_{\alpha }}{\rho }{w}_{\ell \alpha }〉=\ell \frac{\tau }{\rho }u{\partial }_x\left(\rho u\right)+\widehat{w}$ allows one to simplify the form of the balance PDEs for the momentum and total energy of the mixture (54) and (55):

$$\begin{array}{c}{\partial }_t\left(\rho u\right)+{\partial }_x\left(\rho (u-w_{\ell })u\right)+{\partial }_{x}p={\partial }_x{\mathrm{\Pi }}_{\ell }+\left(\rho -\ell \tau {\partial }_x\left(\rho u\right)\right)f,\end{array}$$

(61)

$$\begin{array}{c}{\partial }_{t}E+{\partial }_x\left(0.5\rho u^2(u-w_{\ell })+\langle {\rho }_{\alpha }{c}_{p\alpha }\theta (u-w_{\ell \alpha })\rangle \right)={\partial }_x(-q_{\ell }+{\mathrm{\Pi }}_{\ell }u)+\rho (u-w_{\ell })f+Q.\end{array}$$

(62)

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 ${\overline{\omega }}_h$ on $[-X,X]$, with the nodes $x_i=-X+ih$, $0⩽i⩽N$, and the step $h=\frac{2X}{N}$. Let ${\omega }_h={\overline{\omega }}_h\setminus \{-X,X\}$ be its internal part. Define also an auxuliary mesh ${\omega }_h^{*}$ with the nodes $x_{i+1/2}=(i+1/2)h$, $0⩽i⩽N-1$.

Let $H\left(\omega \right)$ be the space of functions defined on a mesh $\omega$. 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

$$\begin{array}{c}{\left[v\right]}_{i+1/2}=0.5(v_i+v_{i+1}),\delta v_{i+1/2}=\frac{v_{i+1}-v_i}{h},{\left[y\right]}_i^{*}=0.5(y_{i-1/2}+y_{i+1/2}),{\delta }^{*}{y}_i=\frac{y_{i+1/2}-y_{i-1/2}}{h}\end{array}$$

on functions $v\in H\left({\overline{\omega }}_h\right)$ and $y\in H\left({\omega }_h^{*}\right)$, where $v_i=v\left(x_i\right)$ and $y_{i+1/2}=y\left(x_{i+1/2}\right)$.

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

$$\begin{array}{c}{\partial }_t{\rho }_{\alpha }+{\delta }^{*}\left({\left[{\rho }_{\alpha }\right]}_{ln}(\left[u\right]-w_{\ell \alpha })+d_{\alpha }\right)=0,\alpha =\overline{1,K},\end{array}$$

(63)

$$\begin{array}{c}{\partial }_t\left(\rho u\right)+{\delta }^{*}\left(\langle {\left[{\rho }_{\alpha }\right]}_{ln}(\left[u\right]-w_{\ell \alpha })\rangle \left[u\right]+\left[p\right]\right)={\delta }^{*}{\mathrm{\Pi }}_{\ell },\end{array}$$

(64)

$$\begin{array}{c}{\partial }_{t}E+{\delta }^{*}\left\{〈({\left[E_{\alpha }\right]}_2+\left[p_{\alpha }\right])(\left[u\right]-w_{\ell \alpha })〉-0.25h^2\delta u·\delta p\right\}={\delta }^{*}(-q_{\ell }+{\mathrm{\Pi }}_{\ell }\left[u\right])+{\left[Q\right]}^{*}\end{array}$$

(65)

on ${\omega }_h\times [0,T]$. The main sought functions ${\rho }_1>0,\dots ,{\rho }_K>0$, u and $\theta >0$ together with the functions $p_{\alpha }$, ${\epsilon }_{\alpha }$ and $E_{\alpha }$ are defined in space on the main mesh ${\overline{\omega }}_h$. In the equations, the above expressions (4) for $p_{\alpha }$, ${\epsilon }_{\alpha }$ and $E_{\alpha }$ as well as (5) for $\rho$, p, $\epsilon$ and E, with ${\left|\mathbf{u}\right|}^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):

$$\begin{array}{c}{w}_{\ell \alpha }=\ell \frac{\tau }{\left[{\rho }_{\alpha }\right]}\left[u\right]\delta \left({\rho }_{\alpha }u\right)+\widehat{w},\alpha =\overline{1,K},\widehat{w}=\frac{\tau }{\left[\rho \right]}(\left[\rho \right]\left[u\right]\delta u+\delta p),\end{array}$$

(66)

as well as of the viscous stress, and heat flux and diffusion fluxes (57)–(60):

$$\begin{array}{c}{\mathrm{\Pi }}_{\ell }=\nu \delta u+{\mathrm{\Pi }}_{\ell }^{\tau },{\mathrm{\Pi }}_{\ell }^{\tau }=\left[u\right]\left[\rho \right]\widehat{w}+\ell \tau \left\{\left[u\right]\delta p+\langle {\gamma }_{\alpha }{\left[p_{\alpha }\right]}_1\rangle \delta u-\langle {\gamma }_{\alpha }{Q}_{\alpha }\rangle +Q\right\},\end{array}$$

(67)

$$\begin{array}{c}{q}_{\ell }=-\varkappa \delta \theta +q^d+\ell q^{\tau },\end{array}$$

(68)

$$\begin{array}{c}-q^{\tau }=\tau \left\{\left(\langle c_{V\alpha }\left[{\rho }_{\alpha }\right]\rangle \delta \theta -\left[\theta \right]\delta \langle R_{\alpha }{\rho }_{\alpha }\rangle \right){\left[u\right]}^2-Q\left[u\right]\right\}=\tau \left\{\left(\left[c_V\rho \right]\delta \theta -\left[\theta \right]\delta \left(R\rho \right)\right){\left[u\right]}^2-Q\left[u\right]\right\},\end{array}$$

(69)

$$\begin{array}{c}-d_{\alpha }={〈d_{\alpha \beta }\left(\delta (G_{\alpha }-G_{\beta })+(e_{\alpha }-e_{\beta })\delta \theta \right)〉}_{\beta },q^d=〈\left(\left[G_{\alpha }\right]+e_{\alpha }\left[\theta \right]\right)d_{\alpha }〉.\end{array}$$

(70)

The functions $w_{\ell \alpha }$, $\widehat{w}$, ${\mathrm{\Pi }}_{\ell \alpha }$, $q_{\ell \alpha }$, $\tau$, ${\nu }_{\alpha }$, ${\varkappa }_{\alpha }$ and $Q_{\alpha }$ are defined in space on the auxiliary mesh ${\omega }_h^{*}$. Moreover, $G_{\alpha }$, the Gibbs potential of the component $\alpha$, see Formula (15), is defined in space on ${\overline{\omega }}_h$, whereas the functions $d_{\alpha },e_{\alpha }$ and $q^d$ are defined in space on ${\omega }_h^{*}$.

Here, we apply nonstandard averages of ${\rho }_{\alpha }$, $p_{\alpha }$, $E_{\alpha }$ and ${\epsilon }_{\alpha }$ of the form [41,47]

$$\begin{array}{c}{\left[{\rho }_{\alpha }\right]}_{ln}=\frac{1}{ln({\rho }_{\alpha -};{\rho }_{\alpha +})},{\left[p_{\alpha }\right]}_1=R_{\alpha }\left[{\rho }_{\alpha }\right]\left[{\theta }_{\alpha }\right],{\left[E_{\alpha }\right]}_2=0.5{\left[{\rho }_{\alpha }\right]}_{ln}{u}_{-}{u}_{+}+{\left[{\rho }_{\alpha }\right]}_{ln}{\left[{\epsilon }_{\alpha }\right]}^{ln},\\ [{\epsilon }_{\alpha }{]}^{ln}=c_{V\alpha }{\left[\theta \right]}^{ln},{\left[\theta \right]}^{ln}:=ln\left({\frac{1}{\theta }}_{-};{\frac{1}{\theta }}_{+}\right)={\theta }_{-}{\theta }_{+}ln({\theta }_{-};{\theta }_{+})\end{array}$$

that exploit the divided difference for the logarithmic function

$$ln(a;b)=\frac{lnb-lna}{b-a}\mathrm{for}a\ne b,ln(a;a)=\frac{1}{a},a>0,b>0.$$

Consequently, we have $\langle {\gamma }_{\alpha }{\left[p_{\alpha }\right]}_1\rangle =\langle \left[{\rho }_{\alpha }\right]R_{\alpha }{\gamma }_{\alpha }\rangle \left[\theta \right]$ in expression (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 $\tau =\mathcal{T}\left(\rho ,u,\theta \right)$, one can set, in particular, $\tau =\mathcal{T}\left(\left[\rho \right],\left[u\right],\left[\theta \right]\right)$ or $\tau =\left[\mathcal{T}(\rho ,u,\theta )\right]$ in space on ${\omega }_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_{\ell \alpha }$ (approximating formulas (9), not (7), in the 1D case) and the much more general expression (70) for $d_{\alpha }$ 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

$$\begin{array}{c}\delta G_{\alpha }=-\left[\theta \right]\delta s_{\alpha }+(c_{p\alpha }-\left[s_{\alpha }\right])\delta \theta ,\delta s_{\alpha }=-R_{\alpha }ln({\rho }_{\alpha -};{\rho }_{\alpha +})\delta {\rho }_{\alpha }+c_{V\alpha }ln({\theta }_{-};{\theta }_{+})\delta \theta ,\end{array}$$

we obtain $\left[G_{\alpha }\right]+e_{\alpha }\left[\theta \right]=(c_{p\alpha }+e_{\alpha })\left[\theta \right]-\left[s_{\alpha }\theta \right]$ and

$$\begin{array}{c}\delta G_{\alpha }+e_{\alpha }\delta \theta =R_{\alpha }\frac{\left[\theta \right]}{{\left[{\rho }_{\alpha }\right]}_{ln}}\delta {\rho }_{\alpha }+\left(c_{p\alpha }-c_{V\alpha }\frac{\left[\theta \right]}{{\left[\theta \right]}_{ln}}-\left[s_{\alpha }\right]+e_{\alpha }\right)\delta \theta .\end{array}$$

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_{\alpha }$, for $K⩾3$; this form is applicable in [41] as well.

 Theorem 3.

Let $d_{\alpha \beta }=d_{\beta \alpha }⩾0$ for any $\alpha \ne \beta$. For the 1D semi-discrete method (63)–(70), the balance equation for the mixture entropy holds

$$\begin{array}{c}{\partial }_t\left(\rho s\right)+{\delta }^{*}〈j_{\ell \alpha }\left[s_{\alpha }\right]〉={\delta }^{*}\left((\varkappa \delta \theta -\ell q^{\tau })\left[\frac{1}{\theta }\right]-\langle e_{\alpha }{d}_{\alpha }\rangle \frac{{\left[\theta \right]}^2}{{\theta }_{-}{\theta }_{+}}+B_h^{\left(d\right)}\right)+{\left[{\mathcal{P}}_h^{NS}+{\mathcal{P}}_h^{\tau }\right]}^{*}\end{array}$$

(71)

on ${\omega }_h\times [0,T]$, with the component mass fluxes $j_{\ell \alpha }={\left[{\rho }_{\alpha }\right]}_{ln}(\left[u\right]-w_{\ell \alpha })$ and the terms

$$\begin{array}{c}{B}_h^{\left(d\right)}:=〈R_{\alpha }{j}_{\ell \alpha }\left(1-\frac{\left[{\rho }_{\alpha }\right]}{{\left[{\rho }_{\alpha }\right]}_{ln}}\right)+c_{V\alpha }{j}_{\ell \alpha }\left(1-{\left[{\epsilon }_{\alpha }\right]}^{ln}\left[\frac{1}{{\epsilon }_{\alpha }}\right]\right)〉\\ -0.25h^2\left({\mathrm{\Pi }}_{\ell }\delta u-\langle d_{\alpha }\delta G_{\alpha }\rangle +\langle w_{\ell \alpha }\delta p_{\alpha }\rangle +Q\right)\delta \frac{1}{\theta },\\ {\mathcal{P}}_h^{NS}:=\frac{1}{{\theta }_{-}{\theta }_{+}}\left\{\varkappa {\left(\delta \theta \right)}^2+\nu \left[\theta \right]{\left(\delta u\right)}^2+0.5\left[\theta \right]{〈d_{\alpha \beta }{\left(\delta G_{\alpha }+e_{\alpha }\delta \theta -(\delta G_{\beta }+e_{\beta }\delta \theta )\right)}^2〉}_{\alpha ,\beta }\right\}⩾0,\\ {\mathcal{P}}_h^{\tau }:=\frac{1}{{\theta }_{-}{\theta }_{+}}\{\left[\theta \right]\left[\rho \right]\frac{1}{\tau }{\widehat{w}}^2+\ell \tau {\left[\theta \right]}^2〈\frac{R_{\alpha }}{\left[{\rho }_{\alpha }\right]}{\left(\delta \left({\rho }_{\alpha }u\right)\right)}^2〉\\ +\ell \tau 〈c_{V\alpha }\left[{\rho }_{\alpha }\right]{\left(\left[u\right]\delta \theta +({\gamma }_{\alpha }-1)\left[\theta \right]\delta u-\frac{Q_{\alpha }}{2c_{V\alpha }\left[{\rho }_{\alpha }\right]}\right)}^2〉+\left[\theta \right]〈Q_{\alpha }\left(1-\ell \frac{\tau ({\gamma }_{\alpha }-1)Q_{\alpha }}{4{\left[p_{\alpha }\right]}_1}\right)〉\}.\end{array}$$

The term ${\left[{\mathcal{P}}_h^{NS}+{\mathcal{P}}_h^{\tau }\right]}^{*}$ in Equation (71) is the semi-discrete entropy production. The first three terms of ${\mathcal{P}}_h^{\tau }$ are non-negative, and the last term is non-negative for $\ell =0$, as well as for $\ell =1$ under the condition $\tau 〈\frac{({\gamma }_{\alpha }-1)Q_{\alpha }^2}{4{\left[p_{\alpha }\right]}_1}〉⩽Q.$ This condition is certainly true provided that $\tau ({\gamma }_{\alpha }-1)Q_{\alpha }⩽4{\left[p_{\alpha }\right]}_1$, $\alpha =\overline{1,K}$.

 Proof. 

The following semi-discrete balance equations for the total mass and kinetic and internal energies of the mixture hold

$$\begin{array}{c}{\partial }_t\rho +{\delta }^{*}{j}_{\ell }=0,0.5{\partial }_t\left(\rho u^2\right)+0.5{\delta }^{*}\left(j_{\ell }{u}_{-}{u}_{+}\right)+\left({\delta }^{*}\left[p\right]\right)u=\left({\delta }^{*}{\mathrm{\Pi }}_{\ell }\right)u,\\ {\partial }_t\left(\rho \epsilon \right)+{\delta }^{*}\langle j_{\ell \alpha }{\left[{\epsilon }_{\alpha }\right]}^{ln}\rangle =-{\delta }^{*}{q}_{\ell }-\langle p_{\alpha }{\delta }^{*}(\left[u\right]-w_{\ell \alpha })\rangle +{[{\mathrm{\Pi }}_{\ell }\delta u+\langle w_{\ell \alpha }\delta p_{\alpha }\rangle +Q]}^{*}\end{array}$$

(72)

on ${\omega }_h\times [0,T]$, with $j_{\ell }:=\langle j_{\ell \alpha }\rangle$. 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_{\ell \alpha }$, ${\mathrm{\Pi }}_{\ell }$ and $q_{\ell }$, 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

$$\begin{array}{c}{\partial }_t\left(\rho s\right)+{\delta }^{*}〈j_{\ell \alpha }\left[s_{\alpha }\right]〉={\delta }^{*}\left((\varkappa \delta \theta -\ell q^{\tau })\left[\frac{1}{\theta }\right]-\langle e_{\alpha }{d}_{\alpha }\rangle \frac{{\left[\theta \right]}^2}{{\theta }_{-}{\theta }_{+}}+B_h^{\left(d\right)}\right)\\ +{\left[\frac{\varkappa {\left(\delta \theta \right)}^2}{{\theta }_{-}{\theta }_{+}}+\frac{\nu \left[\theta \right]{\left(\delta u\right)}^2}{{\theta }_{-}{\theta }_{+}}+q^d\delta \frac{1}{\theta }-〈d_{\alpha }\delta \frac{G_{\alpha }}{\theta }〉+\frac{\langle A_{\ell \alpha }\rangle }{{\theta }_{-}{\theta }_{+}}\right]}^{*},\end{array}$$

(73)

where we have taken into account that $e_{\alpha }$ plays the role of $K^{-1}{b}_{\alpha }$ 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_{\ell \alpha }$ does not matter in this derivation. Here, the following term and its decomposition

$$\begin{array}{c}{A}_{\ell }=-\ell q^{\tau }\delta \theta +\left({\mathrm{\Pi }}_{\ell }^{\tau }\delta u+\langle w_{\ell \alpha }\delta p_{\alpha }\rangle +Q\right)\left[\theta \right]=\langle A_{\ell \alpha }\rangle \end{array}$$

are involved, with

$$\begin{array}{c}{A}_{\ell \alpha }=-\ell q_{\alpha }^{\tau }\delta \theta +\left({\mathrm{\Pi }}_{\ell \alpha }^{\tau }\delta u+w_{\ell \alpha }\delta p_{\alpha }+Q_{\alpha }\right)\left[\theta \right],-q_{\alpha }^{\tau }=\tau \left\{{\left[u\right]}^2\left(c_{V\alpha }\left[{\rho }_{\alpha }\right]\delta \theta -R_{\alpha }\left[\theta \right]\delta {\rho }_{\alpha }\right)-Q_{\alpha }\left[u\right]\right\},\\ {\mathrm{\Pi }}_{\ell \alpha }^{\tau }=\left[u\right]\left[{\rho }_{\alpha }\right]\widehat{w}+\ell \tau \left\{\left[u\right]\delta p_{\alpha }+{\gamma }_{\alpha }{\left[p_{\alpha }\right]}_1\delta u-({\gamma }_{\alpha }-1)Q_{\alpha }\right\},\end{array}$$

where only the term $\left[u\right]\left[{\rho }_{\alpha }\right]\widehat{w}$ is written differently, but the corresponding average $\langle \left[u\right]\left[{\rho }_{\alpha }\right]\widehat{w}\rangle =\left[u\right]\left[\rho \right]\widehat{w}$ is the same term of ${\mathrm{\Pi }}_{\ell }^{\tau }$ in the above expression for $A_{\ell }$.

Let us transform the difference of the third and fourth terms under the sign ${[·]}^{*}$ on the right in Equation (73). We need the elementary formulas

$$\begin{array}{c}\delta \frac{G_{\alpha }}{\theta }=\left[G_{\alpha }\right]\delta \frac{1}{\theta }+\left(\delta G_{\alpha }\right)\left[\frac{1}{\theta }\right],\delta \frac{1}{\theta }=-\frac{\delta \theta }{{\theta }_{-}{\theta }_{+}},\left[\frac{1}{\theta }\right]=\frac{\left[\theta \right]}{{\theta }_{-}{\theta }_{+}}.\end{array}$$

Applying them and then identity (16), we obtain

$$\begin{array}{c}{q}^d\delta \frac{1}{\theta }-〈d_{\alpha }\delta \frac{G_{\alpha }}{\theta }〉=\left(q^d-\langle d_{\alpha }\left[G_{\alpha }\right]\rangle \right)\delta \frac{1}{\theta }-\langle d_{\alpha }\delta G_{\alpha }\rangle \left[\frac{1}{\theta }\right]\\ =-\frac{\left[\theta \right]}{{\theta }_{-}{\theta }_{+}}〈d_{\alpha }(e_{\alpha }\delta \theta +\delta G_{\alpha })〉=\frac{\left[\theta \right]}{2{\theta }_{-}{\theta }_{+}}{〈d_{\alpha \beta }{\left(\delta G_{\alpha }+e_{\alpha }\delta \theta -(\delta G_{\beta }+e_{\beta }\delta \theta )\right)}^2〉}_{\alpha ,\beta }.\end{array}$$

Next, using expressions (66), we can extract from $A_{\ell \alpha }$ the term $\left[\theta \right]A_{\ell \alpha }^{\prime }$ such that

$$\begin{array}{c}{A}_{\ell \alpha }^{\prime }:=\left[{\rho }_{\alpha }\right]\left[u\right]\widehat{w}\delta u+w_{\ell \alpha }\delta p_{\alpha }=(\left[{\rho }_{\alpha }\right]\left[u\right]\delta u+\delta p_{\alpha })\widehat{w}+\ell \frac{\tau }{\left[{\rho }_{\alpha }\right]}\left[u\right]\left(\delta \left({\rho }_{\alpha }u\right)\right)\delta p_{\alpha },\end{array}$$

and, according to ([41], Appendix A) or [47], the following formula holds

$$\begin{array}{c}{A}_{\ell \alpha }-\left[\theta \right](\left[{\rho }_{\alpha }\right]\left[u\right]\delta u+\delta p_{\alpha })\widehat{w}=\ell \tau {\left[\theta \right]}^2\frac{R_{\alpha }}{\left[{\rho }_{\alpha }\right]}{\left(\delta \left({\rho }_{\alpha }u\right)\right)}^2\\ +\ell \tau c_{V\alpha }\left[{\rho }_{\alpha }\right]{\left(\left[u\right]\delta \theta +({\gamma }_{\alpha }-1)\left[\theta \right]\delta u-\frac{Q_{\alpha }}{2c_{V\alpha }\left[{\rho }_{\alpha }\right]}\right)}^2+\left[\theta \right]Q_{\alpha }\left(1-\ell \frac{\tau ({\gamma }_{\alpha }-1)Q_{\alpha }}{4{\left[p_{\alpha }\right]}_1}\right).\end{array}$$

Applying the operation $\langle ·\rangle$ to it and accomplishing the transformations

$$\begin{array}{c}〈\left[\theta \right](\left[{\rho }_{\alpha }\right]\left[u\right]\delta u+\delta p_{\alpha })\widehat{w}〉=\left[\theta \right](\left[\rho \right]\left[u\right]\delta u+\delta p)\widehat{w}=\left[\theta \right]\left[\rho \right]\frac{1}{\tau }{\widehat{w}}^2,\end{array}$$

(74)

we complete the proof. □

As in the differential case, the entropy production remains non-negative for $\tau ⩾0$, where one should pass to another form for the first relaxation term in ${\mathcal{P}}_h^{\tau }$ inside the curly brackets: $\left[\theta \right]\left[\rho \right]\frac{1}{\tau }{\widehat{w}}^2=\tau \frac{\left[\theta \right]}{\left[\rho \right]}{(\left[\rho \right]\left[u\right]\delta u+\delta p)}^2.$

At the end of the section, following [41,47], we generalize the constructed semi-discrete 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

$${\left[{\rho }_{*\ell }f\right]}^{*},〈{\left[\left[{\rho }_{\alpha }\right](\left[u\right]-w_{\ell \alpha })f\right]}^{*}〉+0.25h^2{\left[{\rho }_{*\ell }\left(\delta u\right)\left(\delta \frac{1}{\theta }\right)f\right]}^{*}\theta$$

to their right-hand sides, with the functions ${\rho }_{*\ell }:=\left[\rho \right]-\ell \tau \delta \left(\rho u\right)$ and f defined in space on ${\omega }_h^{*}$. We also generalize the expression for $\widehat{w}$ as $\widehat{w}=\frac{\tau }{\left[\rho \right]}(\left[\rho \right]\left[u\right]\delta u+\delta p-\left[\rho \right]f)$.

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):

$$\mathsf{\Psi }:=〈{\left[\left[{\rho }_{\alpha }\right](\left[u\right]-w_{\ell \alpha })f\right]}^{*}〉-{\left[{\rho }_{*\ell }f\right]}^{*}u+0.25h^2{\left[{\rho }_{*\ell }\left(\delta u\right)\left(\delta \frac{1}{\theta }\right)f\right]}^{*}\theta .$$

To derive the semi-discrete balance equation for the mixture entropy (71), one should multiply this term by $\frac{1}{\theta }$ and transform the result. The required transformation was accomplished in [41] for general $w_{\ell \alpha }=\ell \frac{\tau }{\left[{\rho }_{\alpha }\right]}\left[u\right]\delta \left({\rho }_{\alpha }u\right)+{\widehat{w}}_{\alpha }$ with any ${\widehat{w}}_{\alpha }$, and in our case where ${\widehat{w}}_{\alpha }=\widehat{w}$ is independent of $\alpha$, it leads to the formulas

$$\begin{array}{c}\frac{\mathsf{\Psi }}{\theta }=-{\left[\langle \left[{\rho }_{\alpha }\right]\rangle \widehat{w}f\left[\frac{1}{\theta }\right]\right]}^{*}+{\delta }^{*}{C}_h,C_h:=0.25h^2\left(\left[\rho \right]\widehat{w}\delta \frac{1}{\theta }+{\rho }_{*\ell }\left(\delta u\right)\left[\frac{1}{\theta }\right]\right)f.\end{array}$$

As a result, in the preliminary entropy balance equation (73), the additional term $-\left[\theta \right]\left[{\rho }_{\alpha }\right]\widehat{w}f$ appears in $A_{\ell \alpha }$, and the term $C_h$ should be added to $B_h^{\left(d\right)}$. Thus, Formula (74) takes the form

$$〈\left[\theta \right](\left[{\rho }_{\alpha }\right]\left[u\right]\delta u+\delta p_{\alpha }-\left[{\rho }_{\alpha }\right]f)\widehat{w}〉=\left[\theta \right](\left[\rho \right]\left[u\right]\delta u+\delta p-\left[\rho \right]f)\widehat{w}=\left[\theta \right]\left[\rho \right]\frac{1}{\tau }{\widehat{w}}^2.$$

With the given f-dependent extensions, Theorem 3 remains valid.

5. 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_{\ell \alpha }$ in the balance PDEs (53)–(55):

$$\begin{array}{c}{w}_{\ell \alpha }=\ell \frac{\tau }{{\rho }_{\alpha }}u{\partial }_x\left({\rho }_{\alpha }u\right)+\widehat{w},\widehat{w}=\tau \left(u{\partial }_{x}u+\frac{1}{\rho }{\partial }_{x}p\right)\left(A\right);\\ w_{\ell \alpha }=\ell \frac{\tau }{{\rho }_{\alpha }}u{\partial }_x\left({\rho }_{\alpha }u\right)+{\widehat{w}}_{\alpha },{\widehat{w}}_{\alpha }=\tau \left(u{\partial }_{x}u+\frac{1}{{\rho }_{\alpha }}{\partial }_x{p}_{\alpha }\right)\left(B\right);w_{\ell \alpha }=\frac{\tau }{\rho }u{\partial }_x\left(\rho u\right)+\widehat{w}\left(C\right),\end{array}$$

with $\ell =1$ and $\alpha =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_{\ell \alpha }$ is independent of ℓ and $\alpha$. We discretize these expressions, respectively, according to Formula (66) as well as

$$\begin{array}{c}{w}_{\ell \alpha }=\ell \frac{\tau }{\left[{\rho }_{\alpha }\right]}\left[u\right]\delta \left({\rho }_{\alpha }u\right)+{\widehat{w}}_{\alpha },{\widehat{w}}_{\alpha }=\tau \left(\left[u\right]\delta u+\frac{1}{\left[{\rho }_{\alpha }\right]}\delta p_{\alpha }\right)\mathrm{and}{w}_{\ell \alpha }=\frac{\tau }{\left[\rho \right]}\left[u\right]\delta \left(\rho u\right)+\widehat{w}.\end{array}$$

We consider four test examples known in the literature. In Examples 1–3, we take the following piecewise constant initial data $({\rho }_1,{\rho }_2,p,u){|}_{t=0}=({\rho }_1^0,{\rho }_2^0,p^0,u^0)$ (a Riemann problem) and piecewise constant physical parameters

$$\begin{array}{c}({\rho }_1^0,{\rho }_2^0,p^0,u^0,\gamma ,c_V)\left(x\right)=\left\{\begin{array}{cc}({\rho }_{1l},{\rho }_{2l},p_l,u_l,{\gamma }_1,c_{V1}),\hfill & -X⩽x<0\hfill \\ ({\rho }_{1r},{\rho }_{2r},p_r,u_r,{\gamma }_2,c_{V2}),\hfill & 0⩽x⩽X\hfill \end{array}\right..\end{array}$$

Moreover, ${\rho }_{1r}={\rho }_{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_{fin}$ can be found in

Table 1. The initial parameters to the left and right of the discontinuity between two gases and the final time of computations.

Example	$\mathit{\rho }$	p	u	$\mathit{\gamma }$	${\mathit{c}}_{\mathit{V}}$	${\mathit{t}}_{\mathbf{fin}}$
(1) left	0.138	1	0.5	1.67	3.11	0.2
(1) right	1	1	0.5	1.4	0.72
(2) left	14.54903	1.943$\times {10}^7$	0	$5/3$	2420	0.0002
(2) right	1.16355	${10}^5$	0	1.4	732
(3) left	1	500	0	1.4	1	0.011
(3) right	1	0.2	0	1.6	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 $\tau =\nu =\varkappa =0$, there exists a closed Euler-type system of PDEs for the sought functions $\rho$, $\gamma$, u and $\theta$, see details in Appendix A), and we have chosen them ourselves. The initial temperature ${\theta }^0$ is defined in accordance with Formulas (4) and (5): $p^0=R{\rho }^0{\theta }^0=(({\gamma }_1-1)c_{V1}{\rho }_1+({\gamma }_2-1)c_{V2}{\rho }_2){\theta }^0.$ The boundary values of the sought functions at $x=\pm 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: $\frac{p_l}{p_r}=1,194.3,2500.$

Recall that in order to avoid loss of accuracy in computation of ${\left[{\rho }_{\alpha }\right]}_{ln}$ and ${\left[{\epsilon }_{\alpha }\right]}^{ln}$, one can apply the trapezoidal or midpoint rule to the integral representation of $ln(a;b)$ in the case $\frac{b}{a}\approx 1$:

$$ln(a;b)={\int }_0^1\frac{1}{(1-r)a+br}dr\approx \frac{1}{2a}+\frac{1}{2b},\frac{2}{a+b}.$$

We apply them for $\epsilon$ and ${\rho }_{\alpha }$, respectively, that leads to the formulas ${\left[{\rho }_{\alpha }\right]}_{ln}\approx \left[{\rho }_{\alpha }\right]$ and ${\left[{\epsilon }_{\alpha }\right]}^{ln}\approx \left[{\epsilon }_{\alpha }\right]$.

We introduce a non-uniform mesh in time $0=t_0<t_1<\dots <t_{\overline{m}}=t_{fin}$, 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 $\nu =\frac{4}{3}\mu +\lambda$:

$${\tau }_i^m=\frac{ah}{c_{si}^m+i_{\tau }\left|u_i^m\right|},\nu =\tau (a_{S1}{p}_1+a_{S2}{p}_2),\varkappa =\tau a_{Pr}({\gamma }_1{c}_{V1}{p}_1+{\gamma }_2{c}_{V2}{p}_2).$$

Here, $0<a<1$ is a parameter, $h=\frac{2X}{N}$, $i_{\tau }=0$ in Examples 1 and 3 or $i_{\tau }=1$ in Examples 2 and 4, and, for example, $u_i^m=u(x_i,t_m)$. In addition, $a_{S1}=a_{S2}=\frac{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

$$h_{tm}=\frac{\beta h}{{max}_i(c_{si}^{m-1}+|u_i^{m-1}\left|\right)},1⩽m⩽\overline{m}-1,h_{t\overline{m}}=t_{fin}-t_{\overline{m}-1}⩽\frac{\beta h}{{max}_i\left(c_{si}^{\overline{m}-1}+\left|u_i^{\overline{m}-1}\right|\right)},$$

where $\beta$ 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 $\beta$. If the values of ${\rho }_1$ or ${\rho }_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$, $\beta =0.7$, $N=251$ (similarly to [54]) and 1001 are given in Figure 1; note that the scales for p, u and $\theta$ are enlarged there. In this and other examples, we exhibit graphs of six functions: ${\rho }_1$, ${\rho }_2$, $\rho$, p, u and $\theta$. The deviations from constant values near the contact discontinuity are very small in ${\rho }_1$ and p, slightly larger for u and $\theta$. They diminish as N grows and, for $N=1001$, they disappear for p, almost disappear in ${\rho }_1$ and u and become very small in $\theta$. 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 $\beta$ to obtain suitable results: $a=0.2$ and $\beta =0.1$. However, for $N=251$, the results are still not so nice: the behaviour of ${\rho }_1$, ${\rho }_2$ and $\rho$ is too smooth near the contact discontinuity, the deviations from constant values near the contact discontinuity are very small for p, small for $\theta$, but rather large in u. For $N=1001$, the results become better, see Figure 2.

In case C, the results for $a=0.2$, $\beta =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).

 Example 2.

(a version of the Sod problem) from ([52], Table I). In the original paper, $t_{fin}$ 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 ${\rho }_1$, ${\rho }_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 ${\rho }_1$ and p are non-increasing, whereas ${\rho }_2$, ρ, u and θ are non-monotone, with the maximal values of ${\rho }_2$, u and θ in front of the shock. In addition, ${\rho }_2$ is piecewise constant. The final maximal Mach number is $M_{max}:={max}_i\frac{|u_i^{\overline{m}}|}{c_{si}^{\overline{m}}}\approx 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 $\beta =0.4$. Notice that, in this example, we can take zero artificial viscosity $\nu =0$ without essential changing the results.

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 $\rho$ at the contact discontinuity. Notice that scaling u and $\theta$ by the natural divisors $u_{*}={\theta }_{*}=1000$ (with no scaling of ${\rho }_1$, ${\rho }_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 ${\rho }_2$ is slightly more smoothed and the graph of $\rho$ 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 ${\rho }_1$, p and $\theta$ are computed rather accurately, the graphs of ${\rho }_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 $\rho$ 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 ${\rho }_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 ${\rho }_2$ is slightly more smoothed and the graph of $\rho$ has a slightly deeper hollow near the contact discontinuity. In case C, the results are already poor in general: though ${\rho }_1$, p and $\theta$ are computed rather accurately, the graphs of ${\rho }_2$ and u have visible “fingers” at the point of contact discontinuity, and the graph of $\rho$ 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_{fin}$ 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 ${\rho }_2$ is narrower, θ is non-increasing and the jumps in the values of ${\rho }_2$, ρ and θ are high. Furthermore, $M_{max}\approx 1.44$, so the flow is partially supersonic once again. We take $a=0.25$ and $\beta =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 ${\rho }_2$ and $\rho$ are more smoothed near the contact discontinuity, as well as u and $\theta$ 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 $\beta =0.4$ and N fails. For four times smaller $\beta =0.1$, the results are not bad (see them for $N=4001$ in ([41], Example 3)), but now the graph of ${\rho }_1$ acquires an additional rather high, though very narrow, false step, and the graphs of ${\rho }_2$ and $\rho$ are slightly more smoothed, both near the point of contact discontinuity.

In case C, the results are specific, see Figure 6. Now, for $\beta =0.2$ and the same N, the graph of ${\rho }_1$ has a large “finger” at the contact discontinuity which height is about twice the value of ${\rho }_1$ to the left of that point. Surprisingly, the rest of graphs look rather accurate including even $\rho$ (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 $\rho$, p, u and $\theta$ for the mixture does not guarantee the same concerning both ${\rho }_1$ and ${\rho }_2$ for the components (the latter graphs are sometimes omitted from the numerical results).

 Example 4.

(Shock-bubble interaction problem) from ([55], Test 3.4). In this example, the structure of the initial parameters is more complicated than in previous examples:

$$({\rho }_1,{\rho }_2,u,1.4p,\gamma ){|}_{t=0}=\left\{\begin{array}{cc}(0,\phantom{.3764}0.1819,1.22,\phantom{00}1,\phantom{.5698}1.648),\hfill & |x+4|<\frac{1}{\sqrt{2}}\hfill \\ (1,\phantom{.3764}0,\phantom{.1819}1.22,\phantom{00}1,\phantom{.5698}1.4),\hfill & x<-3\hfill \\ (1.3764,0,\phantom{.1819}0.8864,1.5698,1.4),\hfill & \mathrm{otherwise}\hfill \end{array}\right.,$$

thus, in gas 1 with ${\gamma }_1=1.4$, there is “the bubble” of gas 2 with ${\gamma }_2=1.648$ moving to the right. Here, $t_{fin}=4$, and $M_{min}:={min}_i\frac{|u_i^{\overline{m}}|}{c_{si}^{\overline{m}}}\approx 0.32$ and $M_{max}\approx 1.22$, so the flow is transonic. In this problem, shock waves and contact discontinuities interact, that complicates computations. Let $a=0.5$ and $\beta =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 $\rho$ and p are very close to those only presented in [55]. The very small ledges in the graphs of ${\rho }_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\alpha }$ and $a_{Pr}$ as in previous examples. On the other hand, ${\rho }_1$, ${\rho }_2$, $\rho$ and $\theta$ are too smooth near the two mentioned points, especially for $N=1001$.

In case C, the graphs of ${\rho }_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 $\theta$ 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.