Mixtures: sequential stability of variational entropy solutions

The purpose of this work is to analyze the mathematical model governing motion of $n$-component, heat conducting reactive mixture of compressible gases. We prove sequential stability of weak variational entropy solutions when the state equation essentially depends on the species concentration and the species diffusion fluxes depend on gradients of partial pressures. Of crucial importance for our analysis is the fact that viscosity coefficients vanish on vacuum and the source terms enjoy the admissibility condition dictated by the second law of thermodynamics.

In (1), t denotes the time, t ∈ (0, T ) and the length of time interval T is assumed to be arbitrary large, but finite. The space domain Ω is a periodic box T 3 . The vectors belonging to the physical space R 3 as well as the tensors are printed in the boldface style. The species mass conservation equations can be equivalently written in terms of species mass fractions: where Y k , k ∈ {1, . . . , n} are defined by Y k = ̺ k ̺ and they satisfy: We remark that we will freely switch from one notation to the other using the species unknowns (̺, ̺ 1 , . . . , ̺ n ) or equivalently (̺, Y 1 , . . . , Y n ). The global-in-time existence of solutions for system (1), supplemented with physically relevant constitutive relations, was established by Giovangigli [24,chapter 9,Theorem 9.4.1]. He considered, for instance, generic matrix C kl relating the diffusion deriving forces d k to the species diffusion fluxes F k = − n l=1 C k,l d l + Soret effect, k = 1, ...n, which was singular CY = 0, Y = (Y 1 , . . . , Y n ) T , Y k = ̺ k ̺ , k = 1, . . . , n, and not symmetric in general. The result holds provided the initial data are sufficiently close to an equilibrium state and our main motivation is to extend it for the case of arbitrary large data. It should be however emphasized that, in comparison with [18], many simplifications are adopted in the present model. We concentrate only on diffusion effects due to the mole fraction and pressure gradients and restrict to particular form of the matrix C.
Due to mass conservation, the equations for species must sum to continuity equation which is hyperbolic. Moreover, the highly complex structure of diffusion fluxes causes that any standard approach for parabolic systems can not be applied unhindered. Nevertheless, it has been observed that even for strongly coupled systems of PDEs the concept of entropy can still be used in order to derive appropriate compactness results [7,17]. In our case, it is possible provided the matrix D kl = C kl ̺Y k , k, l = 1, . . . , n is symmetric and coercive on the hyperplanes which do not contain the vector Y . This assumption corresponds to the non-negativity of entropy production rate associated with diffusive process and has been postulated for example by Waldmann [24]. It turns out that such a description captures quite accurately lot of practical applications if only the species deriving forces d k are well-fitting.
The most exhaustively studied approximation is the Fick law, widely used especially in the simplified models: two species kinetics, 1D geometry, irreversible or isothermal reactions (cf. [6,8,9,25,27]). It states that the diffusion flux of a single species depends only on the gradient of its concentration. Hence, it does not take into account the cross-effects that play an important role in the multicomponent flows. As far as the latter are concerned, the issue of global-in-time existence of weak variational entropy solutions was investigated by Feireisl, Petzeltová and Trivisa [14]. They generalized the proof from [12] to the case of chemically reacting flows, when there is no interaction among the species diffusion fluxes and the pressure does not depend on the chemical composition of the mixture. The presence of the species concentration in the state equation would result in the entropy production rate which may fail to be non-negative. This in turn interferes with obtaining the fundamental a-priori estimates being the corner stone of the analysis presented in [14]. A similar problem was apparent in the approximation considered by Frehse, Goj and Málek [15,16] who proved existence and uniqueness of solutions to the steady Stokes-like system for a mixture of two (non-reactive) fluids. In their case, neglecting some nonlinear interaction terms in the source of momentum caused that the basic energy equality was no longer preserved.
Significantly less is known for models with general diffusion. The local-in-time well posedness of the Maxwell-Stefan equations with multicomponent diffusion in the isobaric, isothermal case is due to Bothe [1]. More recently, Mucha, Pokorný and Zatorska [22] proved that the system of n diffusionreaction equations with particular form of diffusion matrix C admits a global-in-time weak solution, provided the additional regularity of the total density is available. This extra information can be deduced from the new kind of estimate found for the Korteweg system by Bresch, Desjardins and Lin [5]. It also holds for the compressible Navier-Stokes equations [4] if the viscosity coefficients µ and ν are density-dependent and the stisfy On the other hand, such an assumption leads automatically to a problem with defining the velocity field on the vacuum zones. It was solved on the level of weak sequential stability of solutions by Mellet and Vasseur [20] who combined the previous results with logarithmic estimate for the velocity. Nevertheless, construction of sufficiently smooth approximate solutions in this framework remains an open problem. A way to overcome this obstacle in the case of heat-conducting fluids was proposed in a recent paper by Bresch and Desjardins [3]. They introduced a modification of the pressure close to the zero Kelvin isothermal curve by which it was possible to show that the sets of presence of vacuum have a zero Lebesgue measure. The same idea was then employed in [26,21] to prove the global-in-time existence of weak solutions for the isothermally reacting mixture of two gaseous components. In the present paper we indicate a possible way of generalization of these results to the case of n-component heat-conducting mixture.
Let us emphasize, that unlike to [3] we use the framework of weak variational solutions, that is to say, the total energy balance is replaced by the entropy inequality and the global total energy balance. These solutions were introduced by Feireisl [10] to study solvability of the NSF system, for which the weak solutions may dissipate more kinetic energy than, if there are any, the classical ones. However, it can be verified that this "missing energy" is equal to 0 provided the weak solution is regular enough and satisfies the total global energy balance. Moreover, in [13] the authors proved the weak-strong uniqueness of such solution, meaning that it coincides with the strong solution, emanating from the same initial data, as long as the latter exists.
The article is organized as follows. In the next section we specify the structural properties for the transport coefficients and postulate several simplifications. In Section 3, we define the notion of weak variational solutions and state the main result of the paper. The key a priori estimates are derived in Section 4 together with some further estimates and positivity of the absolute temperature. Finally, the last step of the proof of Theorem 5 -the limit passage -is performed in Section 5.

