MATHEMATICAL MODELING OF A DISCONTINUOUS SOLUTION OF THE GENERALIZED POISSON–NERNST–PLANCK PROBLEM IN A TWO-PHASE MEDIUM

. In this paper a mathematical model generalizing Poisson–Nernst– Planck system is considered. The generalized model presents electrokinetics of species in a two-phase medium consisted of solid particles and a pore space. The governing relations describe cross-diﬀusion of the charged species together with the overall electrostatic potential. At the interface between the pore and the solid phases nonlinear electro-chemical reactions are taken into account provided by jumps of ﬁeld variables. The main advantage of the generalized model is that the total mass balance is kept within our setting. As the result of the variational approach, well-posedness properties of a discontinuous solution of the problem are demonstrated and supported by the energy and entropy estimates.

1. Introduction. The paper studies mathematical aspects of a generalization of the Poisson-Nernst-Planck (PNP) system with respect to coupling phenomena, the total mass balance, positivity of species concentrations, as well as nonlinear interface reactions and discontinuous field functions.
The PNP models arise in applications related to reaction-diffusion phenomena in biology and electrochemistry, see, for instance, [10,13,24,29]. In particular, we are motivated by the application to an electrolyte solution in Lithium ion batteries, see [2,23]. Because of drawbacks of the classic PNP model, which, for example, violates the mass balance and coupling phenomena, its improvement was suggested in [8,9,14]. The improvement is based on general thermodynamic principles, we refer the interested reader to [6,16,25].
The main features of our approach are given below. We examine the cross-diffusion of multiple charged species which is expressed by mass concentrations and an overall electrostatic potential. In the Fick's law, the corresponding diffusion fluxes are formulated in terms of electro-chemical potentials and assumed to be coupled through diffusivity matrices, whereas the coupling is usually violated in the classic case.
The multi-component solution is considered in a two-phase domain composed of two disjoint parts. One part is the solid phase consisted of particles, which are included in a surrounding pore space. We perceive the differences between equations in the solid and in the pore phases. In the solid phase, in other words inside the particles, we set the Gauss's flux law for the electrostatic potential and we define usual electro-chemical potentials. In the pore phase, the electrostatic potential and the mass concentrations are coupled within the Gauss's flux law. The corresponding quasi-Fermi electro-chemical potentials are derived from a Landau thermodynamic potential (see (5)) and depend on the pressure parameter.
The pressure parameter is determined from the force balance. It is a consequence of the Navier-Stokes equations with the flow velocity assumed to be zero. For physical consistency, the mass concentrations should be within a Gibbs simplex which implies that the total mass balance and positivity constraints hold. In its turn, accounting for the pressure parameter and for the total mass balance (see Theorem 2.1) avoids the main thermodynamic drawback of the classic model.
We aim to describe electro-chemical reactions at the phase interface. For this purpose, the field variables should allow a jump, while the diffusion fluxes and the electric current are assumed to be continuous in this paper. For the modeling of discontinuous fluxes, see [11,12]. The interface reactions lead to inhomogeneous boundary conditions compared to the standard PNP models. Moreover, we show that the boundary fluxes of species cannot be nontrivial constants or linear functions but they should depend nonlinearly on the mass concentrations.
From a mathematical point of view, the existence theorems providing a weak solvability to the PNP problem are based on the variational theory and are given in [17,26,27], and in [4] in the stationary case. We refer to [1,15,28] for homogenization, and to [5,14] for numerical methods suitable for reaction-diffusion models.
In this paper we formulate the generalized PNP problem considering all specialties mentioned above. The main difficulty for its analysis is connected with the strong nonlinearity of the governing PDEs as well as the nonlinearity of the boundary conditions. To describe the two-phase medium, discontinuous solutions allowing jumps through the phase interface are employed. We give a weak formulation of the problem and describe properties of its solution: well-posedness, the total mass balance, the weak maximum principle, and dynamic stability in the sense of Lyapunov, which are supported by the energy and entropy estimates.
In the well-posedness theorems, the matrices of diffusion are assumed to satisfy the so-called weak assumption, see (16). It guarantees the total mass balance and a local existence of the positive solution at least in a small time interval. To ensure global non-negativity of the solution we suppose the stronger assumption, see (17), implying that the diffusion fluxes are decoupled. By this, the nonlinear functions of interface reactions need to satisfy sufficient assumptions on the mass balance, the growth, the continuity, and the positive production rate.
To prove existence of a weak solution, we derive a reduced formulation excluding the constraints on the mass concentrations and apply the Schauder-Tikhonov fixed point theorem as suggested in [19] based on the energy estimate and following [22]. Afterwards, the solution of the reduced problem is examined for properties of the total mass balance and positivity. We show under which conditions the reduced and the original formulations are equivalent. Under an additional assumption of regularity of the function of the electrostatic potential and the function of boundary fluxes, the solution is unique. In the last part we examine the entropy dissipation and prove the entropy estimate following [7].
To begin with, we describe geometry of the problem.
Let Ω be a bounded domain in R d for d = {1, 2, 3} with the Lipschitz boundary ∂Ω. Motivated by a two-phase medium, we split the domain Ω in two non-trivial disjoint parts ω and Q separated by an interface ∂ω. The domain ω called the solid phase represents multiple disjoint particles surrounded by the pore space Q. The thickness of the interface ∂ω is assumed to be small in comparison with the size of particles and we do not take it into account. Let the interface ∂ω be a Lipschitz continuous manifold with the unit normal vector denoted by ν = (ν 1 , . . . , ν d ). We set ν pointed outward to each particle in ω, thus inward to the pore part Q. An example domain Ω = Q ∪ ω ∪ ∂ω is illustrated in Fig. 1.
For each simply connected compact set in ω, which we call a particle, we distinguish the negative face ∂ω − as the boundary of ω, and the positive face ∂ω + as its opposite part corresponding to the boundary of the pore phase Q, see zoom in  the Fick's law of diffusion: where D ij are given d-by-d matrices of diffusion for each indexes i, j = 1, . . . , n.
Here ∇ stands for the gradient vector of the length d, and the upper sign swaps between rows and columns.
In the solid phase (0, T ) × ω, the electrostatic potential ϕ is represented by the Gauss's flux law with the help of a d-by-d matrix A of permittivity: the Gauss's flux law : −div(∇ϕ A) = 0; (2a) while the electro-chemical potentials µ i appeared in the constitutive law (1b) are given by the usual expressions (where i = 1 . . . , n): Here the Boltzmann constant k B , the temperature Θ, and volume factors of species β i are positive physical parameters.
In the pore phase (0, T ) × Q, equations for the electrostatic potential ϕ, electrochemical potentials µ i , i = 1, . . . , n, and the pressure parameter p are coupled as follows: quasi-Fermi electro-chemical potentials: The form of (3c) will be argued from the Landau thermodynamic potential (5) below.
The force balance (3b) came from the Navier-Stokes equations with the zero flow velocity as shown in [8,9]. The pressure parameter p from (3b) enters also the quasi-Fermi electro-chemical potentials in the form of (3c). As the consequence, further it will be possible to keep the total mass balance proved in Theorem 2.1 at the end of this section. This generalization has an advantage over the classic PNP equations.
By the reason of physical consistency we assume the balance of the total mass and positivity for the mass concentrations, which form the so-called Gibbs simplex: the total mass balance: Physically, (4a) follows that ∂ ∂t ( n i=1 ρ i ) = 0 implying the total mass balance when the flow velocity is zero.
The thermodynamic Landau grand potential of the system has the form (see [3]): After variation over admissible variables ρ, ϕ and p this provides two sets of relations in the solid phase ω and in the pore phase Q. The variation of L with respect to ϕ: ∂L ∂ϕ = 0, after integration by parts, follows the Poisson equations (2a) and (3a) supported by the inhomogeneous boundary conditions (10b) following later; the variation with respect to ρ i : ∂L ∂ρi = 0 gives the expressions (2b) and (3c) for µ i ; and the variation with respect to p: ∂L ∂p = 0 implies the total mass balance (4a) in Q. The variables µ i appear in L as Lagrange multipliers to the positivity constraint (4b).
The proper function spaces of the field variables in relations (1)-(5) will be provided further in (25).
The initial conditions in Q ∪ ω for the time derivative in (1a) are: The initial data are such that the conditions in the manner of (4) hold in Q ∪ ω, namely: the total mass balance: At the external boundary (0, T ) × ∂Ω associated with a bath, we set the usual Dirichlet conditions: (8b) We suggest that the Dirichlet boundary data satisfy the following relations for each index i = 1, . . . , n: the total mass balance: and compatibility conditions: At the interface (0, T ) × ∂ω we allow a jump of the mass concentrations and the electrostatic potential, while we assume that the diffusion fluxes and the electric current are continuous and inhomogeneous such that: (10b) The minus sign in front of ν is due to the normal direction pointed inside Q.
The first condition in (10b) implies the continuity of the electric current, while the second condition describes the charge exchange between the two phases. We emphasize that in (10a) functions g i (ρ,φ) depend nonlinearly on the variableρ, see an example in (13). Compared to the trivial case g i (ρ,φ) ≡ 0, specifying a non-zero flux g i (ρ,φ) at any of the two interface faces implies, generally, a boundary reaction, for example, a phase-boundary catalysis. More discussion of the interface conditions (10) is available in [11,12]. In particular, discontinuous fluxes are modeled in [11].

