Convergence of an implicit Euler Galerkin scheme for Poisson–Maxwell–Stefan systems

A fully discrete Galerkin scheme for a thermodynamically consistent transient Maxwell–Stefan system for the mass particle densities, coupled to the Poisson equation for the electric potential, is investigated. The system models the diffusive dynamics of an isothermal ionized fluid mixture with vanishing barycentric velocity. The equations are studied in a bounded domain, and different molar masses are allowed. The Galerkin scheme preserves the total mass, the nonnegativity of the particle densities, their boundedness and satisfies the second law of thermodynamics in the sense that the discrete entropy production is nonnegative. The existence of solutions to the Galerkin scheme and the convergence of a subsequence to a solution to the continuous system is proved. Compared to previous works, the novelty consists in the treatment of the drift terms involving the electric field. Numerical experiments show the sensitive dependence of the particle densities and the equilibration rate on the molar masses.


Introduction
The Maxwell-Stefan equations describe the dynamics of a fluid mixture in the diffusive regime. They have numerous applications, for instance, in sedimentation, dialysis, electrolysis, and ion exchange. While Maxwell-Stefan models have been investigated since several decades from a modeling and simulation viewpoint in the engineering literature (e.g., [17]), the mathematical and numerical analysis started more recently [1,20]. The global existence of weak solutions under natural conditions was proved in [8,26] for neutral mixtures. In the case of ion transport, the electric charges and the self-consistent electric potential need to be taken into account. To our knowledge, no mathematical results are available in the literature for such Poisson-Maxwell-Stefan models. In this paper, we prove the existence of a weak solution to a structure-preserving fully discrete Galerkin scheme and its convergence to the continuous problem. This provides, for the first time, a global existence result for Poisson-Maxwell-Stefan systems.

Model equations
We consider an ionized fluid mixture consisting of n components with the partial mass density ρ i , partial flux J i , and molar mass M i of the ith species. The evolution of the particle densities ρ i is governed by the partial mass balance equations where r i are the production rates satisfying k ij (ρ j J i − ρ i J j ) = D i := ∇x i + (z i x i − (z · x)ρ i )∇ , i = 1, . . . , n. (2) where k ij = k ji are the rescaled (reciprocal) Maxwell-Stefan diffusivities, D i is the driving force, z i the electric charge of the ith component, and the electric potential. We refer to Section 2 for details on the modeling. These equations are coupled to the (scaled) Poisson equation where λ is the scaled permittivity, and f (y) is a fixed background charge. The equations are solved in a bounded domain ⊂ R d (d ≥ 1) and supplemented by the boundary conditions J i · ν = 0 on ∂ , i = 1, . . . , n, where D models the electric contacts, N = ∂ \ D is the union of insulating boundary segments, and ν denotes the exterior unit normal vector to ∂ . This means that the mixture cannot leave the container and an electric field is applied at the contacts D . The initial conditions are given by We assume that the total mass is constant initially, n i=1 ρ 0 i = 1, which implies from (1) that the total mass is constant for all times, n i=1 ρ i (t) = 1, expressing total mass conservation.
Observe that (2) defines a linear system in the diffusion fluxes. Since n i=1 D i = 0, the kernel of that system is nontrivial, and we need to invert the relation between the fluxes J i and the driving forces D i on the orthogonal component of the kernel. It was shown in [26, Section 2] that we can write (2) as D = −A 0 J , where D = (D 1 , . . . , D n−1 ), J = (J 1 , . . . , J n−1 ), and A 0 ∈ R (n−1)×(n−1) is invertible; see Section 3.1 for details. The nth components are recovered from D n = − n−1 i=1 D i and J n = − n−1 i=1 J i . Thus, (1) can be written compactly as the cross-diffusion system [1,26] ∂ t ρ − div(A −1 0 D ) = r (x), where ρ = (ρ 1 , . . . , ρ n−1 ). However, A −1 0 is not positive definite. To obtain a positive definite diffusion matrix, we need to transform the system. With the so-called entropy variables we may formulate (1) as where B = (B ij ) ∈ R (n−1)×(n−1) is symmetric and positive definite; see Section 3.1 for details. Here, ρ and x are interpreted as (invertible) functions of w and . This transformation 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, introducing the entropy a formal computation shows that if D is constant and ∇w : B∇w = n−1 i,j =1 B ij ∇w i · ∇w j . (A discrete analog is shown in Theorem 1 below.) Thus, if the right-hand side is nonpositive, the entropy t → H (ρ(t)) is a Lyapunov functional and we may obtain suitable 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 [8,26] but without electric force terms. The derivation of gradient estimates is more delicate in the presence of the electric potential; see Lemma 8. Second, the densities ρ i = ρ i (w) are automatically positive and bounded and it holds that n i=1 ρ i (w) = 1; see Corollary 7. This property is inherent of the transformation and it holds without the use of a maximum principle and independent of the functional setting.
The aim of this paper is to extend the global existence result of [8,26] to Maxwell-Stefan systems with electric forces and to suggest a fully discrete Galerkin scheme that preserves the structure of the system, namely the nonnegativity of the particle densities, the L ∞ bound n i=1 ρ i = 1, and a discrete analog of the entropy production inequality (10).