Constitutive relations and main hypothesis
Let us now supplement system (1) with a set of expressions determining the form of thermodynamical functions and transport fluxes in terms of macroscopic variable gradients in the spirit of [18], Chapter 2.
Thermal equation of state. We consider the pressure π = π(̺, ϑ, Y ), which can be decomposed into π = π m + π c , where the latter component depends solely on the density and it corresponds to the barotropic process of viscous gas. It is the only non-vanishing component of the pressure when temperature goes to zero, thus will be termed a "cold" pressure. We assume that π c is a continuous function satisfying the following growth conditions for positive constants c 1 , c 2 and γ − , γ + > 1. It was proposed in [3] to encompass plasticity and elasticity effects of solid materials, for which low densities may lead to negative pressures. By this modification the compactness of velocity at the last level of approximation can be obtained without requiring more a priori regularity than expected from the usual energy approach [20].
The first component π m = π m (̺, ϑ, Y ) is the classical molecular pressure of the mixture which is determined through the Boyle law as a sum of partial pressures p k : where e st k = const . is the formation energy of the k-th species, while c v is the constant-volume specific heat, which is assumed to be the same for each of the species. The "cold" components of the internal energy e c = e c (̺) and pressure π c are related through the following equation of state: The last relation is a consequence of the second law of thermodynamics which postulates the existence of a state function called the entropy.
The entropy equation. The entropy of a thermodynamical system is defined (up to an additive constant) by the differentials of energy, total density, and species mass fractions via the Gibbs relation: where D denotes the total derivative with respect to the state variables {̺, ϑ, Y }; whereas g k are the Gibbs functions Here, h k = h k (ϑ) denotes the specific enthalpy and s k = s k (ϑ, ̺ k ) is the specific entropy of the k-th species h k (ϑ) = e st k + c pk ϑ, where s st k = const. denotes the formation entropy of k-th species. The constant-volume (c vk ) and constant-pressure (c pk ) specific heats for the k-th species are constants related by In accordance with (6), (7) and (8), the specific entropy of the mixture can be expressed as a weighted sum of the species specific entropies and is governed by the following equation where σ is the entropy production rate By virtue of the second law of thermodynamics, the entropy production rate must be non-negative for any admissible process.
The stress tensor. The viscous part of the stress tensor obeys the Newton rheological law, namely: where D(u) = 1 2 ∇u + (∇u) T and µ = µ(̺), ν = ν(̺) are C 1 (0, ∞) functions related by which is strictly a mathematical constraint allowing to obtain better regularity of ̺, [2]. In addition, they enjoy the following bounds for some positive constant µ ′ ,

