Global existence of weak solutions and weak–strong uniqueness for nonisothermal Maxwell–Stefan systems

The dynamics of multicomponent gas mixtures with vanishing barycentric velocity is described by Maxwell–Stefan equations with mass diffusion and heat conduction. The equations consist of the mass and energy balances, coupled to an algebraic system that relates the partial velocities and driving forces. The global existence of weak solutions to this system in a bounded domain with no-flux boundary conditions is proved by using the boundedness-by-entropy method. A priori estimates are obtained from the entropy inequality which originates from the consistent thermodynamic modelling. Furthermore, a conditional weak–strong uniqueness property is shown by using the relative entropy method.


Introduction
The dynamics of multicomponent gaseous mixtures with vanishing barycentric velocity and constant temperature can be described by the Maxwell-Stefan equations [23,26].The existence of local-in-time smooth and global-in-time weak solutions to these systems has been proved in [2,13,16,21].The analysis of nonisothermal gas mixtures is, however, incomplete.The existence of local-in-time solutions was shown in [19], while [15] investigated a special nonisothermal case.In this paper, we prove the existence of global-in-time weak solutions and the weak-strong uniqueness property for a rather general nonisothermal Maxwell-Stefan system.The novelty of our approach is the consistent thermodynamic modeling.
1.1.Model equations.The evolution of the mass densities ρ i (x, t) of the ith gas component and the temperature θ(x, t) of the mixture is described by the mass and energy balances (ρ j e j + p j )u j in Ω, t > 0, (2) where Ω ⊂ R 3 is a bounded Lipschitz domain, J i and J e are the diffusion and energy fluxes, respectively, u i are the diffusional velocities, ρ = n i=1 ρ i is the total mass density, p i the partial pressure with the total pressure p = n i=1 p i , ρ i e i the partial internal energy ρ i e i with the total energy ρe = n i=1 ρ i e i , and κ(θ) is the heat conductivity.Equations ( 1)-( 2) are supplemented with the boundary and initial conditions where ν is the exterior unit normal vector to ∂Ω, θ 0 > 0 is the given background temperature, and λ > 0 is a relaxation constant.The boundary conditions mean that the gas components cannot leave the domain, while heat exchange through the boundary is possible and proportional to the difference between the gas and background temperatures.To close the model, we need to determine u i , ρ i e i , and p i .
The velocities u i are computed from the constrained algebraic Maxwell-Stefan system where the constant coefficients b ij = b ji > 0 model the interaction between the ith and jth components.The driving force d i is given by ( 6) where µ i is the chemical potential.The constraint (7) ∇p = 0 in Ω, t > 0, is needed in order for our system to be thermodynamically consistent.We refer to Section 2 for details.
The internal energies ρ i e i and chemical potentials µ i are determined from the Helmholtz free energy (see (16)), and the pressure is computed from the Gibbs-Duhem relation.As shown in Section 2, these quantities are explicitly given by ( 8) where ρ i η i is the entropy density of the ith component and c w > 0 is the heat capacity.Then the driving force d i and energy flux J e simplify to (9) The matrix M associated to the algebraic system (5) is singular (since n i=1 d i = 0) and thus not positive definite.However, we recall in Section 3.1 that it is positive definite on the subspace L = {y = (y 1 , . . ., y n ) ∈ R n : √ ρ • y = 0} (here, √ ρ is the vector with components √ ρ i ).Therefore, the Bott-Duffin inverse of M, denoted by M BD = M BD (ρ), exists and is symmetric and positive definite on L.Moreover, we show in Section 3.3 below that the fluxes can be expressed as a linear combination of the entropy variables (or thermo-chemical potentials) µ/θ = (µ 1 /θ, . . ., µ n /θ) and −1/θ, (10) , where Q(ρ, θ) = A B B T a , and A = (A ij ) ∈ R n×n , B = (B i ) ∈ R n , a > 0 are given by ( 11) A ij m i m j .
Here, variables in bold font are n-dimensional vectors.The Onsager matrix Q turns out to be positive semidefinite (see (33)), which reveals the parabolic structure of equations ( 1)-(2).
1.2.State of the art.The isothermal Maxwell-Stefan equations can be derived from the multispecies Boltzmann equations in the diffusive approximation [6].The high-friction limit in Euler (-Korteweg) equations reveals a formal gradient-flow form of the Maxwell-Stefan equations [17], leading to Fick-Onsager diffusion fluxes instead of (5).In fact, it is shown in [5] that the Fick-Onsager and generalized Maxwell-Stefan approaches are equivalent.A formal Chapman-Enskog expansion of the stationary nonisothermal model was given in [27].Another nonisothermal Maxwell-Stefan system was derived in [1], but with a different energy flux than ours.Maxwell-Stefan systems with nonvanishing barycentric velocities can be formulated in the framework of hyperbolic-parabolic systems, which allows one to perform a local-in-time existence analysis [13].Global-in-time regular solutions around the constant equilibrium state were found to exist in [14].An existence analysis for Maxwell-Stefan systems coupled to the Navier-Stokes equations for the barycentric velocity can be found in [8] for the incompressible case and in [4] for the compressible situation.For steady-state problems, we refer to, e.g., [7,24].
When the barycentric velocity vanishes, the (isothermal) Maxwell-Stefan equations can be solved by generalized parabolic theory.The existence of local-in-time classical solutions was proved in [2], while the existence of global-in-time weak solutions with general initial data was shown in [21].Concerning the nonisothermal equations, we refer to [15], where an existence analysis for global-in-time weak solutions was presented.However, this model has some modeling deficiencies explained below.Therefore, our first aim is to prove the global existence for a thermodynamically consistent nonisothermal model.
The uniqueness of strong solutions to the isothermal Maxwell-Stefan equations was shown in [2,16,19], but the uniqueness of weak solutions for general coefficients b ij is still unsolved.A very special case (the coefficients b ij have two degrees of freedom only) was investigated in [9].It was shown in [18] that strong solutions are unique in the class of weak solutions, which is known as the weak-strong uniqueness property.Our second aim is to prove this property for the nonisothermal case.
Let us detail the main differences of our work compared to [15]: (i) The most important difference is the lack of validity of the Onsager reciprocity relations in the model of [15].The relations imply the symmetry of the coefficients of the Onsager matrix; see (10).The choice in [15] leads to a cancelation in the entropy inequality, thus simplifying the estimation.Our results do not rely on this simplification; see Remark 6 for further details.(ii) The constraint (7) on the pressure is not taken into account in [15].This condition is not necessary mathematically, but its lack creates an inconsistency with the assumption of vanishing barycentric velocity.Indeed, a difference in pressure induces a force difference, which can result in an acceleration according to Newton's second law, if there is no additional force to balance it.(iii) According to Onsager's reciprocity relations, the Onsager matrix Q in (10) has to be positive semidefinite.We show that Q is in fact positive definite on the subspace L = {y ∈ R n : y • √ ρ = 0}.In [15], is is assumed that this subspace equals {y ∈ R n : y • 1 = 0}.This is not consistent with the thermodynamic modeling.(iv) We consider different molar masses m i , while they are assumed to be the same in [15].When we assume equal molar masses, the cross-terms cancel, and we end up with the simple heat flux J e = −κ∇θ (see (9) and the constraint in ( 5)), thus decoupling the equations.
1.3.Main results.We impose the following assumptions: (A1) Domain: Ω ⊂ R 3 is a bounded domain with Lipschitz boundary, and T > 0. We set Ω T = Ω × (0, T ) and The lower bound for the total mass density ρ is needed to derive uniform estimates for the temperature.The proof of Lemma 10 in [18] shows that M BD ij (ρ) is bounded for all ρ ∈ R n + .The growth condition for the heat conductivity is used to derive higher integrability bounds for the temperature, which are needed to derive a uniform estimate for the discrete time derivative of the temperature.We may also assume reaction terms R i in (1) with the properties that the total reaction rate n i=1 R i vanishes and the vector of reaction rates R i is derived from a convex, nonnegative potential [11,Section 2.2].
The first main result is the existence of solutions.
Our second main result concerns the weak-strong uniqueness property.
By a strong solution, we understand a solution that has sufficient regularity to satisfy the entropy equality stated in Lemma 14; see Section 5. Observe that we require the boundedness of the temperature θ, which is not proved in Theorem 1.The proof of Theorem 2 is based on the relative entropy, defined by where E = c w ρθ and Ē = c w ρ θ are the internal energy densities.The idea is to compute the time derivative: where c > 0 is some constant and C > 0 depends on the L ∞ (Ω T ) norms of θ, ūi , and ∇ log θ, i = 1, . . ., n.The difficulty is to estimate the expressions arising from the time derivative of the relative entropy in such a way that only ūi and θ need to be bounded.Thanks to the positive lower bound for θ, we can bound the right-hand side in terms of the relative entropy, Then Gronwall's lemma shows that H((ρ, θ)(t)|( ρ, θ)(t)) = 0 for t > 0 and hence (ρ, θ)(t) = ( ρ, θ)(t).Compared to [18], we include the temperature terms and combine them with the entropy variables w i in such a way that the positive semidefiniteness of M BD can be exploited.
The paper is organized as follows.We detail the thermodynamic modeling of equations ( 1)-( 8) in Section 2. The inversion of the Maxwell-Stefan system (5), the definition of the (relative) entropy variables, and the formulations of the fluxes in terms of the relative entropy variables, as well as the corresponding weak formulation is presented in Section 3. Section 4 is concerned with the proof of Theorem 1, and Theorem 2 is proved in Section 5.