2.2.
Properties of the problem. For the well-posedness analysis of the generalized PNP problem we formulate next the assumptions on the data.
We suggest that the nonlinear functions g i (ρ,φ) in (10a) satisfy the following properties for every species i = 1, . . . , n: the growth conditions with the uniform bounds γ i the mass balance: the positive production rate : and the Lipschitz continuity with a constant K L > 0 : for all ρ and r = (r 1 , . . . , r n ) satisfying n k=1 ρ k = n k=1 r k = C on ∂ω ± , wherer = (r| ∂ω + , r| ∂ω − ). In (11c), the partition in the positive and the negative parts is determined by (12b) We note that a nontrivial constant g i cannot fulfill the zero product in (11c), and a linear function g i (ρ) cannot be uniformly bounded as in (11a). To give an illustrative example, the nonlinear functions where G j (ρ) := ρ j n k=1 ρ k such that |G j (ρ)| 1 and G j (ρ + )ρ − j = 0, fulfill all the conditions (11) with γ i 2 = 0 and γ i 1 = |∂ω||h i | 2 when the numbers h i ∈ R are chosen such that n i=1 h i = 0. In particular, g 2 (ρ,φ) = −g 1 (ρ,φ), and g i (ρ,φ) = 0 for other i = 3, . . . , n, are suitable. We note that the sign of h i in (13) defines the sign of the boundary flux J i ν according to the boundary condition (10a). This means that every species can flow either only into or only outside of the solid phase which is determined by the underlying electro-chemical reaction.
The electric permittivity matrix A and diffusitity matrices D and D ij are required to satisfy the following assumptions below.