Remark 1
The assumption µ ′ ≤ µ ′ (̺) is not optimal, but it makes the proof much easier, however the bound from above is essential in order to get integrability of several important quantities. For further discussion on this topic we refer to [20].
The species diffusion fluxes. Following [18], Chapter 2, Section 2.5.1, we postulate that the diffusion flux of the k-th species is given by where C 0 , C kl are multicomponent flux diffusion coefficients and d k is the species k diffusion force specified, in the absence of external forces, by the following relation We restrict to an exact form of the flux diffusion matrix C compatible with the set of mathematical assumptions postulated in [18], Chapters 4, 7 (see also references therein). The prototype example is the same as studied in [22], namely where Concerning the diffusion coefficient C 0 , we assume that it is a continuously differentiable function of ϑ and ̺ and that there exist positive constants C 0 , C 0 such that Remark 2 One of the main consequences of (20) is that Remark 3 Note that (20) also implies that the vector of species diffusion forces d = (d 1 , . . . , d n ) T is an eigenvector of the matrix C corresponding to the eigenvalue 1 and we recast that The heat flux. It is a sum of two components where the first term describes transfer of energy due to the species molecular diffusion, whereas the second term represents the Fourier law with the thermal conductivity coefficient In the above formulas κ 0 , κ 0 , α are positive constants and α ≥ 2.
The species production rates. We assume that ω k are continuous functions of Y bounded from below and above by the positive constants ω and ω we also suppose that Another restriction is dictated by the second law of thermodynamics, which asserts in particular that ω k must enjoy the following condition Initial data. The choice of quantities describing the initial state of system (1) is dictated by the weak formulation of the problem, which is specified in Definition 4 below. We take In addition, we assume that ̺ 0 is a nonnegative measurable function such that and the initial densities of species satisfy Further, we call the initial temperature ϑ 0 a measurable function such that and the following compatibility condition is satisfied Finally, we require that the initial distribution of the momentum is such that and the global total energy at the initial time is bounded

Weak formulation and main result
In what follows we define a notion of weak variational entropy solutions to system (1) and then we formulate our main result.

The balance of momentum
holds for any smooth test function φ(t, x) such that φ(T, ·) = 0.

