Convergence of an implicit Euler Galerkin scheme for Poisson-Maxwell-Stefan systems
aa r X i v : . [ m a t h . NA ] S e p CONVERGENCE OF AN IMPLICIT EULER GALERKIN SCHEMEFOR POISSON–MAXWELL–STEFAN SYSTEMS
ANSGAR J ¨UNGEL AND OLIVER LEINGANG
Abstract.
A fully discrete Galerkin scheme for a thermodynamically consistent transientMaxwell–Stefan system for the mass particle densities, coupled to the Poisson equationfor the electric potential, is investigated. The system models the diffusive dynamics ofan isothermal ionized fluid mixture with vanishing barycentric velocity. The equationsare studied in a bounded domain, and different molar masses are allowed. The Galerkinscheme preserves the total mass, the nonnegativity of the particle densities, their bound-edness, and satisfies the second law of thermodynamics in the sense that the discreteentropy production is nonnegative. The existence of solutions to the Galerkin scheme andthe convergence of a subsequence to a solution to the continuous system is proved. Com-pared to previous works, the novelty consists in the treatment of the drift terms involvingthe electric field. Numerical experiments show the sensitive dependence of the particledensities and the equilibration rate on the molar masses. Introduction
The Maxwell–Stefan equations describe the dynamics of a fluid mixture in the diffusiveregime. They have numerous applications, for instance, in sedimentation, dialysis, electrol-ysis, and ion exchange. While Maxwell–Stefan models have been investigated since severaldecades from a modeling and simulation viewpoint in the engineering literature (e.g. [13]),the mathematical and numerical analysis started more recently [1, 16]. The global exis-tence of weak solutions under natural conditions was proved in [6, 21] for neutral mixtures.In case of ion transport, the electric charges and the self-consistent electric potential needto be taken into account. To our knowledge, no mathematical results are available in theliterature for such Poisson–Maxwell–Stefan models. In this paper, we prove the existenceof a weak solution to a structure-preserving fully discrete Galerkin scheme and its conver-gence to the continuous problem. This provides, for the first time, a global existence resultfor Poisson–Maxwell–Stefan systems.1.1.
Model equations.
We consider an ionized fluid mixture consisting of n componentswith the partial mass density ρ i , partial flux J i , and molar mass M i of the i th species. The Date : September 5, 2018.2000
Mathematics Subject Classification.
Key words and phrases.
Maxwell–Stefan systems, cross diffusion, ionized fluid mixtures, entropymethod, finite-element approximation, Galerkin method, numerical convergence.The authors acknowledge partial support from the Austrian Science Fund (FWF), P27352, P30000,F65, and W1245. evolution of the particle densities ρ i is governed by the partial mass balance equations(1) ∂ t ρ i + div J i = r i ( x ) , i = 1 , . . . , N, where r i are the production rates satisfying P ni =1 r i ( x ) = 0 and P ni =1 J i = 0. The molarconcentrations are defined by c i = ρ i /M i and x i = c i /c are the molar fractions, where c tot = P ni =1 c i denotes the total concentration and we have set x = ( x , . . . , x n ). Thepartial fluxes J i and the gradients of the molar fractions x i are related by the (scaled)Maxwell–Stefan equations(2) − N X j =1 k ij ( ρ j J i − ρ i J j ) = D i := ∇ x i + ( z i x i − ( ρ · x ) ρ i ) ∇ Φ , i = 1 , . . . , n. where k ij = k ji are the rescaled (reciprocal) Maxwell–Stefan diffusivities, D i is the drivingforce, z i the electric charge of the i th component, and Φ the electric potential. We refer toSection 2 for details on the modeling. These equations are coupled to the (scaled) Poissonequation(3) − λ ∆Φ = n X i =1 z i c i + f ( y ) , where λ is the scaled permittivity and f ( y ) is a fixed background charge. The equationsare solved in a bounded domain Ω ⊂ R d ( d ≥
1) and supplemented by the boundaryconditions J i · ν = 0 on ∂ Ω , i = 1 , . . . , n, (4) Φ = Φ D on Γ D , ∇ Φ · ν = 0 on Γ N , (5)where Γ D models the electric contacts, Γ N = ∂ Ω \ Γ D is the union of insulating boundarysegments, and ν denotes the exterior unit normal vector to ∂ Ω. This means that themixture cannot leave the container Ω and an electric field is applied at the contacts Γ N .The initial conditions are given by(6) ρ i ( · ,
0) = ρ i in Ω , i = 1 , . . . , n. We assume that the total mass is constant initially, P ni =1 ρ i = 1, which implies from(1) that the total mass is constant for all times, P ni =1 ρ i ( t ) = 1, expressing total massconservation.Observe that (2) defines a linear system in the diffusion fluxes. Since P ni =1 D i = 0,the kernel of that system is nontrivial, and we need to invert the relation between thefluxes J i and the driving forces D i on the orthogonal component of the kernel. It wasshown in [21, Section 2] that we can write (2) as D ′ = − A J ′ , where D ′ = ( D , . . . , D n − ), J ′ = ( J , . . . , J n − ), and A ∈ R ( n − × ( n − is invertible; see Section 3.1 for details. The n th components are recovered from D n = − P n − i =1 D i and J n = − P n − i =1 J i . Thus, (1) canbe written compactly as the cross-diffusion system [1, 21] ∂ t ρ ′ − div( A − D ′ ) = r ′ ( x ) , GALERKIN SCHEME FOR POISSON–MAXWELL–STEFAN SYSTEMS 3 where ρ ′ = ( ρ , . . . , ρ n − ). However, A − is not positive definite. To obtain a positivedefinite diffusion matrix, we need to transform the system. With the so-called entropyvariables(7) w i = log x i M i − log x n M n + (cid:18) z i M i − z n M n (cid:19) Φ , i = 1 , . . . , n − , we may formulate (1) as(8) ∂ t ρ ′ − div( B ∇ w ) = r ′ ( x ) , where B = ( B ij ) ∈ R ( n − × ( n − is symmetric and positive definite; see Section 3.1 fordetails. Here, ρ ′ and x are interpreted as (invertible) functions of w and Φ. This transfor-mation is well known in nonequilibrium thermodynamics, where w i is called the electro-chemical potential and B is the mobility or Onsager matrix.The transformation to entropy variables has two important advantages. First, introduc-ing the entropy(9) H ( ρ ) = Z Ω h ( ρ ) dy, h ( ρ ) = c tot n X i =1 x i log x i + λ |∇ (Φ − Φ D ) | , a formal computation shows that(10) dHdt + Z Ω ∇ w : B ∇ wdy = Z Ω n X i =1 r i ( x ) ∂h∂ρ i dy, if Φ D is constant, where A : B denotes the Frobenius matrix product between matrices A and B . (A discrete analog is shown in Theorem 1 below.) Thus, if the right-handside is nonpositive, the entropy t H ( ρ ( t )) is a Lyapunov functional and we may obtainsuitable estimates for w i . The entropy production (the diffusion term) is nonnegative,which expresses the second law of thermodynamics. This technique has been used in [6, 21]but without electric force terms. The derivation of gradient estimates is more delicate inthe presence of the electric potential; see Lemma 8. Second, the densities ρ i = ρ i ( w ) areautomatically positive and bounded and it holds that P ni =1 ρ i ( w ) = 1; see Corollary 7.This property is inherent of the transformation and it holds without the use of a maximumprinciple and independent of the functional setting.The aim of this paper is to extend the global existence result of [6, 21] to Maxwell–Stefansystems with electric forces and to suggest a fully discrete Galerkin scheme that preservesthe structure of the system, namely the nonnegativity of the particle densities, the L ∞ bound P ni =1 ρ i = 1, and a discrete analog of the entropy production inequality (10).1.2. State of the art.
Before presenting our main results, we briefly review the stateof the art of Maxwell–Stefan models. They were already derived in the 19th century byMaxwell using kinetic gas theory [25] and Stefan using continuum mechanics [32]. A moremathematical derivation from the Boltzmann equation can be found in [4, 15], includinga non-isothermal setting [19]. An advantage of the Maxwell–Stefan approach is that thedefinition of the driving forces can be adapted to the present physical situation, leading tovery general and thermodynamically consistent models [2].
A. J ¨UNGEL AND O. LEINGANG
When electrolytes are considered, we need to take into account the electric force. Usually,this is done in the context of Nernst–Planck models [27, 29], where the diffusion flux J i only depends on the density gradient of the i th component, thus without any cross-diffusion effects. Duncan and Toor [12] showed that cross-diffusion terms need to be takeninto account in a ternary gas. Dreyer et al. [11] outline some deficiencies of Nernst–Planck models and propose thermodynamically consistent Maxwell–Stefan type models.A numerical comparison between Nernst–Planck and Maxwell–Stefan models can be foundin [30].The first global-in-time existence result to the Maxwell–Stefan equations (1)-(2) withoutPoisson equation was proved by Giovangigli and Massot [16] for initial data around theconstant equilibrium state. The local-in-time existence of classical solutions was shownby Bothe [1]. The entropy structure of the Maxwell–Stefan system was revealed in [21],and a general global existence theorem could be shown. Further global existence resultscan be found in [18, 24]. The Maxwell–Stefan system was coupled to the heat equation[20] and to the incompressible Navier–Stokes equations [6]. In [15, Theorem 9.7.4] and [18,Theorem 4.3], the large-time asymptotics for initial data close to equilibrium was analyzed.The convergence to equilibrium for any initial data was investigated in [6, 21] withoutproduction terms and in [7] with production terms for reversible reactions. Salvarani andSoares proved a relaxation limit of the Maxwell–Stefan system to a system of linear heatequations [31].Surprisingly, there are not many papers concerned with numerical schemes which pre-serve the properties of the solution like conservation of total mass, nonnegativity, andentropy production. Many approximation schemes can be found in the engineering lit-erature, for instance finite-difference [22, 23] or finite-element [5] discretizations. In themathematical literature, finite-volume [28] and mixed finite-element [26] schemes as wellas explicit finite-difference schemes with fast solvers [14] were proposed. The existence ofdiscrete solutions was shown in [26], but only for ternary systems and under restrictionson the diffusion coefficients. The schemes of [3, 28] conserve the total mass, while thoseof [3, 8] also preserve the L ∞ bounds. The result of [8] is based on maximum principlearguments. Note that we are able to show the L ∞ bounds without the use of a maximumprinciple, as a result of the formulation in terms of entropy variables, and that we do notimpose any restrictions on the diffusivities (except positivity).All the cited results are concerned with the Maxwell–Stefan equations for neutral fluids,i.e. without electric effects. In this paper, we analyze for the first time Poisson–Maxwell–Stefan systems and show a discrete entropy production inequality. The cross-diffusionterms cause some mathematical difficulties which are not present in Nernst–Planck models.1.3. Main results.
Let ( θ ( k ) ) be an orthonormal basis of H D (Ω) and ( v ( k ) ) be an or-thonormal basis of H (Ω; R n − ) such that v ( k ) ∈ L ∞ (Ω; R n − ). We introduce the Galerkinspaces P N = span { u (1) , . . . , u ( N ) } , V N = span { v (1) , . . . , v ( N ) } . Furthermore, let
T > N ∈ N and set τ = T /N >
0. We impose the followingassumptions:
GALERKIN SCHEME FOR POISSON–MAXWELL–STEFAN SYSTEMS 5 (A1) Domain: Ω ⊂ R d is a bounded domain with Lipschitz boundary ∂ Ω = Γ D ∪ Γ N , whereΓ D ∩ Γ N = ∅ , Γ N is open in ∂ Ω, and meas(Γ D ) > ρ = ( ρ , . . . , ρ n ) is nonnegative and measurablesatisfying R Ω P ni =1 ρ i log ρ i dy < ∞ , ρ n = 1 − P n − i =1 ρ i ≥
0. The boundary dataΦ D ∈ H (Ω) ∩ L ∞ (Ω) solves − λ ∆Φ D = f in Ω and ∇ Φ D · ν = 0 on Γ N . Furthermore,let f ∈ L ∞ (Ω).(A3) Diffusion matrix: For any given ρ ∈ [0 , ∞ ) n , the transpose of the matrix A = ( A ij ) ∈ R n × n , defined by(11) A ij = (cid:26) P nℓ =1 , ℓ = i k iℓ ρ ℓ for i = j, − k ij ρ i for i = j, has the kernel ker( A ⊤ ) = span { } , where = (1 , . . . , ∈ R n .(A4) Production rates: The functions r i ∈ C ([0 , n ; R ) satisfy P ni =1 r i ( x ) log x i /M i ≤ x ∈ (0 , n , i = 1 , . . . , n .Assumptions (A1) and (A2) are rather natural. The condition ρ i log ρ i ∈ L (Ω) isneeded to apply the entropy method. By definition of A , it holds that ker( A ⊤ ) ⊂ span { } .If k ij > ρ j > { } = ker( A ⊤ ). For the generalcase k ij ≥
0, this property cannot be guaranteed and needs to be assumed. This explainsAssumption (A3). Assumption (A4) is needed to derive the entropy production inequality(10). It is satisfied for reversible reactions; see [7, Lemma 6].We consider the implicit Euler Galerkin scheme1 τ Z Ω (cid:0) ρ ′ ( u k + w D , Φ k ) − ρ ′ ( u k − + w D , Φ k − ) (cid:1) · φdy + ε Z Ω u k · φdy + Z Ω ∇ φ : B ( u k + w D , Φ k ) ∇ ( u k + w D ) dy = Z Ω r ′ ( x ( u k + w D , Φ k )) · φdy, (12) λ Z Ω ∇ Φ k · ∇ θdy = Z Ω (cid:18) n X i =1 z i c i ( u k + w D , Φ k ) + f ( y ) (cid:19) dy (13)for φ ∈ V N , θ ∈ P N , ε >
0, and we have defined(14) w D = ( w D, , . . . , w D,n − ) , w D,i = (cid:18) z i M i − z n M n (cid:19) Φ D . The discrete entropy variables are given by w k = u k + w D , and we used the notation c i ( w k , Φ k ) = ρ i ( w k , Φ k ) /M i , x i ( w k , Φ k ) = c i ( w k , Φ k ) /c k tot for i = 1 , . . . , n , and c k tot = P ni =1 ρ i ( w k , Φ k ) /M i .At time k = 0, we assume that ρ i ≥ η > w via definition(7). The condition can be removed by performing the limit η → ∈ H (Ω) ∩ L ∞ (Ω) be the unique solution to − λ ∆Φ = n X i =1 z i ρ i M i + f ( y ) in Ω , ∇ Φ · ν = 0 on Γ N , Φ = Φ D on Γ D . This defines ( w , Φ ). A. J ¨UNGEL AND O. LEINGANG
Theorem 1 (Existence for the Galerkin scheme) . Let Assumptions (A1)-(A4) hold. Thenthere exists a weak solution ( w k , Φ k ) ∈ V N × P N to (12) - (13) with w k = u k + w D , satisfying • preservation of L ∞ bounds: < ρ ki < for i = 1 , . . . , n ; • conservation of total mass: P ni =1 ρ ki = 1 in Ω ; • discrete entropy production inequality: H ( ρ k ) + τ Z Ω ∇ ( w k − w D ) : B ( w k , Φ k ) ∇ w k dy + ετ Z Ω | w k − w D | dy ≤ τ Z Ω n X i =1 z i M i r i ( x k )(Φ k − Φ D ) dy + H ( ρ k − ) , (15) where ρ k = ρ ( w k , Φ k ) . Theorem 1 is proved by using a fixed-point argument in the entropy variables. Using w k − w D as a test function in the fully discrete version of (8), we show in Section 4 that H ( ρ k ) + τ K Z Ω n X i =1 |∇ ( x ki ) / | dy + ετ Z Ω | w k − w D | dy ≤ τ K + H ( ρ k − ) , where K > ε is needed to conclude a uniform L estimate for w k , which is sufficientto apply the Leray-Schauder fixed-point theorem in the finite-dimensional Galerkin space.The ε -independent gradient estimate for x ki cannot be used since it does not give an estimatefor w ki (see (7)). It is possible to analyze system (12)-(13) for ε = 0 – see Step 2 of theproof of Theorem 3 –, but we lose the information about w k and obtain a solution in termsof ρ k . The term involving ε is technical and not essential for the numerical simulations (orthe structure preservation). However, we are not able to prove an existence result in termsof the entropy variable without such a regularization. Remark 2 (Conservation of partial mass) . When r i = 0, we have from (1) conservationof the partial mass k ρ i k L (Ω) . This conservation property does not hold exactly on thediscrete level because of the ε -regularization. It holds that for any δ >
0, there exists ε > < ε < ε ( ε is the value in (12)), (cid:12)(cid:12) k ρ ki k L (Ω) − k ρ i k L (Ω) (cid:12)(cid:12) ≤ δ k ρ i k L (Ω) , i = 1 , . . . , n − , (cid:12)(cid:12) k ρ kn k L (Ω) − k ρ n k L (Ω) (cid:12)(cid:12) ≤ δ n − X i =1 k ρ i k L (Ω) . The proof is the same as in [21, Theorem 4.1]. As δ > (cid:3)
Theorem 3 (Convergence of the Galerkin solution) . Let Assumptions (A1)-(A4) hold.Let ( ρ k , Φ k ) be a solution to (12) - (13) and set ρ τi ( y, t ) = ρ ki ( y ) , x τi ( y, t ) = x ki ( y ) , c τi ( y, t ) = c ki ( y ) , Φ τ ( y, t ) = Φ k ( y ) GALERKIN SCHEME FOR POISSON–MAXWELL–STEFAN SYSTEMS 7 for y ∈ Ω , t ∈ (( k − τ, kτ ] , i = 1 , . . . , n and introduce the shift operator ( σ τ ρ τi )( y, t ) = ρ k − i ( y ) for y ∈ Ω and t ∈ (( k − τ, kτ ] . Then there exist subsequences (not relabeled) suchthat, as ε → , N → ∞ , and τ → , ρ τi → ρ i strongly in L p (0 , T ; L p (Ω)) for any p < ∞ ,x τi ⇀ x i , Φ τ ⇀ Φ weakly in L (0 , T ; H (Ω)) ,τ − ( ρ τi − σ τ ( ρ τi )) ⇀ ∂ t ρ weakly in L (0 , T ; H (Ω) ′ ) , i = 1 , . . . , n, and the limit ( ρ, Φ) satisfies for all φ ∈ L (0 , T ; H (Ω; R n − )) and θ ∈ H D (Ω) , Z T h ∂ t ρ ′ , φ i dt + Z T Z Ω ∇ φ : A − ( ρ ) D ′ dydt = Z T Z Ω r ′ ( x ) · φdydt, (16) λ Z Ω ∇ Φ · ∇ θdy = Z Ω (cid:18) n X i =1 z i ρ i M i + f ( y ) (cid:19) θdy, (17) where D i = ∇ x i + ( z i x i − ( z · x ) ρ i ) ∇ Φ , ρ i = c tot M i x i , and c tot = P ni =1 ρ i /M i . Moreover, ρ n = 1 − P n − i =1 ρ i . In Theorem 3, h· , ·i denotes the duality bracket between H (Ω; R n − ) ′ and H (Ω; R n − ).The difficult part of the proof is the estimate of the diffusion term because of the contri-bution of the electric field. We show in Lemma 8 that Z Ω ∇ w k : B ∇ w k dy ≥ K Z Ω n X i =1 M / i | D ki | x ki dy ≥ K Z Ω n X i =1 |∇ ( x ki ) / | dy − K holds for some constants K , K , K >
0, which are independent of ε , N , and τ . Then theuniform L ∞ bound for x ki gives a uniform H (Ω) bound for x ki and consequently for ρ ki .Weak compactness allows us to pass to the limits ε → N → ∞ , and the limit τ → J i and the inversionof the map ρ w are presented in Section 3. Sections 4 and 5 are devoted to the proof ofthe main theorems. Finally, some numerical experiments are shown in Section 6.2. Modeling
We consider an isothermal electrolytic mixture of n fluid components in the boundeddomain Ω ⊂ R d ( d ≥
1) with boundary ∂ Ω. We assume that the mixture is not moving,so the barycentric velocity vanishes. The thermodynamic state of the mixture is describedby the partial mass densities ρ , . . . , ρ n and the electric field E . We suppose the quasi-static approximation E = −∇ Φ, where Φ is the electric potential. The evolution of themass densities ρ i = M i c i with the molar masses M i and molar concentrations (or numberdensities) c i is governed by the partial mass balances [10, (4)] ∂ t ρ i + div J i = r i ( x ) in Ω , t > , i = 1 , . . . , n, A. J ¨UNGEL AND O. LEINGANG where x = ( x , . . . , x n ) is the vector of molar fractions x i = ρ i / ( c tot M i ), c tot = P ni =1 c i isthe total concentration, J i the diffusion flux, and r i ( x ) the mass production rate of the i thspecies. We assume that the total flux and the total production vanishes, n X i =1 J i = 0 , n X i =1 r i ( x ) = 0 , which are necessary constraints to achieve total mass conservation, ∂ t P ni =1 ρ i = 0. Wesuppose that the total initial mass is constant in space, P ni =1 ρ i = ρ tot >
0, which impliesthat the total mass is constant in space and time, P ni =1 ρ i ( t ) = ρ tot for t > − ε (1 + χ )∆Φ = F n X i =1 z i c i + f ( y ) in Ω , where ε is the dielectric constant, χ the dielectric susceptibility, F the Faraday constant, z i the charge number of the i th species, and f ( y ) with y ∈ Ω models the charge of fixedbackground ions.The basic assumption of the Maxwell–Stefan theory is that the difference in speed andmolar fractions leads to a diffusion flux. They are implicitly given by the driving forces d i according to [2, (200)] − n X j =1 x j ( J i /M i ) − x i ( J j /M j ) c tot D ij = d i , i = 1 , . . . , n, where the numbers D ij = D ji are the Maxwell–Stefan diffusivities. Inserting the definition x i = ρ i / ( c tot M i ), we find that(18) − n X j =1 ρ j J i − ρ i J j c M i M j D ij = d i . In the present situation, the driving force is given by two components, the variation of thechemical potential µ i and the contribution of the body forces b i [2, (211)]: d i = c i M i RT ∇ µ i − ρ i RT ( b i − b tot ) , i = 1 , . . . , n, where R is the gas constant and T the (constant) temperature. Since ( D ij ) is symmetric,summing (18) from i = 1 , . . . , n leads to P ni =1 d i = 0. Furthermore, P ni =1 ∇ µ i vanishestoo; see below. This shows that b tot = ρ − P ni =1 ρ i b i . We assume that the only force is dueto the electric field (i.e., we neglect effects of gravity), b i = − ( z i /M i ) F ∇ Φ [30, (3)].It remains to determine the chemical potential. We define it by µ i = ∂h mix /∂ρ i , where h mix ( ρ ) = c tot RT ( P ni =1 x i log x i + 1) is the mixing free energy density [10, (23)]. Then µ i = 1 c tot M i ∂h mix ∂x i = RTM i (log x i + 1) , GALERKIN SCHEME FOR POISSON–MAXWELL–STEFAN SYSTEMS 9 and the driving force becomes d i = c i ∇ log x i + ρ i FRT M i (cid:18) z i − ρ tot n X j =1 z j ρ j M j (cid:19) ∇ Φ= c tot (cid:18) ∇ x i + FRT (cid:18) z i x i − ( z · x ) ρ i ρ tot (cid:19) ∇ Φ (cid:19) , (19)where z = ( z , . . . , z n ) and x = ( x , . . . , x n ). The Gibbs-Duhem equation n X i =1 ρ i ∂h mix ∂ρ i − h mix ( ρ ) = RT n X i =1 ρ i log x i + 1 M i − c tot RT (cid:18) n X i =1 x i log x i + 1 (cid:19) = 0shows that the pressure vanishes, which is consistent with our choice of the driving force (see[2, (211)]). The driving force in [30, (7)] contains a non-vanishing pressure that is relatedto our expression for the total body force. The resulting driving force (19), however, is thesame.We summarize the model equations: ∂ t ρ i + div J i = r i ( x ) , i = 1 , . . . , n, (20) − ε (1 + χ )∆Φ = F n X i =1 z i c i + f ( y ) , (21) − n X j =1 ρ j J i − ρ i J j c M i M j D ij = d i c tot = ∇ x i + FRT (cid:18) z i x i − ( z · x ) ρ i ρ tot (cid:19) ∇ Φ , (22)and the relations c i = ρ i M i , x i = ρ i c tot M i , c tot = n X i =1 c i . Equations (1)-(3) are obtained from (20)-(22) after setting λ = ε (1 + χ ) /F , k ij =1 / ( c M i M j D ij ), and D i = d i /c tot and after nondimensionalization. In particular, wescale the particle densities by ρ tot (then the scaled quantities satisfy P ni =1 ρ i = 1) and theelectric potential by F/ ( RT ). 3. Auxiliary results
We collect some auxiliary results needed for the existence analysis. The starting point isthe relation (2) below. Observe that the coefficients k ij depend on ρ i via c tot = P ni =1 ρ i /M i .This dependency does not complicates the analysis since the results in Section 3 holdpointwise for any given ρ i and c tot is uniformly bounded from above and below by1max i =1 ,...,n M i ≤ c tot = n X i =1 ρ i M i ≤ i =1 ,...,n M i . Expressions for the diffusion fluxes.
We review three different expressions forthe diffusion fluxes following [6, 21] and extend the formulas to electro-chemical potentials.We reformulate (2):(23) D i = − X j = i k ij ( ρ j J i − ρ i J j ) = X j = i k ij ρ i ρ j (cid:18) J i ρ i − J j ρ j (cid:19) . The symmetry of ( k ij ) implies that P ni =1 D i = 0. Compactly, we may write D = − AJ ,where D = ( D , . . . , D n ) ⊤ , J = ( J , . . . , J n ) ⊤ , and A = ( A ij ) with(24) A ij = (cid:26) P nℓ =1 , ℓ = i k iℓ ρ ℓ for i = j, − k ij ρ i for i = j. By Assumption (A3), it holds that im( A ) = ker( A ⊤ ) ⊥ = span { } ⊥ , where = (1 , . . . , ⊤ ∈ R n . We conclude from [21, Lemma 2.2] that all eigenvalues of e A := A | im( A ) are positiveuniformly in ρ ∈ [0 , n and that e A is invertible. Since P ni =1 J i = 0, each row of J =( J , . . . , J n ) is an element of im( A ), so the linear system D = − e AJ can be inverted,yielding J = − e A − D .We obtain another formulation by inverting the system in the first n − D ′ = ( D , . . . , D n − ) and J ′ = ( J , . . . , J n − ), we can write D ′ = − A J ′ , where the matrix A = ( A ij ) ∈ R ( n − × ( n − is defined by A ij = (cid:26) P n − ℓ =1 , ℓ = i ( k iℓ − k in ) ρ ℓ + k in if i = j, − ( k ij − k in ) ρ i if i = j. It is shown in [6, Lemma 4] that A is invertible and A − is bounded uniformly in ρ ∈ [0 , n .Thus, J ′ = − A − D ′ .Finally, we invert the relations (23). Using J n = − P n − i =1 J i , these relations (or theequivalent form D i = − P nj =1 A ij J j ) can be written as(25) D i ρ i − D n ρ n = − n − X j =1 C ij J j , where C ij = A ij ρ i − A in ρ i − A nj ρ n + A nn ρ n = − Y ij ρ i ρ j + Y in ρ i ρ n + Y nj ρ n ρ j − Y nn ρ n ,Y ij = (cid:26) P nℓ =1 , ℓ = i k iℓ ρ i ρ ℓ for i = j, − k ij ρ i ρ j for i = j. The matrix − Y = ( − Y ij ) ∈ R n × n is symmetric (since ( k ij ) is symmetric), quasi-positive,irreducible, and it has the strictly positive eigenvector with eigenvalue zero. Hence, bythe Perron-Frobenius theorem, the spectral bound of ( − Y ij ) is a simple eigenvalue (withvalue zero) and the spectrum of ( Y ij ) consists of numbers with positive real part and zero.Thus, Y is positive semidefinite. GALERKIN SCHEME FOR POISSON–MAXWELL–STEFAN SYSTEMS 11
We claim that the matrix C = ( C ij ) ∈ R ( n − × ( n − is positive definite on span { } ⊥ .Indeed, let y ∈ span { ρ } ⊥ . Then y · ρ = 0. Since · ρ = 1, we have y span { } = ker( Y ) andconsequently, span { ρ } ⊥ ⊂ ker( Y ) c . This means that − Y is negative definite on span { ρ } ⊥ .A computation shows that for any vector w = ( w , . . . , w n − ) ∈ R n − , it holds that n − X i,j =1 C ij w i w j = − n X i,j =1 Y ij ρ i ρ j e w i e w j where e w i = w i for i = 1 , . . . , n − e w n = − P n − i =1 w i . Then e w = ( e w , . . . , e w n ) ∈ span { } ⊥ . Since − Y is negative definite on span { ρ } ⊥ , we infer that ( − Y ij / ( ρ i ρ j )) is neg-ative definite on span { } ⊥ . Therefore, C is positive definite on span { } ⊥ . Its inverse B := c tot C − with B = ( B ij ) exists, only depends on the mass density vector ρ , and ispositive definite uniformly for all ρ ∈ [0 , n satisfying P ni =1 ρ i = 1 [6, Lemma 10]. Wededuce from (25) and (2) that J i = − n − X j =1 B ij (cid:18) D j ρ j − D n ρ n (cid:19) = − n − X j =1 B ij (cid:18) ∇ log x j M j − ∇ log x n M n + (cid:18) z j M j − z n M n (cid:19) ∇ Φ (cid:19) = − n − X j =1 B ij ∇ w j (26)for i = 1 , . . . , n − J n = − P n − i =1 J i , recalling definition (7) of w i . We summarize: Lemma 4 (Formulations of J i ) . Equations (23) can be written equivalently as J = − e A − D, J ′ = − A − D ′ , J ′ = − B ∇ w. The last expression for J i shows that the partial mass balances (1) can be formulated as ∂ t ρ ′ − div( B ∇ w ) = r ′ ( ρ ) , where ρ = ρ ( w ) and B = B ( ρ ( w )). By Definition (7), w is a function of ρ (and Φ). Theinverse relation ρ ( w ) is discussed in the following subsection.3.2. Inversion of ρ w . Definition (7) defines, for given Φ ∈ R , a mapping x w . Weclaim that this mapping can be inverted. If the molar masses are all the same, M := M i ,this can be done explicitly:(27) ρ i ( w ) = exp( M w i − ( z i − z n )Φ)1 + P n − j =1 exp( M w j − ( z j − z n )Φ) , i = 1 , . . . , n − , and ρ n = 1 − P n − i =1 ρ i . Unfortunately, when the molar masses are different, we cannotderive an explicit formula. Instead we adapt first Lemma 6 in [6]. Lemma 5 (Inversion of w and x ) . Let Φ ∈ R and define the function W Φ : (cid:26) x = ( x , . . . , x n ) ∈ (0 , n : n X i =1 x i = 1 (cid:27) → R n − by W Φ ( x ) = ( w ( x ) , . . . , w n − ( x )) , where w i ( x ) = log x i M i − log x n M n + (cid:18) z i M i − z n M n (cid:19) Φ , i = 1 , . . . , n − . Then W Φ is invertible and we can define x ′ ( w, Φ) := W − ( w ) and x n ( w, Φ) := 1 − P n − i =1 x i ,where x ′ ( w, Φ) = ( x , . . . , x n − ) .Proof. The proof is similar to that one of [6, Lemma 6]. Let w = ( w , . . . , w n − ) ∈ R n − and Φ ∈ R be given. Define the function f : [0 , → [0 , ∞ ) by f ( s ) = n − X i =1 (1 − s ) M i /M n exp (cid:20) M i w i − M i (cid:18) z i M i − z n M n (cid:19) Φ (cid:21) , s ∈ [0 , . Then f is continuous, strictly decreasing, and 0 = f (1) < f ( s ) < f (0) for s ∈ (0 , s ∈ (0 ,
1) such that f ( s ) = s . We define(28) x i = (1 − s ) M i /M n exp (cid:20) M i w i − M i (cid:18) z i M i − z n M n (cid:19) Φ (cid:21) > , i = 1 , . . . , n − . By definition, we have P n − i =1 x i = f ( s ) = s <
1. We set x n = 1 − s > P ni =1 x i = 1. Moreover, (28) can be written equivalently aslog x i M i + log(1 − s ) M n + (cid:18) z i M i − z n M n (cid:19) Φ = w i , and since 1 − s = x n , this shows that W − ( w ) = x ′ is the inverse mapping. (cid:3) Given ρ ∈ [0 , n , we know that x i = ρ i / ( c tot M i ) for i = 1 , . . . , n and P ni =1 x i = 1. Thisrelation can be inverted too. We recall [6, Lemma 7]: Lemma 6 (Inversion of ρ and x ) . Let x ′ ∈ (0 , n − and x n = 1 − P n − i =1 x i > be givenand define for i = 1 , . . . , n , ρ i ( x ′ ) = ρ i := c tot M i x i , where c tot = (cid:18) n X j =1 M j x j (cid:19) − . Then ρ = ( ρ , . . . , ρ n ) is the unique vector satisfying ρ n = 1 − P n − i =1 ρ i > , x i = ρ i / ( c tot M i ) for i = 1 , . . . , n , and c tot = P ni =1 ρ i /M i . Combining Lemmas 5 and 6, we conclude as in [6] that the mapping ρ w can beinverted. In fact, we just have to define ρ ′ = ρ ′ ( x ′ ( w, Φ)).
GALERKIN SCHEME FOR POISSON–MAXWELL–STEFAN SYSTEMS 13
Corollary 7 (Inversion of ρ and w ) . Let w = ( w , . . . , w n − ) ∈ R n − and Φ ∈ R be given.Then there exists a unique vector ρ = ( ρ , . . . , ρ n ) ∈ (0 , n satisfying P ni =1 ρ i = 1 suchthat (7) holds for ρ n = 1 − P n − i =1 ρ i and x i = ρ i / ( c tot M i ) with c tot = P ni =1 ρ i /M i . Themapping ρ ′ : R n − → (0 , n − , ρ ′ ( w, Φ) = ( ρ , . . . , ρ n − ) , is bounded. Proof of Theorem 1
Step 1: existence of solutions.
The idea is to apply the Leray-Schauder fixed-pointtheorem. We need to define the fixed-point operator. For this, let χ ∈ L ∞ (Ω; R n − ) and σ ∈ [0 , k − Φ D ∈ P N to the linear finite-dimensionalproblem λ Z Ω ∇ Φ k · ∇ θdy = Z Ω (cid:18) n X i =1 z i c i ( χ + w D , Φ k ) + f ( y ) (cid:19) θdy for all θ ∈ P N . In particular, Φ k ∈ L ∞ (Ω). Next, we wish to solve the linear finite-dimensional problem(29) a ( u, φ ) = σF ( φ ) for all φ ∈ V N , where a ( u, φ ) = Z Ω ∇ φ : B ( χ + w D , Φ k ) ∇ udy + ε Z Ω u · φdy,F ( φ ) = − τ Z Ω (cid:0) ρ ′ ( χ + w D , Φ k ) − ρ ′ ( u k − + w D , Φ k − ) (cid:1) dy + Z Ω r ′ ( x ( χ + w D , Φ k )) · φdy − Z Ω ∇ φ : B ( χ + w D , Φ k ) ∇ w D dy for u , φ ∈ V N . Since χ + w D ∈ L ∞ (Ω; R n − ) and Φ k ∈ L ∞ (Ω), Corollary 7 shows that ρ ( χ + w D , Φ k ) is bounded. We know from Section 3.1 that the matrix B = B ( χ + w D , Φ k )is positive definite and its elements are bounded. We deduce that the forms a and F arecontinuous on V N . Exploiting the equivalence of the norms in the finite-dimensional space V N , we find that a ( u, u ) ≥ ε k u k L (Ω) ≥ εK N k u k H (Ω) for some constant K N >
0, which implies that a is coercive on V N . By the Lax–Milgramlemma, there exists a unique solution u ∈ V N ⊂ L ∞ (Ω; R n − ) to (29) satisfying(30) εKN k u k L ∞ (Ω) ≤ a ( u, u ) = σF ( u ) ≤ K F k u k H (Ω) , and the constants K N and K F are independent of τ and σ . This defines the fixed-pointoperator S : L ∞ (Ω; R n − ) × [0 , → L ∞ (Ω; R n − ), S ( χ, σ ) = u . Standard arguments showthat S is continuous. Since V N is finite-dimensional, S is also compact. Furthermore, S ( χ,
0) = 0. Estimate (30) provides a uniform bound for all fixed points of S ( · , σ ). Thus,by the Leray-Schauder fixed-point theorem, there exists u k ∈ V N such that S ( u k ,
1) = u k ,and w k := u k + w D , Φ k solve (12)-(13). Step 2: proof of the discrete entropy production inequality (15) . We use the test function τ ( w k − w D ) ∈ V N in (12) and set ρ k := ρ ′ ( w k , Φ k ): Z Ω ( ρ k − ρ k − ) · ( w k − w D ) dy + τ Z Ω ∇ ( w k − w D ) : B ( w k , Φ k ) ∇ w k dy + ετ Z Ω | w k − w D | dy ≤ τ Z Ω r ′ ( x k ) · ( w k − w D ) dy. We claim that the first term on the left-hand side is the difference of the entropies at timesteps k and k −
1. To show this, we split the entropy density into two parts, h ( ρ k ) = h ( ρ k ) + h ( ρ k ), where h ( ρ k ) = c k tot n X i =1 x ki log x ki , h (Φ k ) = λ |∇ (Φ k − Φ D ) | , where we recall that x ki = ρ ki / ( c k tot M i ) and c k tot = P ni =1 ρ ki /M i . By the convexity of h , wehave h ( ρ k ) − h ( ρ k − ) ≤ ∂h ∂ρ ′ ( ρ k ) · ( ρ k − ρ k − ) = n X i =1 ( ρ ki − ρ k − i ) log x ki M i . Therefore, using ρ kn − ρ k − n = − P n − i =1 ( ρ ki − ρ k − i ), Z Ω (cid:0) h ( ρ k ) − h ( ρ k − ) (cid:1) dx ≤ Z Ω (cid:18) n − X i =1 ( ρ ki − ρ k − i ) log x ki M i + ( ρ kn − ρ k − n ) log x kn M n (cid:19) dy = Z Ω n − X i =1 ( ρ ki − ρ k − i ) (cid:18) log x ki M i − log x kn M n (cid:19) dy. (31)For the estimate of h , we first observe that n − X i =1 ( ρ ki − ρ k − i ) (cid:18) z i M i − z n M n (cid:19) = n − X i =1 ( ρ ki − ρ k − i ) z i M i + ( ρ kn − ρ k − n ) z n M n = n X n =1 ( ρ ki − ρ k − i ) z i M i . We infer from the Poisson equation (13) and Young’s inequality that Z Ω n − X i =1 ( ρ ki − ρ k − i ) (cid:18) z i M i − z n M n (cid:19) (Φ k − Φ D ) dy = Z Ω n X i =1 ( ρ ki − ρ k − i ) z i M i (Φ k − Φ D ) dy = Z Ω n X i =1 z i ( c ki − c k − i )(Φ k − Φ D ) dy = λ Z Ω ∇ (cid:0) (Φ k − Φ D ) − (Φ k − − Φ D ) (cid:1) (Φ k − Φ D ) dy GALERKIN SCHEME FOR POISSON–MAXWELL–STEFAN SYSTEMS 15 ≥ λ Z Ω |∇ (Φ k − Φ D ) | dy − λ Z Ω |∇ (Φ k − − Φ D ) | dy = Z Ω (cid:0) h (Φ k ) − h (Φ k − ) (cid:1) dy. (32)Taking into account the property r n ( ρ k ) = − P n − i =1 r i ( ρ k ), definition (7) of w ki , and As-sumption (A4), we compute Z Ω r ′ ( x k ) · ( w k − w D ) dy = Z Ω n − X i =1 r i ( x k ) (cid:18) log x ki M i − log x kn M n (cid:19) dy + Z Ω n − X i =1 r i ( x k ) (cid:18) z i M i − z n M n (cid:19) (Φ k − Φ D ) dy = Z Ω n X i =1 r i ( x k ) log x ki M i dy + Z Ω n X i =1 r i ( x k ) z i M i (Φ k − Φ D ) dy ≤ Z Ω n X i =1 r i ( x k ) z i M i (Φ k − Φ D ) dy, (33)Combining (31)-(33) gives the conclusion.5. Proof of Theorem 3
Let ( w k , Φ k ) be a weak solution to scheme (12)-(13) and define ρ k = ρ ( w k , Φ k ). Step 1: uniform estimates.
We derive estimates for ρ k and Φ k independent of ε , τ , and N . The starting point is the discrete entropy production inequality (15), and the maintask is to estimate the diffusion part. Lemma 8 (Estimate of the diffusion part) . There exist constants K > and K > ,both independent of ε , τ , and N , such that Z Ω ∇ ( w k − w D ) : B ∇ w k dy ≥ K n X i =1 k∇ ( x ki ) / k L (Ω) − K . Proof.
We drop the superindex k in the proof to simplify the notation. Recall that e A = A | im( A ) , where im( A ) = span { } ⊥ . We introduce as in the proof of Lemma 12 in [6]the symmetrization e A S = P − / e AP / , where P / = M / X / and M / := diag( M / ,. . . , M / n ), X / := diag( x / , . . . , x / n ). Then e A − S = P − / e A − P / is a self-adjoint en-domorphism whose smallest eigenvalue is bounded from below by some positive constantwhich depends only on ( k ij ).Since 0 = P ni =1 J i = P ni =1 ( B ∇ w ) i , we can express the last component in terms of theother components, ( B ∇ w ) n = − P n − i =1 ( B ∇ w ) i . Then ∇ w : B ∇ w = n − X i =1 (cid:26) ∇ log x i M i − ∇ log x n M n + (cid:18) z i M i − z n M n (cid:19) ∇ Φ (cid:27) · ( B ∇ w ) i = n − X i =1 M i ∇ (log x i + z i Φ) · ( B ∇ w ) i − M n ∇ (log x n + z n Φ) n − X i =1 ( B ∇ w ) i = n X i =1 M i ∇ (log x i + z i Φ) · ( B ∇ w ) i . To simplify the notation, we set Ψ i = ∇ (log x i + z i Φ) /M i , and Ψ = (Ψ , . . . , Ψ n ). ByLemma 4, B ∇ w = e A − D = P / e A − S P − / D . Hence, ∇ w : B ∇ w = Ψ : B ∇ w = Ψ : M / X / e A − S X − / M − / D = n X i,j =1 Ψ i M / i x / i ( e A − S ) ij x − / j M − / j D i = n X i,j =1 (cid:0) ∇ x / i + z i x / i ∇ Φ (cid:1) M − / i ( e A − S ) ij M − / j × (cid:0) ∇ x / j + ( z j x / j − ( x · z ) ρ j x − / j ) ∇ Φ (cid:1) . (34)In view of P ni =1 ( B ∇ w ) i = 0, it follows that n X i,j =1 (cid:0) M − / i x − / i ( z · x ) ρ i ∇ Φ (cid:1) ( e A S ) − ij M − / j (cid:0) ∇ x / j + ( z j x / j − ( x · z ) ρ j x − / j ) ∇ Φ (cid:1) = n X i,j =1 (cid:0) c ( z · x ) ∇ Φ (cid:1) e A − ij (cid:0) ∇ x j + ( z j x j − ( x · z ) ρ j ∇ Φ (cid:1) = (cid:0) c ( z · x ) ∇ Φ (cid:1) · n X i =1 ( B ∇ w ) i = 0 . Adding this expression to (34), we find that ∇ w : B ∇ w = n X i,j =1 M − / i (cid:0) ∇ x / i + ( z i x / i − ( z · x ) ρ i x − / i ∇ Φ (cid:1) ( e A S ) − ij M − / j × (cid:0) ∇ x / j + ( z j x / j − ( z · x ) ρ j x − / j ∇ Φ (cid:1) . The matrix e A − S is positive definite on im( e A S ) = span { ρ / } . As the vector (2 ∇ x / i +( z i x / i − ( x · z ) ρ i x − / i ∇ Φ) ni =1 lies in span { ρ / } , we obtain ∇ w : B ∇ w ≥ K B n X i =1 M − i (cid:12)(cid:12) ∇ x / i + ( z i x / i − ( x · z ) ρ i x − / i ∇ Φ (cid:12)(cid:12) ≥ K n X i =1 |∇ x / i | − K n X i =1 (cid:12)(cid:12) ( z i x / i − ( x · z ) ρ i x − / i ∇ Φ (cid:12)(cid:12) , GALERKIN SCHEME FOR POISSON–MAXWELL–STEFAN SYSTEMS 17 where K > K > M , . . . , M n . Since x i and ρ i x − / i = ρ / i / ( c tot M i )are bounded, the previous inequality becomes(35) ∇ w : B ∇ w ≥ K n X i =1 |∇ x / i | − K |∇ Φ | , where K depends on K and z i .In the following, let K > ε , n , and τ . We estimatethe expression involving the boundary term ∇ w D : B ∇ w = ∇ w D : A − D ′ = n − X i,j =1 ( A − ) ij (cid:18) z i M i − z n M n (cid:19) ∇ Φ D · (cid:0) ∇ x i + ( z i x i − ( z · x ) ρ i ) ∇ Φ (cid:1) ≤ Kδ + δ n − X i =1 (cid:12)(cid:12) ∇ x i + ( z i x i − ( z · x ) ρ i ) ∇ Φ (cid:12)(cid:12) , where K > ∇ Φ D , z i , M i , and A − . Since 0 ≤ x i ≤
1, we have |∇ x i | =4 x i |∇ x / i | ≤ |∇ x / i | and therefore,(36) ∇ w D : B ∇ w ≤ Kδ + 4 δ |∇ x / i | + δK |∇ Φ | . We infer from (35) and (36) that Z Ω ∇ ( w − w D ) : B ∇ wdy ≥ ( K − δ ) n X i =1 k∇ x / i k L (Ω) − K k∇ Φ k L (Ω) − Kδ .
By the boundedness of c i , the elliptic estimate for the Poisson equation gives(37) k Φ k H (Ω) ≤ K (1 + k c i k L (Ω) ) ≤ K. This proves the lemma. (cid:3)
Combining the discrete entropy inequality (15) and the estimate of Lemma 8 and sum-mation over k leads to the following result. Corollary 9.
There exist constants K > and K > , both independent of ε , n , and τ ,such that (38) H ( ρ k ) + τ K k X j =1 n X i =1 k∇ ( x ki ) / k L (Ω) + ετ k X j =1 k w j − w D k L (Ω) ≤ τ kK + H ( ρ ) . Step 2: limit ε → . For a fixed time step k , let ( w ε , Φ ε ) be a solution to (12)-(13) with ρ ε = ρ ( w ε , Φ ε ) and x εi = ρ εi / ( c ε tot M i ). Estimates (37) and (38) yield the following uniformbounds: k ρ εi k L ∞ (Ω) + k x εi k L ∞ (Ω) ≤ , i = 1 , . . . , n, (39) k x εi k H (Ω) + k Φ ε k H (Ω) + ε / k w εi k L (Ω) ≤ K, (40) where K > ε and N . The bound for x εi in H (Ω) is a consequence of thebound for ( x εi ) / in H (Ω) from (38) and the uniform L ∞ bound for x εi from (39). It followsthat c ε tot = P ni =1 ρ εi /M i is uniformly bounded in L ∞ (Ω). Moreover, because of P ni =1 ρ εi = 1, c ε tot ≥ (max i M i ) − > ρ εi = c ε tot M i x εi is uniformlybounded in H (Ω). Oberserving that the embedding H (Ω) ֒ → L (Ω) is compact, thereexist subsequences, which are not relabeled, such that as ε → x εi → x i , ρ εi → ρ i , Φ ε → Φ strongly in L (Ω) ,x εi ⇀ x i , ρ εi ⇀ ρ i , Φ ε ⇀ Φ weakly in H (Ω) ,εw εi → L (Ω) . In view of the L ∞ bounds for ( x εi ) and ( ρ εi ), the strong convergences for these (sub-)sequences hold in L p (Ω) for any p < ∞ . Consequently, c ε tot → c tot := P ni =1 ρ i /M i stronglyin L (Ω), and we can identify ρ i = c tot M i x i for i = 1 , . . . , n . Furthermore, c εi = ρ εi /M i → c i := ρ i /M i strongly in L (Ω) , i = 1 , . . . , n. Recalling definition (2) of D i , we have(41) D εi = ∇ x εi + ( z i x εi − ( z · x ε ) ρ εi ) ∇ Φ ε ⇀ D i := ∇ x i + ( z i x i − ( z · x ) ρ i ) ∇ Φweakly in L q (Ω) for any q < i = 1 , . . . , n . Since ( D εi ) is bounded in L (Ω), there existsa subsequence which converges to some function e D i weakly in L (Ω). By the uniquenessof the weak limits, we can identify e D i = D i . This shows that the convergence (41) holds in L (Ω). We deduce from the strong convergence of ( x εi ), the boundedness of ( x εi ) in L ∞ (Ω),and the continuity of r i that r i ( x ε ) → r i ( x ) strongly in L (Ω).We know from Lemma 4 that B ( w ε ) ∇ w ε = A − ( ρ ε )( D ε ) ′ . As A − ( ρ ) is uniformlybounded for ρ ∈ [0 , n and ( ρ ε ) converges strongly to ρ , we infer that A − ( ρ ε ) → A − ( ρ )strongly in L (Ω); the convergence holds even in every L p (Ω) for p < ∞ . Then, becauseof (41),(42) A − ( ρ ε )( D ε ) ′ ⇀ A − ( ρ ) D ′ weakly in L q (Ω) for all q < . In fact, since A − ( ρ ε )( D ε ) ′ ) is bounded in L (Ω) and thus (up to a subsequence) weaklyconverging in L (Ω), the convergence holds in L (Ω).These convergences are sufficient to perform the limit ε → ρ k , Φ k ) := ( ρ, Φ) solves1 τ Z Ω (cid:0) ( ρ k ) ′ − ( ρ k − ) ′ (cid:1) · φdy + Z Ω ∇ φ : A − ( ρ k ) ∇ ρ k dy = Z Ω r ′ ( x k ) · φdy, (43) λ Z Ω ∇ Φ k · ∇ θdy = Z Ω (cid:18) n X i =1 z i c ki + f ( y ) (cid:19) θdy (44)for all φ ∈ V N , θ ∈ P N . Step 3: limit N → ∞ . Let ( ρ N , Φ N ) be a solution to (43)-(44). Estimates (39)-(40)are independent of N . Thus, we can exactly argue as in step 2 and obtain limit functions( x, ρ, Φ) and c i = c tot M i x i for i = 1 , . . . , n as N → ∞ . These functions satisfy (43)-(44) GALERKIN SCHEME FOR POISSON–MAXWELL–STEFAN SYSTEMS 19 for all φ ∈ V N and θ ∈ P N and for all N ∈ N . The union of all V N is dense in H (Ω; R n − )and the union of all P N is dense in H D (Ω). Thus, by a density argument, system (43)-(44)holds for all test functions φ ∈ H (Ω; R n − ) and θ ∈ H D (Ω). Step 4: limit τ → . Let ( ρ k , Φ k ) be a solution to (43)-(44) with test functions φ ∈ H (Ω; R n − ) and θ ∈ H D (Ω). Then ρ ki = c k tot M i x ki and c ki = ρ ki /M i for i = 1 , . . . , n . We set ρ τi ( y, t ) = ρ ki ( y ) , x τi ( y, t ) = x ki ( y ) , c τi ( y, t ) = c ki ( y ) , Φ τ ( y, t ) = Φ k ( y )for y ∈ Ω, t ∈ (( k − τ, kτ ], i = 1 , . . . , n and introduce the shift operator ( σ τ ρ τ )( y, t ) = ρ τ ( y ) for y ∈ Ω and t ∈ (( k − τ, kτ ]. Finally, we set D τi = ∇ x τi + ( z i x τi − ( z · x τ ) ρ τi ) ∇ Φ τ and T = mτ for some fixed m ∈ N . Then we can write system (43)-(44) as1 τ Z T Z Ω (cid:0) ( ρ τ ) ′ − σ τ ( ρ τ ) ′ (cid:1) · φdydt + Z T Z Ω ∇ φ : A − ( ρ τ )( D τ ) ′ dydt = Z t Z Ω r ′ ( x τ ) · φdydt, (45) λ Z Ω ∇ Φ τ · ∇ θdy = Z Ω (cid:18) n X i =1 z i c τi + f ( y ) (cid:19) θdy (46)for all piecewise constant functions φ : (0 , T ) → H (Ω; R n − ) and θ : (0 , T ) → H D (Ω). Theentropy inequality (38), formulated in terms of ( ρ τ , Φ τ ), provides us with further uniformbounds since the right-hand side of (38) does not depend on τ : k ρ τi k L ∞ (Ω T ) + k x τi k L ∞ (Ω T ) ≤ K, (47) k ρ τi k L (0 ,T ; H (Ω)) + k x τi k L (0 ,T ; H (Ω)) + k Φ τ k L (0 ,T ; H (Ω)) ≤ K, (48)where we have set Ω T = Ω × (0 , T ). As a consequence, ( D τi ) is bounded in L (0 , T ; H (Ω)).It remains to derive a uniform estimate for the discrete time derivative of ρ τ . Takinginto account the uniform bound for A − ( ρ τ ), it follows that1 τ (cid:12)(cid:12)(cid:12)(cid:12) Z t Z Ω (cid:0) ( ρ τ ) ′ − σ τ ( ρ τ ) ′ (cid:1) · φdydt (cid:12)(cid:12)(cid:12)(cid:12) ≤ Z T k∇ φ k L (Ω) k A − ( ρ τ ) k L ∞ (Ω) k ( D τ ) ′ k L (Ω) dt + Z T k r ′ ( x τ ) k L (Ω) k φ k L (Ω) dt ≤ C k φ k L (0 ,T ; H (Ω)) . As the piecewise constant functions φ : (0 , T ) → H (Ω; R n − ) are dense in L (0 , T ; H (Ω; R n − )), this estimate also holds for all φ ∈ L (0 , T ; H (Ω; R n − )), and we concludethat τ − (cid:13)(cid:13) ( ρ τ ) ′ − σ τ ( ρ τ ) ′ (cid:13)(cid:13) L (0 ,T ; H (Ω) ′ ) ≤ K, i = 1 , . . . , n − . This estimate also holds for i = n since ρ τn = 1 − P n − i =1 ρ τi .By the Aubin-Lions lemma in the version of [9], there exists a subsequence of ( ρ τ ) whichis not relabeled such that, as τ → ρ τi → ρ i strongly in L (Ω T ) , i = 1 , . . . , n. In view of the L ∞ bound (47) for ρ τ , this convergence also holds in L p (Ω T ) for any p < ∞ .Furthermore, by (48), we have up to subsequences, x τi ⇀ x i , Φ τ ⇀ Φ weakly in L (0 , T ; H (Ω)) ,τ − ( ρ τi − σ τ ( ρ τi )) ⇀ ∂ t ρ i weakly in L (0 , T ; H (Ω) ′ ) . In particular, D τi ⇀ D i weakly in L (Ω T ), and we can identify D i = ∇ x i + ( z i x i − ( z · x ) ρ i ) ∇ Φ. The strong convergence of ( ρ τ ) and the weak convergence of ( D τi ) imply that A − ( ρ τ )( D τ ) ′ ⇀ A − ( ρ ) D ′ weakly in L q (Ω T ) , q < . Again, since ( A − ( ρ τ )( D τ ) ′ ) is bounded in L (Ω T ), this convergence holds in L (Ω T ). Fur-thermore, r ′ ( x τ ) → r ′ ( x ) strongly in L (Ω T ). Therefore, we can pass to the limit τ → ρ i ≥ η > ρ i ≥ η →
0. This is carried out in [6, Section 3.2] and we refer to this reference for details.6.
Numerical experiments
In this section, some numerical experiments based on scheme (12)-(13) in one spacedimension are presented.6.1.
Discretization and iteration procedure.
Let Ω = (0 ,
1) be divided into n p ∈ N uniform subintervals of length h = 1 /n p . We use uniform time steps with time step size τ > w, Φ), the molar fractions x i are computed from thefixed-point problem (see the proof of Lemma 5)(49) f ( s ) = n − X i =1 (1 − s ) M i /M n exp (cid:20) M i w i − M i (cid:18) z i M i − z n M n (cid:19) Φ (cid:21) , s ∈ [0 , , with unique solution s ∈ (0 , x i = (1 − s ) M i /M n exp (cid:20) M i w i − M i (cid:18) z i M i − z n M n (cid:19) Φ (cid:21) , i = 1 , . . . , n − , and x n = 1 − s . Then we set (see Lemma 6) c tot = P ni =1 ( M i x i ) − and ρ i = c tot M i x i for i = 1 , . . . , n .Instead of solving the nonlinear discrete system (12)-(13) by a full Newton method, weemploy a linearized semi-implicit approach, i.e., we linearize ρ ( w, Φ) and use the previoustime step in the diffusion matrix B ( w ). More precisely, let w ∈ V N and Φ ∈ P N be given.We linearize ρ ( w, Φ) by ρ ( w, Φ) + ∇ ( w, Φ) ρ ′ ( w, Φ) · ( w − w, Φ − Φ) . This leads to the problem in the variable ζ = ( w − w, Φ − Φ):(50) L ( ζ , φ ) = F ( φ ) , K ( ζ n , θ ) = G ( θ ) for all φ ∈ V N , θ ∈ P N , GALERKIN SCHEME FOR POISSON–MAXWELL–STEFAN SYSTEMS 21 where L ( ζ , φ ) = Z Ω ∇ ( w, Φ) ρ ′ ( w, Φ) · ( ζ , φ ) dy + τ Z Ω ∂ x φ · B ( w, Φ) ∂ x ζ dy + ετ Z Ω ( ζ − w D ) · φdy,F ( φ ) = − Z Ω (cid:0) ρ ′ ( w, Φ) − ρ ′ ( w k − , Φ k − ) (cid:1) · φdy − τ Z Ω ∂ x φ · B ( w, Φ) ∂ x wdy,K ( ζ n , θ ) = λ Z Ω ∂ x ζ n ∂ x φdy − Z Ω n X i =1 z i M i ∇ ( w, Φ) ρ i ( w, Φ) · ζ θdy,G ( θ ) = − λ Z Ω ∂ x Φ ∂ x θdy + Z Ω (cid:18) n X i =1 z i ρ i ( w, Φ) M i + f ( y ) (cid:19) θdy. The iteration with starting point ( w (0) h , Φ (0) h ) := ( w k − , Φ k − ) is then defined by ( w ( m +1) h , Φ ( m +1) h ) := ( w, Φ) + ζ for m ≥
0. The iteration stops when k ζ k ℓ ∞ < ε tol for some tolerance ε tol > m ≥ m max for a maximal number of iterations. We summarize the scheme inAlgorithm 1. Algorithm 1 (Pseudo-code for the finite-element scheme in entropy variables.) procedure Maxwell-Stefan system in entropy variables Set ( w (0) h , Φ (0) h ) = ( w k − , Φ k − ), ρ (0) h = ρ ′ ( w h , Φ h ), x (0) h = ρ (0) h / ( M i c (0) h ), c (0) h = P ni =1 ( ρ (0) h ) i /M i , m = 0, ε >
0, and m max . while err > ε do Solve linear system (50) with solution ζ . Set ( w ( m +1) h , Φ ( m +1) h ) := ( w mh , Φ mh ) + ζ . Solve the fixed-point problem (49) with solution s . Compute x ( m +1) h and ρ ( m +1) h . Set err := k ( w ( m +1) h , Φ ( m +1) h ) − ( w ( m ) h , Φ ( m ) h ) k ℓ ∞ . ( m + 1) ← ( m ). if m > m max or err < ε then Break end if end while end procedure
The linear system (50) and the fixed-point problem (49) are solved using MATLAB.We choose the numerical parameters h = 10 − , τ = 10 − , ε tol = 10 − , and ε = 2 − ≈ . · − (the scheme works also for ε = 0).We have compared our results with the solutions from a finite-element scheme de-rived from the original system in the variables ρ i and a Picard iteration procedure forthe nonlinear discrete system. It turned out that the results are basically the same, i.e. k ρ i − ρ i ( w, Φ) k L ∞ (Ω) ≤ − . Numerical examples.
In all numerical examples, we neglect reaction terms andchoose the diffusivities according to [3, 13]: D = 0 . D = 0 . D = 0 . n = 3. The charges are given by z = z = 1 and z = 0 and the initial data is definedas in [3]: ρ ( y ) = . y < . , − . − η ) y − . η − (0 . · . . ≤ y < . ,η for 0 . ≤ y ≤ η = 10 − , ρ ( y ) = 0 .
2, and ρ ( y ) = (1 − ρ − ρ )( y ) for y ∈ Ω = (0 , y ) = 0 for y ∈ { , } . The dynamics of the particle densitiesand the electric potential are shown in Figure 1. The solution at time t = 17 is essentiallystationary and, in fact, in equilibrium. Because of the choice of the parameters, thestationary solution is symmetric around x = .The situation changes drastically when the molar masses are different (example 2). Fig-ure 2 shows the stationary solutions with the same parameters as in the previous exampleexcept M = 6. Here, the discrete relative entropy is defined by H ∗ ( ρ kh ) = Z (cid:18) c k tot ,h n X i =1 ( x kh ) i log ( x kh ) i ( x ∞ h ) i + λ |∇ (Φ kh − Φ ∞ h ) | (cid:19) dy, where ( ρ kh , Φ kh ) is the finite-element solution at time kτ and ( x ∞ h , Φ ∞ h ) is the stationarysolution. The integral and gradients are computed by the trapezoidal and gradient routinesof MATLAB. The semi-logarithmic plot of the relative entropy shows that the entropyconverges to zero exponentially fast.For example 3, we choose the same initial conditions and parameters as before, butwe take non-equilibrium boundary data Φ(0) = 10, Φ(1) = 0. The solutions at time t = 8 for various molar masses M are displayed in Figure 3. Since ρ and ρ have bothpositive charge and the potential on the left boundary is positive, both species avoid theleft boundary and move to the right.In example 4, we interchange the roles of M and M , i.e., we choose M = 1 and M ∈ { , , } . We observe in Figure 4 that the first species is more concentrated at theright boundary while in the previous example, this holds true for the second species.The previous examples show that the convergence rate to equilibrium strongly dependson the ratio of the molar masses. It turns out that this effect is triggered by the drift term,and without electric field, the convergence rates are similar for different molar masses. Thisbehavior can be observed in Figure 5 (example 5), where we have taken the same parametersas in the previous example but neglect the electric field. In this situation, the steady stateis constant in space and explicitly computable; indeed, we have ρ ∞ i = mean(Ω) − k ρ i k L (Ω) .Note that the steady state in the previous examples is not constant.Finally, we compute the numerical convergence rate when the grid size tends to zero forthe situation of example 3 (non-equilibrium boundary conditions for the potential). Wechoose the time t = 0 .
01 and the time step size τ = 10 − . The solutions are computed on GALERKIN SCHEME FOR POISSON–MAXWELL–STEFAN SYSTEMS 23 . . . . . y ρ t = 0 t = 0 . t = 17 0 0 . . . . . y ρ t = 0 t = 0 . t = 170 0 . . . . . y ρ t = 0 t = 0 . t = 17 . y Φ t = 0 t = 0 . t = 17 Figure 1.
Example 1: Particle densities ρ i and electric potential for molarmasses M = M = M = 1 versus position at various times. The boundaryconditions for the electric potential are in equilibrium.nested meshes with grid sizes h ∈ { . , . , . , . , . } and compared to thereference solution, computed on a very fine mesh with 25601 elements ( h ≈ · − ). Asexpected, we observe a second-order convergence in space; see Figure 6. References [1] D. Bothe. On the Maxwell-Stefan equations to multicomponent diffusion. In: J. Escher et al. (eds).
Parabolic Problems. Progress in Nonlinear Differential Equations and their Applications , pp. 81-93.Springer, Basel, 2011.[2] D. Bothe and W. Dreyer. Continuum thermodynamics of chemically reacting fluid mixtures.
ActaMech.
226 (2015), 1757–1805.[3] L. Boudin, B. Grec, and F. Salvarani. A mathematical and numerical analysis of the Maxwell–Stefandiffusion equations.
Discrete Cont. Dyn. Sys. B
17 (2012), 1427-1440. . . . . . y ρ M = 1 M = 6 0 0 . . . . . y ρ M = 1 M = 60 0 . . . . . y ρ M = 1 M = 6 − − t l og H ∗ ( ρ k h ) M = 1 M = 6 Figure 2.
Example 2: Particle densities ρ i at time t = 4 versus position andrelative entropy (bottom right) for molar masses M = 6 and M = M = 1.The boundary conditions for the electric potential are in equilibrium. [4] L. Boudin, B. Grec, M. Pavi´c, and F. Salvarani. Diffusion asymptotics of a kinetic model for gaseousmixtures. Kinetic Related Models
Intern. J. Numer. Meth. Engrg.
73 (2008), 665-684.[6] X. Chen and A. J¨ungel. Analysis of an incompressible Navier–Stokes–Maxwell–Stefan system.
Com-mun. Math. Phys.
340 (2015), 471–497[7] E. Daus, A. J¨ungel, and B.-Q. Tang. Exponential time decay of solutions to reaction-cross-diffusionsystems of Maxwell–Stefan type. Submitted for publication, 2018. arXiv:1802.10274.[8] K. Dieter-Kisling, H. Marschall, and D. Bothe. Numerical method for coupled interfacial surfactanttransport on dynamic surface meshes of general topology.
Computers & Fluids
109 (2015), 168-184.[9] M. Dreher and A. J¨ungel. Compact families of piecewise constant functions in L p (0 , T ; B ). Nonlin.Anal.
75 (2012), 3072-3077.
GALERKIN SCHEME FOR POISSON–MAXWELL–STEFAN SYSTEMS 25 . . . . . y ρ M = 2 M = 4 M = 6 0 0 . . . . . y ρ M = 2 M = 4 M = 60 0 . . . . . y ρ M = 2 M = 4 M = 6 − − t l og H ∗ ( ρ k h ) M = 2 M = 4 M = 6 Figure 3.
Example 3: Particle densities ρ i at time t = 8 versus position andrelative entropy (bottom right) for various molar masses M . The boundaryconditions for the electric potential are not in equilibrium. [10] W. Dreyer, P.-E. Druet, P. Gajewski, and C. Guhlke. Analysis of improved Nernst-Planck-Poissonmodels of compressible isothermal electrolytes. Part I: Derivation of the model and survey of theresults. WIAS Berlin, Germany, preprint no. 2395, 2017.[11] W. Dreyer, C. Guhlke, and R. M¨uller. Overcoming the shortcomings of the Nernst-Planck-Poissonmodel. Phys. Chem. Chem. Phys.
15 (2013), 7075-7086.[12] J. Duncan and H. Toor. An experimental study of three component gas diffusion.
AIChE J.
Fluid Dynam.
27 (1993), 590-596. Translated from em Izv. Ross. Akad. Nauk Mekh.Zhidk. Gaza 4 (1992), 178-185 (Russian).[14] J. Geiser. Iterative solvers for the Maxwell-Stefan diffusion equations: Methods and applications inplasma and particle transport.
Cogent Math.
Multicomponent Flow Modeling . Birkh¨auser, Basel, 1999. . . . . . y ρ M = 2 M = 4 M = 6 0 0 . . . . . y ρ M = 2 M = 4 M = 60 0 . . . . . y ρ M = 2 M = 4 M = 6 − − t l og H ∗ ( ρ k h ) M = 2 M = 4 M = 6 Figure 4.
Example 4: Particle densities ρ i at time t = 8 versus position andrelative entropy (bottom right) for various molar masses M . The boundaryconditions for the electric potential are not in equilibrium. [16] V. Giovangigli and M. Massot. The local Cauchy problem for multicomponent flows in full vibrationalnon-equilibrium. Math. Meth. Appl. Sci.
21 (1998), 1415-1439.[17] V. Giovangigli and M. Massot. Asymptotic stability of equilibrium states for multicomponent reactiveflows.
Math. Models Meth. Appl. Sci.
Nonlin. Anal.
159 (2017), 264–284.[19] H. Hutridurga and F. Salvarani. Maxwell–Stefan diffusion asymptotics for gas mixtures in non-isothermal setting.
Nonlin. Anal.
159 (2017), 285-297.[20] H. Hutridurga and F. Salvarani. Existence and uniqueness analysis of a non-isothermal cross-diffusionsystem of Maxwell-Stefan type.
Appl. Math. Lett.
75 (2018), 108-113.[21] A. J¨ungel and I. Stelzer. Existence analysis of Maxwell-Stefan systems for multicomponent mixtures.
SIAM J. Math. Anal.
45 (2013), 2421-2440.[22] E. Leonardia and C. Angeli. On the Maxwell–Stefan approach to diffusion: a general resolution in thetransient regime for one-dimensional systems.
J. Phys. Chem. B
114 (2010), 151-164.
GALERKIN SCHEME FOR POISSON–MAXWELL–STEFAN SYSTEMS 27 − − t l og H ∗ ( ρ k h ) M = 2 M = 4 M = 6 Figure 5.
Example 5: Semi-logarithmic plot of the relative entropy H ∗ ( ρ kh )versus time, without electric potential and for different molar masses. [23] J.-P. Loos, P. Verheijen, and J. Moulin. Numerical simulation of the generalized Maxwell–Stefanmodel for multicomponent diffusion in microporous sorbents. Collect. Czech. Chem. Commun.
J. Math. Pures Appl.
104 (2015), 102-138.[25] C. Maxwell. On the dynamical theory of gases.
Phil. Trans. Roy. Soc. London
157 (1866), 49-88.[26] M. McLeod and Y. Bourgault. Mixed finite element methods for addressing multi-species diffusionusing the Maxwell-Stefan equations.
Comput. Meth. Appl. Mech. Engrg.
279 (2014), 515-535.[27] W. Nernst. Die elektromotorische Wirksamkeit der Ionen.
Z. Physikalische Chemie
J.Comput. Phys.
230 (2011), 3525-3537.[29] M. Planck. ¨Uber die Potentialdifferenz zwischen zwei verd¨unnten L¨osungen bin¨arer Electrolyte.
An-nalen der Physik
276 (1890), 561-576.[30] S. Psaltis and T. Farrell. Comparing charge transport predictions for a ternary electrolyte using theMaxwell–Stefan and Nernst–Planck equations.
J. Electrochem. Soc.
158 (2011), A33-A42.[31] F. Salvarani and J. Soares. On the relaxation of the Maxwell–Stefan system to linear diffusion. Sub-mitted for publication, 2018. hal-01791067.[32] J. Stefan. ¨Uber das Gleichgewicht und Bewegung, insbesondere die Diffusion von Gasgemengen.
Sitzungsberichte Kaiserl. Akad. Wiss. Wien
63 (1871), 63-124.[33] J. Wesselingh and R. Krishna.
Mass Transfer in Multicomponent Mixtures . Delft University Press,Delft, 2000.
Institute for Analysis and Scientific Computing, Vienna University of Technology,Wiedner Hauptstraße 8–10, 1040 Wien, Austria
E-mail address : [email protected] Institute for Analysis and Scientific Computing, Vienna University of Technology,Wiedner Hauptstraße 8–10, 1040 Wien, Austria
E-mail address : [email protected] − − − − − − . h i L − e rr o r ρ Order 1Order 2 t = 0 . − − − − − − h i L − e rr o r ρ Order 1Order 2 t = 0 . − − − − − − h i L − e rr o r ρ Order 1Order 2 t = 0 . − − − − − − h i L − e rr o r Φ Order 1Order 2 t = 0 . Figure 6.
Discrete L -error relative to the reference solution for the densi-ties and the potential (bottom right) at time t = 0 ..