Modeling
We consider the following system of equations modeling the dynamics of a nonisothermal gas mixture of n components with mass diffusion and heat conduction: Besides of the variables introduced in the introduction, v denotes the barycentric velocity of the mixture.The quantities ρ i b i are the body forces, where ρb = n i=1 ρ i b i is the total force exerted on the mixture, and ρr is the total heat supply due to radiation.The diffusional velocities u i , the partial internal energy densities ρ i e i , and the partial pressures p i are determined from the free energy; see below.
Equations ( 13)-( 15) correspond to a so-called class-I model.They can be derived either via an entropy invariant model reduction [3] or in the high-friction limit [12] from a class-II model, in which each component has its own velocity v i .Equations ( 13) are the partial mass balances, ( 14) is the momentum balance, and (15) the energy balance.As proved in [12], system ( 13)-( 15) and ( 5) fits into the general theory of hyperbolic-parabolic compositetype systems introduced in [22] and further explored in [25].
As mentioned in the introduction, system (1)-( 2) and ( 7) is supplemented by the constrained Maxwell-Stefan system (5) for the velocities u i .These equations can be derived from a class-II model in the diffusion approximation [3, Section 14, (210)] or in the highfriction limit [12, Section 2, (2.50)] with the driving forces where µ i is the chemical potential of the ith component.Since the pressure is uniform in space, ∇p = 0, and we have neglected external forces, the driving force becomes (6).Then equations ( 1)-( 2) and ( 7) are obtained by setting v = 0 and r = b i = 0.The internal energy densities ρ i e i , partial pressures p i , and the chemical potential µ i are determined from the Helmholtz free energy.We assume that the gas is a simple mixture, which implies that these quantities can be calculated from the partial free energy densities ψ i (ρ i , θ), i = 1, . . ., n.We have where ρ i η i is the entropy density of the ith component and the equation for p i is called the Gibbs-Duhem relation.Defining the partial Helmholtz free energy as ( 16) the thermodynamic quantities are given by (8).Moreover, the driving force d i and enthalpy h i := ρ i e i + p i read as ( 17) This corresponds to equations (9).