VICTOR A. KOVTUNENKO AND ANNA V. ZUBKOVA
Let A ∈ R d×d be a symmetric and positive definite matrix, i.e. there exist positive numbers 0 < a ā such that The following properties of the diffusivity matrices are related to the respective constraints (4a) and (4b) on the mass concentrations. We assume that either the weak assumption: or the strong assumption: holds, where δ ij is the Kronecker delta such that δ ij = 1 for i = j and zero otherwise. Further we assume that the following strong ellipticity condition holds for m i D ij ∈ R d×d , i.e. there exist positive numbers 0 < d d such that and D ∈ R d×d in (16) and (17) is a symmetric and positive definite matrix implying existence of 0 < d 1 d 1 which satisfy the following inequalities: We note that the assumption (16) follows straightforwardly from (17), and the bounds in (18) and (19) coincide such that d 1 = d andd 1 =d when (17)  Theorem 2.1 (Balance of total mass). If the constraint (4a) and the assumption on the initial data (7a) hold, then the total mass n i=1 ρ i = n i=1 ρ in i = C is constant over time, thus implies the conservation ∂ ∂t ( n i=1 ρ i ) = 0. Together with the weak assumption on the diffusivity matrices (16) this is equivalent to the flux balance: . Conversely, to show the flux balance, first, we substitute the electro-chemical potentials µ i from (2b) in the constitutive law (1b) in the solid phase, and we put the quasi-Fermi electro-chemical potentials µ i from (3c) in (1b) in the pore phase. The gradient of the pressure parameter is taken from the force balance (3b). As the result we get the following expressions for the diffusion fluxes J i in the two phases, respectively: Second, equations (20) for J i are summed over i = 1, . . . , n. Due to the weak assumption on the diffusivity matrices (16) and the total mass balance (4a) we derive in the solid phase (0, T ) × ω that and similarly we get in the pore phase (0, T ) × Q: The proof is completed.
Using the assumptions of this section, next we formulate well-posedness results.
3. Well-posedness analysis. Existence and uniqueness of the weak solution of the PNP problem were proved rigorously in [20,21]. Therefore, here we formulate the main results omitting proofs. We start with eliminating the entropy variables p and µ from the problem (1)-(3).
Here in (26) we have used the notation 1 Q for the indicator function of the domain Q such that 1 Q = 1 in Q and zero otherwise (that is 1 Q = 0 in ω), and the notations Υ j and Υ for short: The system (25)- (27) can be extended with identities (2b), (3b), and (3c) to describe entropy variables: the pressure parameter p in Q and the electro-chemical potentials µ 1 , . . . , µ n in Q and ω. Indeed, after solving the variational equations (26) for ρ and ϕ, the pressure parameter p in the pore space Q can be recovered from the force balance (3b). Taking the divergence on both sides of (3b) leads to the Poisson equation −∆p = div Υ(ρ)∇ϕ which can be treated by standard tools.
Avoiding the constraints (4) in the system, we will simplify the problem for analysis.
The following theorem is based on the proof given in [20,21].  (6) conditions. The solution satisfies the a-priori estimates The replacement of this expression in the variational equation (26b) implies (31b). Analogously, for Υ j andῩ j it holds: Consequently, replacing Υ j (ρ) in (26a) byῩ j (ρ + ) implies (31a). In return, (31) follows (26). Indeed, according to the conclusion of Theorem 3.1, if the conditions ρ i > 0 and n i=1 ρ i = C are satisfied, then the equalities (32) and (33) hold again, thus proving the equivalence which results in the following corollary. 2) If additionally the strong assumption on diffusivity matrices (17) holds, then this ensures that the solution is non-negative globally in (0, T ) × (Q ∪ ω) for an arbitrary final time T > 0.
Uniqueness of the solution is guaranteed under specific assumptions as follows.
and let the nonlinear functionsρ i → g i ( · ,φ) of boundary fluxes be injective for all i = 1, . . . , n and satisfy the following assumption: there exists a number M 0, such that it holds the estimate for all ϕ (1) , ϕ (2) and ρ (1) In this case, the weak solution (25) of the problem (26) is assured to be unique.
If g i is not injective with respect to ρ i , then we can give an example of possible non-uniqueness of the solution in the particles. Indeed, assuming g i ≡ 0, we have that ρ 3.2. Entropy and its dissipation. In this section we examine the entropy-production rate which describes a dynamic stability of the system in the sense of Lyapunov with respect to the increasing time t → ∞.
First, we rewrite the Landau grand potential L in (5) by avoiding a constant value. Indeed, according to the total mass balance (4a) we have Therefore, for the mass concentrations ρ satisfying the total mass balance and positivity (4), from (5) we define the entropy S = − ∂L ∂Θ in the pore phase Q as follows: We note that this function is non-positive, since ξ ln ξ ξ − 1 holds for all ξ > 0, then if we assume in the total mass balance (4a) the density C = n i=1 1/β i . Second, we introduce the function of dissipation D expressing the entropy-production rate. Based on the formula for S and using n i=1 ∂ρ i /∂t = ∂C/∂t = 0 we define: According to the second law of thermodynamics the entropy-production rate should be non-negative. With the help of the following assumptions: scalar diffusivity matrices: scalar permittivity matrix: boundary charge neutrality: constant boundary concentrations: where I stands for the d-by-d identity matrix, below we calculate the function of dissipation D.
Theorem 3.3 (Entropy-production rate). Under the assumptions on data (37), for the mass concentrations ρ i satisfying the total mass balance and positivity (4), the entropy dissipation in the pore phase Q defined in (36) can be expressed equivalently as follows: Here, D 1 0, and the dissipation inequality D 0 can be assured by D 2 such that D 2 −D 1 . In particular, the entropy principle holds for α, g, and g i (ρ,φ) sufficiently small such that D 1 + D 2 0.
Proof. For a fixed time t ∈ (0, T ), we choose the test functionρ i = ln(β i ρ i ) in Q andρ i = 0 in ω and insert it in (26a), since ln(β i ρ D i ) = ln 1 = 0 at the boundary ∂Ω due to the assumption (37d). As the result we get: where the assumption (37a) on the diffusivity matrices was used. Summing of (39) over i = 1, . . . , n, inserting the expressions (27) for Υ i and using the expression of dS/dt from (36), if follows that The last term in the left-hand side of (40) is zero due to the following identity n i=1 ∇ρ i = ∇ n i=1 ρ i = ∇C = 0. Now we substitute in (26b) the test functionφ = 0 in ω andφ = n i=1 z i ρ i in Q, which is zero at ∂Ω due to the charge neutrality (37c), and obtain the equality where the assumption (37b) on the permittivity matrix was used. Equation (41) is multiplied by the constant factor −d/(aN A ) and the result is added to equation (40). After multiplication with the constant factor k B N A and using the identity |∇ρ i | 2 /ρ i = 2|∇( √ ρ i )| 2 we arrive at the formula (38).
Since D 1 0 and D 2 has no definite sign, either positive or small D 2 suffices the dissipation inequality D = D 1 + D 2 0.
We note that a similar to (38) expression can be guaranteed also for the entropyproduction rate in the solid phase under the charge neutrality assumption such that n i=1 z i ρ i = 0 in ω. Discussion. In physical applications, the generalized Poisson-Nernst-Planck system describes electro-chemical phenomena at the micro level. In the future work we will derive rigorously the averaged model when the size of solid pores is set to be small within the homogenization procedure. The result of the current paper provides the rigorous mathematical formulation supported by the uniform a-priori estimates for the inhomogeneous problem.