The entropy equation
is satisfied for any smooth function φ(t, x), such that φ ≥ 0 and φ(T,
In this definition the usual weak formulation of the energy equation is replaced by the weak formulation of the entropy inequality and the global total energy balance (36)+(37). Note however that the entropy production rate has now only a meaning of non-negative measure which is bounded from below by the classical value of σ. Nevertheless, a simple calculation employing the Gibbs formula (8) shows that whenever the solution specified above is sufficiently regular, both formulations are equivalent. In particular the entropy inequality (36) changes into equality, see f.i. Section 2.5 in [11].
We are now in a position to formulate the main result of this paper.
. . , n is a sequence of smooth solutions to (1) with the initial data satisfying (30-33), moreover Then, up to a subsequence, {̺ N , u N , ϑ N , ̺ k,N } converges to the weak solution of problem (1) in the sense of Definition 4.
The proof of this theorem can be divided into two main steps. The first one is dedicated to derivation of the energy-entropy estimates which are obtained under assumption that all the quantities are sufficiently smooth. The next step is the limit passage, it consists of various integrability lemmas combined with condensed compactness arguments.

A priori estimates
In this section we present the a priori estimates for which is a sequence of smooth functions solving (1). As mentioned above, assuming smoothness of solutions, we expect that all the natural features of the system can be recovered. The following estimates are valid for each N = 1, 2, . . . but we skip the subindex when it does not lead to any confusion.

Estimates based on the maximum principle
To begin, observe that the total mass of the fluid is a constant of motion, meaning Moreover if solution is sufficiently smooth, a classical maximum principle can be applied to the continuity equation in order to show that ̺ N (t, x) > c(N ) ≥ 0, more precisely in particular ̺ > 0.
Next, by a very similar reasoning we can prove non-negativity of ϑ on [0, T ] × Ω.
Lemma 6 Assume that ϑ = ϑ N is a smooth solution of (1), then Proof. Any solution to (1) which is sufficiently smooth is automatically a classical solution of the system ∂ t ̺ + div(̺u) = 0, where we replaced the total energy balance by the internal energy balance. Let us hence reformulate it in order to obtain the equation describing the temperature. Subtracting from the third equation of (44) the component corresponding to the formation energy we obtain where we used the species mass balance equations. Next renormalizing the continuity equation and employing the relation between e c and π c (7) we get the temperature equation where, in accordance with hypotheses (17) the second term on the r.h.s. is nonnegative. Consequently, (43) is obtained by application of the maximum principle to the above equation, recalling that ess inf x∈Ω ϑ 0 N (x) > 0. ✷ The analogous result for partial masses is stated in the following lemma.
Lemma 7 For any smooth solution of (1) we have Proof. We integrate each of equations of system (1) over the set {̺ k < 0}. Assuming that the boundary i.e. {̺ k = 0} is a regular submanifold we obtain thus |{̺ k < 0}| = 0, for every k = 1, . . . , n. When {̺ k = 0} is not a regular submanifold we construct a sequence {ε l } ∞ l=1 such that ε l → 0 + and {̺ k = ε l } is a regular submanifold and pass with ε l to zero.
The proof of (47) follows by subtracting the sum of species mass balances equations from the continuity equation. The smooth solution of the resulting system must be, due to the initial conditions (31), equal to 0 on [0, T ] × Ω. ✷ As a corollary from this lemma we recover relation (2), moreover we have the following estimate

The energy-entropy estimates
The purpose of this subsection is to derive a priori estimates resulting from the energy and entropy balance equations. The difference comparing to estimates obtained in the previous subsection is that now we look for bounds which are uniform with respect to N . We start with the following Lemma.
Lemma 8 Every smooth solution of (1) satisfies Proof. Integrate the third equation of (1) with respect to the space variable and employ the periodic boundary conditions. ✷ Assuming integrability of the initial conditions (33) the assertion of the above lemma entails several a priori estimates: It is well known that these natural bounds are not sufficient to prove the weak sequential stability of solutions, not even for the barotropic flow. However, taking into account the form of viscosity coefficients (16), (17), further estimates can be delivered.

Lemma 9
For any smooth solution of (1) we have Proof. Multiply the momentum equation by u and integrate over Ω. ✷ The above lemma can not be used to deduce the uniform bounds for the symmetric part of the gradient of u immediately. The reason for that is lack of sufficient information for ϑ, so far we only know (50). However, taking into account the form of viscosity coefficients (16), further estimates can be delivered. The following lemma is a modification of the result proved by Bresch and Desjardins for heat conducting fluids [3].
Lemma 10 Any smooth solution of (1) satisfies the following identity Proof. The rough idea of the proof is the following. The terms from the l.h.s. of this equality can be evaluated by multiplication of the momentum equation by ∇φ(̺) and the continuity equation by |∇φ(̺)| 2 . Then one has to combine these equivalences with the balance of kinetic energy (51) and include (16) to see that some unpleasant terms cancel. For more details we refer to [26] or to the original work of Bresch and Desjardins [3] . ✷ To control the r.h.s. of (51) and (52) one needs i.a. to estimate the gradient of ϑ. To this purpose we take advantage of the entropy balance (13), we have the following inequality Lemma 11 For any smooth solution of (1) we have Proof. Combining the third equation of (1) with the Gibbs relation (8) we derive the entropy equation Integrating it over space and time we obtain where the l.h.s. can be transformed using (9) and (24) into The first two terms on the l.h.s. of (55) have a good sign, the same holds for the last one due to (28). Non-negativity of the third one follows from (22) and (23) Thus, it remains to control the positive part of ̺s(T ) and the negative part of (̺s) 0 . From definition of the entropy (10) we get therefore The two first terms from the r.h.s. are bounded due to (50), whereas to estimate the positive part of the last one we essentially use the assumption that Ω is a bounded domain. Thus, the positive part of −x log x is bounded by a constant, and thus integrable over Ω. ✷ In the rest of Section 4 we show how to use Lemmas 10 and 11 in order to derive uniform estimates for the sequence of smooth solutions {̺ N , u N , ϑ N , ̺ k,N } ∞ N =1 to system (1).

Estimates of the temperature.
One of the main consequences of (53) is that for κ(̺, ϑ) satisfying (25) we have the following a priori estimates for the temperature where s ∈ [0, α 2 ] and α ≥ 2. To control the full norm of ϑ s in L 2 (0, T ; W 1,2 (Ω)) we will apply the following version on the Korn-Poincaré inequality (see f.i. Theorem 10.17 in [23]):

Theorem 12
Let Ω ⊂ R 3 be a bounded Lipschitz domain. Assume that r is a non-negative function such that 0 < M 0 ≤ Ω r dx, Ω r γ dx ≤ K, f or a certain γ > 1. Then for any ξ ∈ W 1,p (Ω).

Estimates following from the Bresch-Desjardin equality
The aim of this subsection is to derive estimates following from (51) and (52). Summing these two expressions we obtain d dt Ω We first need to justify that the terms from the r.h.s. are bounded or have a negative sign so that they can be moved to the l.h.s. The main problem is to control the contribution from the molecular pressure. To this purpose we will employ estimate (53) coupled with the properties of the diffusion matrix C. Denoting where we obtain, for every k-th coordinate k ∈ {1, . . . , n} and every i-th space coordinate i ∈ {1, 2, 3}, the following decomposition Next, multiplying the above expression by m k and summing over k ∈ {1, . . . , n} one gets Returning (62) we can express the full gradients of partial pressures in terms of gradients of temperature, density and the gradient of "known" part of the pressure As was announced, we will use the above expression in order to control the molecular part of the pressure from the r.h.s. of (59).
Estimate of ∇π(̺, ϑ, Y ) · ∇φ. Using definition of φ (see Lemma 10) and (4) we obtain The first term is non-negative due to (5), so it can be considered on the l.h.s. of (59) and we only need to estimate the second one. Since Note that I 2 is non-negative, so we can put it to the l.h.s. of (59). Next, I 1 and I 4 can be estimated in a similar way, we have so for ε sufficiently small, the first term can be controlled by I 2 thanks to (17). Concerning the second integral, from (53) we have Using (23), the integral may be transformed as follows thus, due to (60) and (21), the integral over time of the r.h.s. of (66) is bounded. For I 3 we verify that and the first term is bounded in view of (57) whereas boundedness of the second one follows from the Gronwall inequality applied to (59). Indeed, note that, due to (25), ̺ϑ 2 κ(̺,ϑ) is bounded by some positive constant.
Estimate of π(̺, ϑ, Y ) div u. By virtue of (4) and (7) and the continuity equation Furthermore, by the Cauchy inequality On account of (58), ϑ ∈ L 2 (0, T ; L 6 (Ω)). Moreover, the Sobolev imbedding theorem implies that for 1 ≤ p ≤ 6, hence the Gronwall inequality applied to (59) implies boundedness of (69), whence the term ε µ(̺) div u 2 L 2 (Ω) is then absorbed by Ω S : ∇u dx from the l.h.s. of (59). Resuming, we have proved the following inequality: Uniform estimates. Taking into account all the above considerations, we can complement the so-far obtained estimates as follows Concerning the velocity vector field, in addition to (50) we have

Estimates of species densities
Finally, we can take advantage of the entropy estimate (53) which together with (72) may be used to deduce boudedness of gradients of all species densities.

Lemma 13
For any smooth solution of (1) we have Proof. First, using Remark 3 we may write which is bounded in L 1 ((0, T ) × Ω) on account of (67). Clearly, The r.h.s. of above can be, due to (63), estimated as follows which is bounded thanks to (57), (68) and (71). In consequence, (75) is bounded. Recalling assumptions imposed on C 0 (21) and the form of molecular pressure π m , we deduce that and the r.h.s. is bounded, again by (48) and (57). ✷

Additional estimates.
In this subsection we present several additional estimates based on imbeddings of Sobolev spaces and the simple interpolation inequalities.
Estimate of the velocity vector field. We use the Hölder inequality to write where p = 2γ − γ − +1 , q = 6γ − 3γ − +1 . Therefore, Theorem 12 together with the Sobolev imbedding imply Next, by the similar argument with p ′ = 2γ − , q ′ = 6γ − 3γ − +1 . By a simple interpolation between (80) and (81), we obtain and since γ − > 1, we see in particular that u ∈ L Proof. The above statement is a consequence of the following estimate which can be obtained, again by application of Theorem 12 with ξ = log ϑ N and r = ̺ N . It remains to check that we control the L 1 (Ω) norm of ̺| log ϑ|. By virtue of (36) we have thus substituting the form of ̺s from (56) we obtain and the r.h.s. is bounded on account of (48), (78) and the initial condition. On the other hand, the positive part of the integrant ̺ N log ϑ N is bounded from above by ̺ N ϑ N which belongs to L ∞ (0, T ; L 1 (Ω)) due to (50), so we end up with ess sup which was the missing information in order to apply Theorem 12. This completes the proof of (84). ✷

Passage to the limit
In this section we justify that it is possible to perform the limit passage in the weak formulation of system (1). We remark that we focus only on the new features of the system, i.e. the molecular pressure and multicomponent diffusion, leaving the rest of limit passages to be performed analogously as in [26].

5.1
Strong convergence of the density and passage to the limit in the continuity equation.
The strong convergence of a sequence {̺ N } ∞ N =1 is guaranteed by the following lemma On account of that, it is easy to let N → ∞ in the continuity equation to obtain (34).

Strong convergence of the species densities.
Analogously we show the strong convergence of species densities. We have Lemma 17 Up to a subsequence the partial densities ̺ k,N , k = 1, . . . , n converge strongly in L p (0, T ; L q (Ω)), 1 ≤ p < ∞, 1 ≤ q < 3 to ̺ k . In particular Moreover ̺ k,N → ̺ k in C([0, T ]; L 3 weak (Ω)). Proof. The estimate (74) together with (48) and (78) give the bound for the space gradients of ̺ k,N , Moreover, directly from the equation of species mass conservation we obtain Indeed, the most restrictive term is the diffusion flux, which can be rewritten as Due to (48) we have that which is bounded in L 2 (0, T ; L 3 2 ) on account of(57) and (78). Similarly, thus, according to (74) it remains to control the norm of (ϑ N + 1)̺ N in L p ((0, T ) × Ω) for some p > 2. By (58) and (78) we deduce that ̺ N ϑ N ∈ L α (0, T ; L 3α α+1 (Ω)), so (89) is verified. In this manner we actually proved that the sequence of functions is uniformly bounded and equicontinuous in C([0, T ]), hence, the Arzelá-Ascoli theorem yields Since ̺ k,N is bounded in L ∞ (0, T ; L 3 (Ω)) and due to density argument, this convergence extends to each φ ∈ L 3 2 (Ω). Finally, the Aubin-Lions argument implies the strong convergence of the sequence ̺ k,N to ̺ k in L p (0, T ; L q (Ω)) for p = 2, q < 3, but due to (48) and (78) it can be extended to the case p < ∞. ✷

Strong convergence of the temperature.
From estimate (58) we deduce existence of a subsequence such that ϑ N → ϑ weakly in L 2 (0, T ; W 1,2 (Ω)), however, time-compactness cannot be proved directly from the internal energy equation (45). The reason for this is lack of control over a part of the heat flux proportional to ̺ϑ α ∇ϑ. This obstacle can be overcome by deducing analogous information from the entropy equation (54).
We will first show that all of the terms appearing in the entropy balance (54) are nonnegative or belong to W −1,p ((0, T ) × Ω), for some p > 1. Indeed, first recall that due to (56) whence due to (50), (78) and (84) we deduce that The entropy flux is due to (9) and (24) The first part can be estimated as follows where the most restrictive term can be controlled as follows N |, which is bounded on account of (57) provided ̺ N ϑ α N is bounded in L p ((0, T ) × Ω) for p > 1, uniformly with respect to N . Note that for 0 ≤ β ≤ 1 we have is bounded in L p (0, T ; L q (Ω)), for p and q satisfying 1 p = α−β α , 1 q = β + 1−β 3 + α−β 3α . In particular p, q > 1 provided 0 < β < α 2α−1 .
The remaining part of the entropy flux is equal to where the middle term vanishes due to (22). The worst term to estimate is thus the last one, we rewrite it using (23) in the following way . . , n. Both parts have the same structure, so we focus only on the first one, we have Using (57), (58), (74) and (78) we finally arrive at is bounded in L p ((0, T ) × (Ω)), for 1 < p < 4 3 . (95) We are now ready to proceed with the proof of strong convergence of the temperature. To this end we will need the following variant of the Aubin-Lions Lemma.

Limit in the momentum equation, the species mass balance equations and the global total energy balance
Having proven pointwise convergence of sequences we are ready to perform the limit passage in all the nonlinear terms appearing in the momentum equation, the species mass balance equations and the total global energy balance.

Lemma 19
Let p > 1, then up to a subsequence we have Proof. We already know that ̺ N converges to ̺ a.e. on (0, T ) × Ω. Moreover, due to (80), up to extracting a subsequence, u N converges weakly to u in L p (0, T ; L q (Ω)) for p > 1, q > 3. Therefore, the uniform boundedness of the sequence ̺ N u N in L ∞ (0, T ; L 3 2 (Ω)) implies that Now, we are aimed at improving the time compactness of this sequence. Using the momentum equation, we show that the sequence of functions is uniformly bounded and equicontinuous in C([0, T ]), where φ ∈ C ∞ c (Ω). But since the smooth functions are dense in L 3 (Ω), applying the Arzelà-Ascoli theorem, we show (102).
(iv) Convergence of the diffusion terms. In the proof of Lemma (17) it was shown in particular that {F k,N } ∞ N =1 is bounded in L 4 3 ((0, T ) × (Ω)). By the weak convergence of ∇̺ k,N , ∇ϑ N to ∇̺ k , ∇ϑ, respectively, deduced from (88) and (91) together with (87), (100) and (77) we check that it is possible to let N → ∞ in all terms of (90). In other words, we have The convergence results established above are sufficient to perform the limit passage in the momentum, the total global energy balance and the species mass balance equations and to validate, that the limit quantities satisfy the weak formulation (35),(37) and (38).

Limit in the entropy inequality
In view of (92-95) and the remarks from the previous subsection, it is easy to pass to the limit N → ∞ in all terms appearing in (36), except the entropy production rate σ. However, in accordance with (53) we still have that is bounded in L 2 ((0, T ) × Ω).
Moreover, by virtue of (79), (86) and (100) we deduce weakly in L 2 ((0, T ) × Ω). Evidently, we may treat all the remaining terms in the similar way using the fact that they are linear with respect to the weakly convergent sequences of gradients of u N , ϑ N and ̺ k,N . Thus, preserving the sign of the entropy inequality (36) in the limit N → ∞ follows by the lower semicontinuity of convex superposition of operators.
In order to show (106) we thus need to justify that σ is absolutely continuous with respect to Lebesgue measure on [0, τ ] × Ω. To this end we use in (107) a test function φ = ϑ 0 , we get