Preparations
3.1.Inversion of the Maxwell-Stefan system.We discuss the inversion of the Maxwell-Stefan system (5) following [12] and [18, Section 2].We write (5) equivalently as where the matrix We wish to invert Mv = w, where shows that the kernel of M consists of span{ √ ρ}.Thus, we can invert M only on the subspace We define the projections P L on L and P L ⊥ on L ⊥ by where δ ij is the Kronecker symbol.The matrix we can define the Bott-Duffin inverse of M with respect to L as Hence, we can invert ( 18) by ( 21) The matrix M BD = M BD (ρ) is symmetric and positive definite on L [18, Lemma 4], 3.2.Entropy variables.The mathematical analysis becomes easier when formulating the system in terms of the so-called entropy variables.To this end, we introduce the mathematical entropy density which is the negative of the physical (total) entropy density (8).Summing the mass balances (1) over i = 1, . . ., n and using the constraint n i=1 ρ i u i = 0 from (5), we obtain ∂ t ρ = 0. Thus, the total density is determined by the initial total density, ρ(x, t) = n i=1 ρ 0 i (x) for x ∈ Ω, and is independent of time.This suggests to compute only the first n − 1 mass densities, since the last one can be determined by ρ n = ρ − n−1 i=1 ρ i .Then we interpret the entropy density h as a function of (ρ ′ , θ) := (ρ 1 , . . ., ρ n−1 , θ): The Hessian matrix is positive definite, showing that the entropy is convex.According to thermodynamics [3], the entropy variables equal (µ 1 /θ, . . ., µ n /θ, −1/θ).We set (24) Since the nth partial density is determined by the densities ρ 1 , . . ., ρ n−1 , we prefer to work with the relative entropy variables (25) Setting additionally w = log θ, our new set of variables is (w 1 , . . ., w n−1 , w).The following lemma states that the mapping (ρ 1 , . . ., ρ n , θ) → (w 1 , . . ., w n−1 , w) is invertible.
The previous proof shows that we can formulate the diffusion fluxes in different ways.
Corollary 5.It holds for i = 1, . . ., n that We claim that the Onsager matrix Q ∈ R (n+1)×(n+1) in ( 10) is positive semidefinite.Let a = θ(κ + n i,j=1 A ij /(m i m j )).We compute for ξ ∈ R n+1 : where the nonnegativity follows from the positive semidefiniteness (22) of M BD .This reveals the parabolicity of our system in terms of the entropy variables.
3.4.Weak formulation.The previous subsection shows that we can write our system as the mass and energy balances (1)-( 2) with the fluxes ( 26)- (27).The weak formulation in the relative entropy variables (25) reads as According to (8), the energy is given by E = c w ρθ.Moreover, ρ i , A ij , B i , and E are interpreted as functions of (w 1 , . . ., w n−1 , w).