State of the art
Before presenting our main results, we briefly review the state of the art of Maxwell-Stefan models. They were already derived in the nineteenth century by Maxwell using kinetic gas theory [30] and Stefan using continuum mechanics [37]. A more mathematical derivation from the Boltzmann equation can be found in [4,19], including a non-isothermal setting [22]. An advantage of the Maxwell-Stefan approach is that the definition of the driving forces can be adapted to the present physical situation, leading to very general and thermodynamically consistent models [2].
When electrolytes are considered, we need to take into account the electric force. Usually, this is done in the context of Nernst-Planck models [32,34], where the diffusion flux J i only depends on the density gradient of the ith component, thus without any cross-diffusion effects. Duncan and Toor [15] showed that cross-diffusion terms need to be taken into account in a ternary gas. Dreyer et al. [14] 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 found in [35].
The first global-in-time existence result to the Maxwell-Stefan equations (1)-(2) without Poisson equation was proved by Giovangigli and Massot [20] for initial data around the constant equilibrium state. The local-in-time existence of classical solutions was shown by Bothe [1]. The entropy structure of the Maxwell-Stefan system was revealed in [26], and a general global existence theorem could be shown. Further global existence results can be found in [21,29]. The Maxwell-Stefan system was coupled to the heat equation [23] and to the incompressible Navier-Stokes equations [8]. In [19,Theorem 9.7.4] and [21,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 [8,26] without production terms and in [9] with production terms for reversible reactions. Salvarani and Soares proved a relaxation limit of the Maxwell-Stefan system to a system of linear heat equations [36].
Surprisingly, there are not many papers concerned with numerical schemes which preserve the properties of the solution like conservation of total mass, nonnegativity, and entropy production. Many approximation schemes can be found in the engineering literature, for instance finite-difference [27,28] or finite-element [6] discretizations. In the mathematical literature, finite-volume [33] and mixed finiteelement [31] schemes as well as explicit finite-difference schemes with fast solvers [18] were proposed. The existence of discrete solutions was shown in [31], but only for ternary systems and under restrictions on the diffusion coefficients. The schemes of [3,33] conserve the total mass, while those of [3,11] also preserve the L ∞ bounds. The result of [11] is based on maximum principle arguments. Note that we are able to show the L ∞ bounds without the use of a maximum principle, as a result of the formulation in terms of entropy variables, and that we do not impose 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-diffusion terms cause some mathematical difficulties which are not present in Nernst-Planck models.

Main results
Set H 1 D ( ) = {u ∈ H 1 ( ) : u = 0 on D } and let (θ (k) ) be a basis of H 1 D ( ) and (v (k) ) be a basis of H 1 ( ; R n−1 ) such that v (k) ∈ L ∞ ( ; R n−1 ). We introduce the Galerkin spaces (A3) Diffusion matrix: For any given ρ ∈ [0, ∞) n satisfying n i=1 ρ i = 1, the transpose of the matrix A = (A ij ) ∈ R n×n , defined by has the kernel ker(A ) = span{1}, where 1 = (1, . . . , 1) ∈ R n . (A4) Production rates: The functions r i Assumptions (A1) and (A2) are rather natural. The condition ρ i log ρ i ∈ L 1 ( ) is needed to apply the entropy method. By definition of A, it holds that ker(A ) ⊂ span{1}. If k ij > 0 (and ρ j > 0), a computation shows that span{1} = ker(A ). For the general case k ij ≥ 0, this property cannot be guaranteed and needs to be assumed. This explains Assumption (A3). Assumption (A4) is needed to derive the entropy production inequality (10). It is satisfied for reversible reactions, see [9,Lemma 6] or [10]. It also holds in the context of a tumor-growth model, see [25].
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 where K > 0 only depends on the given data. This is an estimated version of (10). The term involving ε is needed to conclude a uniform L 2 estimate for w k , which is sufficient to apply the Leray-Schauder fixed-point theorem in the finite-dimensional Galerkin space. The ε-independent gradient estimate for x k i cannot be used since it does not give an estimate for w k i (see (7)). It is possible to analyze system (12)-(13) for ε = 0-see Step 2 of the proof of Theorem 3-but we lose the information about w k and obtain a solution in terms of ρ k . The term involving ε is technical and not essential for the numerical simulations (or the structure preservation). However, we are not able to prove an existence result in terms of the entropy variable without such a regularization.
Remark 2 (Conservation of partial mass) When r i = 0, we have from (1) conservation of the partial mass ρ i L 1 ( ) . This conservation property does not hold exactly on the discrete level because of the ε-regularization. It holds that for any δ > 0, there exists ε 0 > 0 such that for any 0 < ε < ε 0 (ε is the value in (12)), The proof is the same as in [26,Theorem 4.1]. As δ > 0 can be chosen arbitrarily small, this shows that the numerical scheme preserves the partial mass approximately.

Theorem 3 (Convergence of the Galerkin solution) Let Assumptions (A1)-(A4) hold.
Let (ρ k , k ) be a solution to (12)-(13) and set for y ∈ , t ∈ ((k − 1)τ, kτ ], i = 1, . . . , n and introduce the shift operator (σ τ ρ τ i )(y, t) = ρ k−1 i (y) for y ∈ and t ∈ ((k − 1)τ, kτ ]. Then, there exist subsequences (not relabeled) such that, in the subsequent limits ε → 0, then N → ∞, and finally τ → 0, where In Theorem 3, ·, · denotes the duality bracket between H 1 ( ; R n−1 ) and H 1 ( ; R n−1 ). The difficult part of the proof is the estimate of the diffusion term because of the contribution of the electric field. We show in Lemma 8 that holds for some constants K, K 1 , K 2 > 0, which are independent of ε, N, and τ . Then, the uniform L ∞ bound for x k i gives a uniform H 1 ( ) bound for x k i and consequently for ρ k i . Weak compactness allows us to pass to the limits ε → 0 and N → ∞, and the limit τ → 0 is performed by means of the Aubin-Lions lemma.
The paper is organized as follows. In Section 2, we detail the thermodynamic modeling of system (1)-(3). Some auxiliary results on the formulation of the fluxes J i and the inversion of the map ρ → w are presented in Section 3. Sections 4 and 5 are devoted to the proof of the main theorems. Finally, some numerical experiments are shown in Section 6.

Modeling
We consider an isothermal electrolytic mixture of n fluid components in the bounded domain ⊂ 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 described by the partial mass densities ρ 1 , . . . , ρ n and the electric field E. The partial mass density represents the mass of the substance per unit volume. In thermodynamics, one introduces also the molar concentrations c i , signifying the amount of substance per unit volume. These quantities are related by ρ i = M i c i , where M i is the molar mass, the mass of the substance, divided by its amount. The total concentration is defined by c tot = n i=1 c i . Furthermore, x i = c i /c tot = ρ i /(c tot M i ) denotes the molar fraction, being the amount of the substance, divided by the total amount of all constituents of the mixture. We suppose the quasi-static approximation E = −∇ , where is the electric potential.
The evolution of the mass densities ρ i is governed by the partial mass balances [13, (4)] is the vector of molar fractions, J i the diffusion flux, and r i (x) the mass production rate of the ith species. We assume that the total flux and the total production vanish, which are necessary constraints to achieve total mass conservation, ∂ t n i=1 ρ i = 0. We suppose that the total initial mass is constant in space, n i=1 ρ 0 i = ρ tot > 0, which implies that the total mass is constant in space and time, n i=1 ρ i (t) = ρ tot for t > 0.
The electric potential is given by the Poisson equation [14, (3) and (25)] where ε 0 is the dielectric constant, χ the dielectric susceptibility, F the Faraday constant, z i the charge number of the ith species, and f (y) with y ∈ models the charge of fixed background ions. The basic assumption of the Maxwell-Stefan theory is that the difference in speed and molar fractions leads to a diffusion flux. They are implicitly given by the driving forces d i according to [2, (200 where the numbers D ij = D ji are the Maxwell-Stefan diffusivities. Inserting the definition In the present situation, the driving force is given by two components, the variation of the chemical potential μ i and the contribution of the body forces b i [2, (211)]: where R is the gas constant and T the (constant) temperature. Since (D ij ) is symmetric, summing (18) We assume that the only force is due to the electric field (i.e., we neglect effects of gravity It remains to determine the chemical potential. We define it by is the mixing free energy density [13, (23)]. Then, and the driving force becomes where z = (z 1 , . . . , z n ) and x = (x 1 , . . . , x n ). The Gibbs-Duhem equation shows that the pressure vanishes, which is consistent with our choice of the driving force (see [2, (211)]). The driving force in [35, (7)] contains a non-vanishing pressure that is related to our expression for the total body force. The resulting driving force (19), however, is the same. We summarize the model equations: and the relations In particular, we scale the particle densities by ρ tot (then the scaled quantities satisfy n i=1 ρ i = 1) and the electric potential by F /(RT ).

Auxiliary results
We collect some auxiliary results needed for the existence analysis. The starting point is relation (2) between the fluxes J i and the gradients ∇x i . Observe that the coefficients k ij depend on ρ i via c tot = n i=1 ρ i /M i . This dependency does not complicate the analysis since the results in Section 3 hold pointwise for any given ρ i and c tot is uniformly bounded by

Lemma 4 (Formulations of J i ) Eqs. 23 can be written equivalently as
The last expression for J i shows that the partial mass balances (1) can be formulated as where ρ = ρ(w) and B = B(ρ(w)). By Definition (7), w is a function of ρ (and ). The inverse relation ρ(w) is discussed in the following subsection.

Inversion of ρ → w
Definition (7) defines, for given ∈ R, a mapping x → w. We claim that this mapping can be inverted. If the molar masses are all the same, M := M i , this can be done explicitly: and ρ n = 1 − n−1 i=1 ρ i . Unfortunately, when the molar masses are different, we cannot derive an explicit formula. Instead, we adapt Lemma 6 in [8]. (w 1 (x), . . . , w n−1 (x)), where

Proof of theorem 1
Step 1: existence of solutions. The idea is to apply the Leray-Schauder fixed-point theorem. We need to define the fixed-point operator. For this, let χ ∈ L ∞ ( ; R n−1 ) and σ ∈ [0, 1]. Since (y, ) → ρ i (ω(y), ) is a bounded function with values in (0, 1), we can use Schauder's fixed-point theorem and standard arguments to show the existence of a solution k − D ∈ P N of the nonlinear finite-dimensional problem for all θ ∈ P N . In particular, the solution is unique and we have k ∈ L ∞ ( ), for → ρ i (ω, ) is Lipschitz continuous and D ∈ L ∞ ( ). Next, we wish to solve the linear finite-dimensional problem where a(u, φ) = ∇φ : B(χ + w D , k )∇udy + ε u · φdy, for u, φ ∈ V N , where we set in case k = 1, ρ (u 0 + ω D ) := (ρ 0 ) . Since χ + w D ∈ L ∞ ( ; R n−1 ) 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 are continuous on V N . Exploiting the equivalence of the norms in the finite-dimensional space V N , we find that for some constant K N > 0, which implies that a is coercive on V N . By the Lax-Milgram lemma, there exists a unique solution u ∈ V N ⊂ L ∞ ( ; R n−1 ) to (29) satisfying for some constant K F , which is independent of τ and σ . Using again the bounds for ρ and r , we find that the constant K F is independent of k . Since all norms are equivalent in the finite-dimensional setting, this provides a uniform L ∞ ( ) bound for u.
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 ): We claim that the first term on the left-hand side is greater than the difference of the entropies at time steps k and k − 1. To show this, we split the entropy density into two parts, where we recall that Therefore, using ρ k n − ρ k−1 For the estimate of h 2 , we first observe that We infer from the Poisson Eq. 13 and Young's inequality that Taking into account the property r n (ρ k ) = − n−1 i=1 r i (ρ k ), definition (7) of w k i , and Assumption (A4), we compute Combining (31)-(33) gives the conclusion.
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 main task is to estimate the diffusion part.
Lemma 8 (Estimate of the diffusion part) There exist constants K 1 > 0 and K 2 > 0, both independent of ε, τ , and N, such that Proof We drop the superindex k in the proof to simplify the notation. Recall that A = A| im(A) , where im(A) = span{1} ⊥ . We introduce as in the proof of Lemma 12 in [8] the symmetrization A S = P −1/2 AP 1/2 , where P 1/2 = M 1/2 X 1/2 and M 1/2 := diag(M  . Then, A −1 S = P −1/2 A −1 P 1/2 is a self-adjoint endomorphism whose smallest eigenvalue is bounded from below by some positive constant which depends only on (k ij ).
In view of n i=1 (B∇w) i = 0, it follows that Adding this expression to (34), we find that The matrix A −1 S is positive definite on im( A S ) = span{ρ 1/2 } ⊥ . A simple computation shows that the vector M −1/2 (2∇x where K 1 > 0 and K 2 > 0 depend on M 1 , . . . , M n . Since x i and ρ i x −1/2 i = (ρ i c tot M i ) 1/2 are bounded, the previous inequality becomes where K 3 depends on K 2 and z i . In the following, let K > 0 be a generic constant independent of ε, n, and τ . We estimate the expression involving the boundary term We infer from (35) and (36) that By the boundedness of c i , the elliptic estimate for the Poisson equation gives This proves the lemma.
Combining the discrete entropy inequality (15) and the estimate of Lemma 8 and summation over k leads to the following result.
is uniformly bounded in H 1 ( ), proving the claim. Observing that the embedding H 1 ( ) → L 2 ( ) is compact, there exist subsequences, which are not relabeled, such that as ε → 0, 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 := n i=1 ρ i /M i strongly in L 2 ( ), and we can identify ρ i = c tot M i x i for i = 1, . . . , n. Furthermore, weakly in L q ( ) for any q < 2 and i = 1, . . . , n. Since (D ε i ) is bounded in L 2 ( ), there exists a subsequence which converges to some function D i weakly in L 2 ( ). By the uniqueness of the weak limits, we can identify D i = D i . This shows that the convergence (41) holds in L 2 ( ). 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 2 ( ).

Numerical experiments
In this section, some numerical experiments based on scheme (12)- (13) in one space dimension are presented. We stress the fact that the experiments just serve as a feasibility study and more effort is necessary to perform two-or three-dimensional simulations showing, for instance, the effects coming from the mixed boundary conditions. In the context of semiconductor simulations, we refer to [16]. The onedimensional setting presented here models liquid electrolytes, which can be used as a simplified version of a model for dye-sensitized solar cells [35].
For the first example, we suppose that all three molar masses are equal to one and that the boundary conditions for the electric potential are in equilibrium, i.e., M 1 = M 2 = M 3 = 1 and (y) = 0 for y ∈ {0, 1}. The dynamics of the particle densities and the electric potential are shown in Fig. 1. The solution at time t = 17 is essentially stationary and, in fact, in equilibrium. Because of the choice of the parameters, the stationary solution is symmetric around x = 1 2 . In order to study the convergence of our iterative linearization procedure, we plot the evolution of the iteration parameter m over time in Fig. 2. It turns out that the  The situation changes drastically when the molar masses are different (example 2). Figure 3 shows the stationary solutions with the same parameters as in the previous example except M 1 = 6. Here, the discrete relative entropy is defined by where (ρ k h , k h ) is the finite-element solution at time kτ and (x ∞ h , ∞ h ) is the stationary solution. The integral and gradients are computed by the trapezoidal and gradient routines of MATLAB. The semi-logarithmic plot of the relative entropy shows that the entropy converges to zero exponentially fast.
For example 3, we choose the same initial conditions and parameters as before, but we take nonequilibrium boundary data (0) = 10, (1) = 0. The solutions at time t = 8 for various molar masses M 1 are displayed in Fig. 4. Since ρ 1 and ρ 2 have  In example 4, we interchange the roles of M 1 and M 2 , i.e., we choose M 1 = 1 and M 2 ∈ {2, 4, 6}. We observe in Fig. 5 that the first species is more concentrated at the right boundary while in the previous example, this holds true for the second species.
The previous examples show that the convergence rate to equilibrium strongly depends on 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. This behavior can be observed in Fig. 6 (example 5), where we have taken the same parameters as in the previous example but neglect the electric field. In this situation, the steady state is constant in space and explicitly computable; indeed, we have ρ ∞ i = mean( ) −1 ρ 0 i L 1 ( ) . 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 for the situation of example 3 (nonequilibrium boundary conditions for the potential). We choose the time t = 0.01 and the time step size τ = 10 −4 . The solutions are computed on nested meshes with grid sizes h ∈ {0.01, 0.005, 0.0025, 0.0006, 0.0001} and compared to the reference solution, computed on a very fine mesh with 25601 elements (h ≈ 4 · 10 −5 ). As expected, we observe a second-order convergence in space; see Fig. 7.