Proof of Theorem 1
The proof follows the lines of [15, Section 3], which is based on the boundedness-byentropy method [20], but some details are different.We approximate equations (34)-( 35) by replacing the time derivative by the implicit Euler scheme and adding a higher-order regularization in w i .The existence of solutions to the approximate system is shown by means of the Leray-Schauder fixed-point theorem, where the compactness of the fixedpoint operator is obtained by the approximate entropy inequality.This inequality yields estimates uniform in the regularization parameters, allowing for the de-regularization limit via the Aubin-Lions compactness lemma.
Let ε ∈ (0, 1), N ∈ N, and τ = T /N.We set w 0 = log θ 0 and w = (w 1 , . . ., w n−1 , w).Let w = ( w1 , . . ., wn−1 , w) ∈ L ∞ (Ω; R n ) be given.We define for test functions φ i ∈ H 2 (Ω), i = 0, . . ., n − 1, the approximate scheme where D 2 w i is the Hesse matrix of w i , the double point ":" denotes the Frobenius matrix product, we recall that E(w) = c w ρθ, and A ij and B i are interpreted as functions of w.The higher-order regularization yields solutions w i , w ∈ H 2 (Ω), and the W 1,4 (Ω) regularization allows us to estimate the higher-order terms when using the test function e −w 0 − e −w (see the estimate of I 11 below).The lower-order regularization (e w 0 − e w )(w − w 0 ) provides an ε-dependent L 2 (Ω) bound for w.

Solution of the linearized approximate problem. Let w
where where we abbreviated ρ * i = ρ i (w * ), ρi = ρ i ( w), E * = c w ρe w * , and Ē = c w ρe w.The bilinear form a is clearly coercive on H 2 (Ω; R n ), and both a and F are continuous on this space.By the Lax-Milgram lemma, there exists a unique solution w ∈ H 2 (Ω; R n ) to (38).4.2.Solution of the approximate problem.The solution w ∈ H 2 (Ω; R n ) to (38) defines the fixed-point operator S : The operator is continuous, compact (because of the compact embedding H 2 (Ω; R n ) ֒→ W 1,4 (Ω; R n )), and it satisfies S(w * , 0) = 0 for all w * ∈ W 1,4 (Ω; R n ).It remains to find a uniform bound for all fixed points of S(•, σ).Let w ∈ H 2 (Ω; R n ) be such a fixed point.Then w solves (38) with w * = w.We choose the test functions φ i = w i for i = 1, . . ., n − 1 and φ 0 = e −w 0 − e −w in (38): We estimate the terms I 1 , . . ., I 11 step by step.First, by the convexity of the entropy and arguing similarly as in [15, Section 3, Step 2], This expression is nonnegative because of the positive semidefiniteness of see (22).Furthermore, since sinh(z)/z ≥ 1 for z ∈ R, z = 0, I 9 = σλ ∂Ω e −w−w 0 (e w − e w 0 ) 2 dx ≥ 0, We claim that there exists m = m(w 0 , σ) > 0 such that for all w ∈ R, where w 0 ∈ R and σ ∈ (0, 1] are given.Indeed, this follows from g(w) → ∞ as |w| → ∞ and g((1 + σ)w 0 ) < 0 (unless w 0 = 0).We conclude that Finally, we can estimate Summarizing these estimates, we find that The right-hand side is bounded since w ∈ L ∞ (Ω; R n ) by assumption, implying that (ρ 1 , . . ., ρn−1 , θ) ∈ L ∞ (Ω; R n ).The first term on the left-hand side is bounded from below since, by definition (23) of h and Ee −w 0 = c w ρθ/θ 0 , Thus, we obtain a uniform bound for w in H 2 (Ω; R n ) and consequently also in W 1,4 (Ω; R n ).We can apply the Leray-Schauder fixed-point theorem to conclude the existence of a fixed point of S(•, 1).This, in turn, shows that w is a weak solution to the approximate problem (36)-(37).
Remark 6 (Treatment of the cross-terms).In the paper [15], the fluxes are given by A multiplication of this equation by ∇(µ/θ, −1/θ) shows that the cross-terms cancel out, since M is assumed to be positive semidefinite in [15].In the present work, we have and the cross-terms do not cancel.This is compensated by the sum n i,j=1 A ij /(m i m j ).Indeed, a computation shows that (also see (40)) since A is positive semidefinite because of (33).

Discrete entropy inequality.
We derive some estimates from (39) with σ = 1, which are uniform in (ε, τ ), by exploiting the sum I 4 + I 5 + I 8 , which we have neglected in (41).Taking into account that the estimate of I 10 becomes for σ = 1 we obtain the discrete entropy inequality where µ > 0 is defined in (22).
We deduce from Assumption (A4) that κ(e w )|∇w| 2 ≥ c κ |∇w| 2 , and in view of (42), this quantity is bounded in L 2 (Ω).Therefore, Lemma 7 yields a gradient bound for √ ρ i in Proof of Lemma 7. It follows from ( 24) and (29) that n i,j=1 and therefore, in view of the definition A ij = M BD ij √ ρ i ρ j and the positive definiteness (22) on the subspace L, n i,j=1 We insert the definition of the projection matrix P L : The last step follows from the pressure constraint (7).Indeed, by ( 8), ( 44) We have shown that n i,j=1 which equals (43) after integration over Ω.
Remark 8. We observe that the sum (44) vanishes even without requiring the constraint (7).Indeed, by (17), The fact that n j=1 d j vanishes is a necessary condition for the invertibility of the linear system (18).
In view of Lemma 7 and the lower bound κ ≥ c κ (1 + θ 2 ), we conclude from (42) the following discrete entropy inequality.
Lemma 9 (Discrete entropy inequality).It holds that Finally, we derive an estimate for the temperature.
Lemma 10.There exists a constant C > 0, only depending on λ, Ω, ∂Ω, and θ 0 such that Proof.We use θ as a test function in the approximate energy equation (37).Observing that We deduce from Young's inequality and Assumption (A4) on κ that Furthermore, J 3 ≥ 0. Definition (28) of B i and A ij as well as the bound ρ j ≤ ρ * show that The integrals J 5 are J 6 are bounded from below since and the dominant term in J 6 is θ 2 log θ, which is bounded from below by a negative constant.Finally, J 7 is nonnegative: Collecting these estimates finishes the proof.
The following lemma can be proved as in [15,Lemma 9].
It remains to pass to the limit (ε, τ ) → 0 in (57).We deduce from the strong convergence of (ρ (τ ) ) and (θ (τ ) ) that We deduce from the strong convergence ρ (τ ) i → ρ i in L q (Ω T ) for any q < ∞ and the boundedness of M BD ij that M BD ij (ρ (τ ) ) → M BD ij (ρ) strongly in any L q (Ω T ).In view of the weak convergences ∇ log θ (τ ) ⇀ ∇ log θ from (55) and ∇(ρ Hence, using ( 21), (ρ weakly in L 2 (Ω T ), where the last identity is the definition of u i .Then, taking into account the boundedness of ρ As the L 2 (Ω T ) norm is weakly lower semicontinuous, ) and, because of the uniform bounds, also in L 2 (Ω T ).Hence, Thus, applying the limit inferior (ε, τ ) → 0 to both sides of (57) yields the result.
The term K 1 can be rewritten as The algebraic system (5) with d i = ∇(ρ i θ)/m i can be formulated as This allows us to rewrite K 2 : Furthermore, it follows from We reformulate K 4 as b ij ρ i ρ j (u i − u j ) • (ū i − ūj )dxds =: K 41 + K 42 + K 43 .
A long but straightforward computation shows that Inserting these expressions into (59), putting K 12 on the left-hand side, and rearranging the terms, we find that H((ρ, θ)(t)|( ρ, θ)(t)) + 1 2 The second term on the left-hand side can be bounded from below.Indeed, it follows from the symmetry of (b ij ), definition (19) of M ij , and the positive definiteness (20) of M on L that where Y j = √ ρ j (u j − ūj ).The norm of the projection is computed according to where we used n i=1 ρ i u i = 0 in the third equality, and C 1 > 0 depends on ρ * and the L ∞ (Ω T ) norms of ūj , j = 1, . . ., n.Consequently, We turn to the estimation of the terms on the right-hand side of (60).By the Lipschitz continuity of κ and Young's inequality, K 11 is estimated as It remains to estimate the right-hand side of (58) in terms of the relative entropy.For this, we observe that, by [18,Lemma 16], Furthermore, for all functions f ∈ C 1 (R) with f ′ (1) = 0, Gronwall's lemma shows that H((ρ, θ)(t)|( ρ, θ)(t)) = 0 and hence ρ(t) = ρ(t) and θ(t) = θ(t) = 0 in Ω for t > 0. This finishes the proof.

1 j=1AB
where we have set θ = e w and θ = e w.Definition (25) of w i , definition(28) of B i , and the relations n−ij (w)∇w j = n j=1 A ij (w)∇q j , i (w)∇q i from (31)-(32) allow us to rewrite the sum I 4 + I 5 + I 8 as (40) I 4 + I 5 + I 8 = σ