Publication Cover
Applicable Analysis
An International Journal
Volume 100, 2021 - Issue 2
1,190
Views
5
CrossRef citations to date
0
Altmetric
Articles

Homogenization of the generalized Poisson–Nernst–Planck problem in a two-phase medium: correctors and estimates

ORCID Icon &
Pages 253-274 | Received 11 Apr 2018, Accepted 22 Mar 2019, Published online: 07 Apr 2019

Abstract

The paper provides a rigorous homogenization of the Poisson–Nernst–Planck problem stated in an inhomogeneous domain composed of two, solid and pore, phases. The generalized PNP model is constituted of the Fickian cross-diffusion law coupled with electrostatic and quasi-Fermi electrochemical potentials, and Darcy's flow model. At the interface between two phases inhomogeneous boundary conditions describing electrochemical reactions are considered. The resulting doubly non-linear problem admits discontinuous solutions caused by jumps of field variables. Using an averaged problem and first-order asymptotic correctors, the homogenization procedure gives us an asymptotic expansion of the solution which is justified by residual error estimates.

COMMUNICATED BY:

AMS Subject Classifications:

1. Introduction

The paper is devoted to the mathematical study of homogenization of a non-linear diffusion model in a two-phase domain.

The Poisson–Nernst–Planck (PNP) model extends the diffusion law due to electro-kinetic phenomena. Namely, we consider cross-diffusion of multiple charged species coupled with an overall electrostatic potential. Motivated by the physical nature, species concentrations satisfy the total mass balance and the positivity conditions. Following [Citation1–4], this approach generalizes the classic PNP model.

The problem under consideration is characterized by the following issues.

We describe a two-phase medium with a micro-structure consisting of solid and pore phases which are separated by a thin interface. The corresponding geometry is represented by a disconnected domain. Therefore, field variables defined in the two-phase domain allow discontinuity with jumps across the interface.

A special interest of our consideration is the interface between the two phases because of electrochemical reactions that occur here. At the interface we state mixed, inhomogeneous Neumann and Robin-type conditions. Diffusion fluxes and the electric current are assumed continuous across the phase interface. The key issue is that the inhomogeneous boundary fluxes are to be described by non-linear functions of the field variables.

From a mathematical point of view, we examine a mixed system of partial differential equations of the parabolic-elliptic type. The governing equations are non-linear, coupled, and differ on the two phases. The non-linearity is due to the presence of electrochemical potentials in the model. The solvability of classic PNP systems was studied in [Citation5,Citation6]. Based on a general approach from [Citation7,Citation8], in the previous works [Citation9–11], we proved existence theorem for the generalized PNP problem and derived a-priori estimates.

Homogenization of diffusion equations is widely studied in the literature, see, for instance [Citation12–17] for adopted approaches. Most of the asymptotic results concern either linear equations, or homogeneous Neumann conditions excluding interface reactions, which are of primary importance in electro-chemistry. For possible transmission conditions stated at the interface we refer to [Citation18–20]. Homogenization of classic PNP equations was studied in [Citation21–23]. A homogenization procedure in a two-phase domain for steady-state Poisson–Boltzmann equations and homogeneous Neumann boundary conditions was investigated in [Citation24]. In the present work we continue this approach to the inhomogeneous conditions in the dynamic case. We rely on hydrostatic setting of the non-stationary problem, which is typical, e.g. for modelling of Li-ion batteries [Citation25]. For homogenization accounting for velocity fields, we refer to [Citation26,Citation27].

The difficulty of the homogenization procedure is caused by the two-phase domain. Typically, homogenization problems are considered in a perforated domain. In contrast, we describe a discontinuous prolongation from the perforated domain inside solid particles following the approach of [Citation28]. In this respect, the two-phase homogenization procedure differs from a perforated domain case. To describe jumps of the field variables across the interface and interface reaction terms, we will specify their suitable asymptotic orders.

To derive an averaged model, typically, the two-scale convergence is applied. As an advantage, we endow our asymptotic expansion with residual error estimates.

As the result of homogenization of the PNP model, we obtain an averaged model consisting of linear parabolic-elliptic equations and supported by first-order correctors. The correctors appear due to oscillating and interface data expressed by solutions of auxiliary cell problems in a unit cell. Respectively, there are three correctors given with respect to:

  • the periodic boundary function of the electric current at the phase interface;

  • the periodic matrix of permittivity;

  • the periodic matrices of diffusivity.

In order to justify cell problems we use the periodic unfolding technique. It is based on the unfolding operator and the averaging operator, which were defined for perforated domains in [Citation29]. We extend the concept of the unfolding operator to a two-phase domain, and we define its extension to a non-periodic boundary according to [Citation30].

The paper has the following structure. Section 2 contains a brief description of the unfolding method: definitions and main properties. In Section 3, we formulate the PNP problem and describe its solution. Section 4 accounts for auxiliary cell problems. In Section 5, a homogenization procedure is introduced and proved rigorously. By this, the averaged problem is formulated and supported by error estimates of the corrector terms.

2. Unfolding technique

Let Ω be a domain in Rd, where dN, with the smooth boundary ∂Ω and the unit normal vector ν, which is outward to Ω. We consider the unit cell Y=(0,1)d consisted of the isolated solid part ω¯Y and the complementary pore part Π:=Yω¯ such that Y=Πωω and ωY=. The interface ω is assumed to be a smooth continuous manifold with a unit normal vector ν. We set ν outward to ω, thus inward to Π.

For a small εR+ every spacial point xRd can be decomposed as follows (1) x=εxε+εxε(1) into the floor part x/εZd and the fractional part {x/ε}Y. There exists a bijection C:ZdN implying a natural ordering, and its inverse is C1:NZd. Based on (Equation1), we can determine a local cell Yεl with the index l=C(x/ε), such that xYεl, and {x/ε}Y are the local coordinates with respect to the cell Yεl.

Let Iε:={lN:YεlΩ} be the set of indexes of all periodic cells contained in Ω, and Ωε:=int(lIεYεl¯) be the union of these cells. For every index lIε, after rescaling y={x/ε}, the local coordinate yω determines the solid particle such that {x/ε}ωεl with the smooth boundary ωεl. Its complement composes the pore Πεl:=Yεlωεl¯ by analogy with Π=Yω¯.

Gathering over all local cells, we define the multi-component domain of periodic particles (the solid phase) denoted by ωε:=lIεωεl with the union of boundaries ωε:=lIεωεl and the unit normal vector ν to each of ωεl. The Hausdorff measure |ωε| of the interface ωε is of the order O(ε1) due to |ωεl|=O(εd1) and the cardinality |Iε|=O(εd). We denote Πε:=Ωεω¯ε, which is a perforated domain. Adding a thin layer ΩΩε, possibly attached to the external boundary ∂Ω, composes the pore phase Qε:=(ΩΩε)Πε.

For fixed ε>0, a two-phase medium associated to the disconnected domain Qεωε with the external boundary ∂Ω and the interface ωε is considered, see an example geometry in Figure .

Figure 1. A two-phase domain consisted of solid particles ωε and the pore space Qε with the phase interface ωε.

Figure 1. A two-phase domain consisted of solid particles ωε and the pore space Qε with the phase interface ∂ωε.

Following [Citation29,Citation30] and based on the decomposition (Equation1), we introduce two linear continuous operators: the unfolding operator f(x)Tε:H1(Qε)×H1(ωε)L2(Ω;H1(Π)×H1(ω)), defined by (2) (Tεf)(x,y)=fεxε+εy,a.e. for xΩεf(x),a.e. for xΩΩεandyΠω,(2) and its left-inverse operator u(x,y)Tε1:L2(Ω;H1(Π)×H1(ω))H1(lIεΠεl)×H1(ωε)×H1(ΩΩε) called the averaging operator: (3) (Tε1u)(x)=1|Y|Πωuεxε+εz,xεdz,a.e. for xΠεωε,1|Y|Πωu(x,y)dy,a.e. for xΩΩε,(3) where |Y| stands for the Hausdorff measure of the set Y in Rd. We note that Tε1u in (Equation3) is discontinuous across Yεl and Ωε. In the homogenization theory, usually x refers to as a macro-variable, y as a micro-variable, and (x,y) as the two-scale variables.

Lemma 2.1

Properties of the operators Tε and Tε1 in the domain

For arbitrary functions f,q,hH1(Qε)×H1(ωε), the following properties hold: (4a) (i)invertibility:(Tε1Tε)f(x)=f(x);(4a) (4b) (TεTε1u)(x,y)=u(y) when u is constant for xQεωε,or a periodic function u(y) of yΠω for xΠεωε;(ii)product rule:Tε(fq)=(Tεf)(Tεq);(4b) (4c) (iii)integration rules in the periodic domain and in the boundary layer:Πεωεf(x)q(x)dx=1|Y|ΩεΠω(Tεf)(x,y)(Tεq)(x,y)dydx,(4c) (4d) ΩΩεf(x)q(x)dx=1|Y|ΩΩεΠω(Tεf)(x,y)(Tεq)(x,y)dydx;(4d) (4e) (iv)boundedness of Tε in the L2-norm and the H1-semi-norm:Qεωεh2(x)dx=1|Y|ΩΠω(Tεh)2(x,y)dydx,(4e) (4f) Qεωε|h|2(x)dx=1ε2|Y|ΩΠω|y(Tεh)|2(x,y)dydx.(4f)

Proof.

(i) For xΩΩε and fL2(ΩΩε), we calculate straightforwardly (Tε1Tε)f(x)=(Tε1(Tεf))(x)=(Tε1f)(x)=f(x). For xΠεωε and fL2(Πε)×L2(ωε) according to (Equation1), the definitions (Equation2) and (Equation3) with zY we have (Tε1(Tεf))(x)=1|Y|Πωfεεxε+εzε+εxεdz=1|Y|Πωf(x)dz=f(x), since x/ε+z=x/ε, hence (Equation4a) holds. The assertion for TεTε1 can be checked.

(ii) The identity (Equation4b) is obvious.

(iii) The proof of (Equation4c) is known (see [Citation29, Section 2]). In the boundary layer, we derive straightforwardly (Equation4d) from (Equation2) and (Equation3).

(iv) Taking first q=f=h, then q=f=h in (Equation4c) and (Equation4d), summing them, and using Tε(f)=(1/ε)y(Tεf) due to the chain rule =(1/ε)y, we arrive at (Equation4e) and (Equation4f). This completes the proof.

A function fH1(Qε)×H1(ωε) given in the two-phase domain allows discontinuity across the interface ωε, see zoom in Figure . In each local cell Yεl we distinguish the negative face (ωεl) as the boundary of the particle ωεl, and the positive face (ωεl)+ as the opposite part of the boundary of the pore Πεl. Gathering over all local cells establishes the positive and negative faces of the interface as ωε±=lIε(ωεl)±. We set the interface jump of f across ωε by (5) [[f]]:=f|ωε+f|ωε,(5) where the corresponding traces of f at ωε± are well defined, see [Citation31, Section 1.4]. Analogously, we define the interface jump for a function u(y)H1(Π)×H1(ω) in the unit cell as [[u]]y:=u|ω+u|ω.

Motivated by the traces, we extend to the interface ωε the unfolding operator f(x)Tε:L2(ωε)L2(Ωε)×L2(ω) by (6) (Tεf)(x,y)=fεxε+εy,a.e. for xΩεandyω,(6) and similarly the averaging operator u(x,y)Tε1:L2(Ωε)×L2(ω)L2(ωε), (7) (Tε1u)(x)=1|Y|Πωuεxε+εz,xεdz,a.e.\ for\ xΩε.(7) Their properties are stated below in the manner of Lemma 2.1.

Lemma 2.2

Properties of the operators Tε and Tε1 at the interface

For arbitrary functions f,qL2(ωε), the following properties hold: (8a) (i)invertibility: (Tε1Tε)f=f;(8a) (8b) (ii)product rule:Tε(fq)=(Tεf)(Tεq);(8b) (8c) (iii)integration rule:ωεf(x)q(x)dSx=1ε|Y|Ωεω(Tεf)(x,y)(Tεq)(x,y)dSydx;(8c) (8d) (iv)boundedness of Tε in the L2-norm:ωεf2(x)dSx=1ε|Y|Ωεω(Tεf)2(x,y)dSydx.(8d)

Proof.

The proof of assertions (i) and (ii) is similar to the proof in Lemma 2.1. The proof of (Equation8c) is known (see [Citation29, Section 4]). Taking q=f in (Equation8c) immediately follows formula (Equation8d) in (iv).

The geometric construction of the operators Tε and Tε1 in this section will be used further for homogenization over Qεωε and ωε as ε0+.

3. Problem formulation

We formulate a generalized Poisson–Nernst–Planck system depending on a fixed parameter ε>0, see [Citation9–11]. We consider the number n of charged species with specific charges zi, molar masses mi>0, volume factors βi>0, and unknown concentrations ciε for i=1,,n and n2. By ϕε we denote the overall electrostatic potential. The two-phase medium introduced in Section 2 will be characterized below separately in the pore phase Qε and the solid phase ωε.

For the time-space variables (t,x) in (0,τ)×(Qεωε) with a fixed final time τ>0, we consider the following governing equations for species i=1,,n: (9a) The Fick's law of diffusion:ciεtdivJiε=0;(9a) (9b) cross-diffusion fluxes:(Jiε)=j=1ncjε(μjε+1QεεκηNACvε)mi(Tε1Dij);(9b) (9c)  electrochemical potentials:μiε=kBΘln(βiciε)+1QεεκNA1Cpε+ziϕε;(9c) (9d) the Darcy flow in pores:ηvε+pε=j=1nzjcjεϕε,div vε=0;(9d) (9e) and the Gauss's flux law:div((ϕε)(Tε1A))=1Qεj=1nzjcjε.(9e) The indicator function 1Qε is equal to 1 in Qε, and 0 in ωε. The Equation (Equation9c) contains the Boltzmann constant kB, the temperature Θ, the Avogadro constant NA, and κ1 in (Equation9c) allows us to average the non-linear diffusion fluxes (see (Equation71)). The fluxes contain the flow velocity following e.g. [Citation3,Citation4], and the dependence of potentials on the fluid pressure is due to the works by Dreyer (see [Citation1,Citation2]). The Equations (Equation9b)–(Equation9d) will be not solved with respect to electro-chemical potentials (μ1ε,,μnε), flow velocity vector vε=(v1ε,,vdε) with the drug coefficient η, and the pressure pε, but rather reduced within a weak formulation (see (Equation22)). Conversely, after finding (c1ε,,cnε) and ϕε, all the entropy variables (μ1ε,,μnε), vε, pε can be restored from the Equations (Equation9c) and (Equation9d) supported by suitable boundary conditions.

In (Equation9e) and (Equation9b) the d-by-d matrices A and Dij for i,j=1,,n imply the electric permittivity and diffusivity, respectively. They can be discontinuous in the two-phase unit cell Πω and satisfy the following assumptions.

  • A(y)Rd×d for yΠω is uniformly bounded and symmetric positive definite (spd) matrix: (10) there exist 0<a_a¯ such that a_|ξ|2ξA(y)ξa¯|ξ|2for ξRd;(10)

  • miDij(y)Rd×d for yΠω are uniformly bounded and elliptic matrices: there exist 0<d_d¯ such that d_i=1n|ξi|2i,j=1nξimiDij(y)ξjd¯i=1n|ξi|2for ξ1,,ξnRd;

  • the mass balance needs a symmetric positive definite (see (Equation13) below) matrix D~(y)Rd×d for yΠω such that: (11) i=1nmiDij(y)=D~(y)for j=1,,n.(11)

It is worth noting that conditions (Equation11) together with (Equation14a) below are sufficient to conserve the mass within the laws (Equation9b)–(Equation9d) as follows: i=1n(Jiε)=j=1ncjεμjε+1QεεκηNACvεTε1D~=kBΘj=1ncjε+1QεεκNA1Cj=1ncjεpε+j=1nzjcjεϕε+ηCj=1ncjεvεTε1D~=0. For homogenization reason, we assume that the diffusivity matrices Dij from (Equation11) admit the asymptotic decomposition as follows (12) miDij(y)=δijD(y)+εD~ij(y)for yΠω,(12) with d-by-d matrices D~ij, i,j=1,,n and a d-by-d uniformly bounded, symmetric positive definite matrix D such that (13) d_|ξ|2ξD(y)ξd¯|ξ|2for ξRd.(13) The oscillating matrices (Tε1Dij)(x)=Dij({x/ε}) and (Tε1A)(x)=A({x/ε}) in the Equations (Equation9b) and (Equation9e) are defined in Ω, and they are periodic in Ωε.

A constant C>0 in (Equation9c) stands for the summary concentration. For the physical consistency, species concentrations need to satisfy in pores (0,τ)×Qε: (14a) the total mass balance:i=1nciε=C;(14a) (14b) the positivity:ciε>0,for i=1,,n.(14b)

The system (Equation9) is supported by the initial condition for ciinH1(Ω): (15) ciε=ciinon Qεωε,(15) where the initial data satisfy the relations in the manner of (Equation14) in pores Qε: (16) i=1nciin=C,ciin>0,for i=1,,n.(16) For given functions ciDH1(0,τ;L2(Ω))L2(0,τ;H1(Ω)) and ϕDL(0,τ;H1(Ω)) the Dirichlet boundary conditions are: (17) ciε=ciD,for i=1,,n,ϕε=ϕDon (0,τ)×∂Ω,(17) with the boundary data satisfying the similar relations and compatibility: (18) i=1nciD=C,ciD>0on (0,τ)×∂Ω;ciD(0,)=ciinin Qεωε.(18) The most delicate part of modelling is the interface conditions on (0,τ)×ωε: (19a) [[Jiε]]ν=0,Jiεν=ε2gi(cˆε,ϕˆε);(19a) (19b) [[(ϕε)(Tε1A)]]ν=0,(ϕε)(Tε1A)ν+αε[[ϕε]]=Tε1g,(19b) where the jump across ωε is defined in (Equation5). The notation cˆε:=(cε|ωε+,cε|ωε) and ϕˆε:=(ϕε|ωε+,ϕε|ωε) implies the pair of traces at the phase interface ωε. The function gL(0,τ;L2(ω)) denotes the electric current through the interface in the unit cell, and (Tε1g)(x)=g({x/ε}) in (Equation19b) is periodic at ωε. The capacitance density α>0. The equality in (Equation19b) implies that the potential jump is asymptotically small [[ϕε]]=O(ε) in the electric double layer. The factor ε2 in (Equation19a) is used in Theorem 5.1 for averaging of the nonlinear, thus non-periodic interface data (see (Equation72)), and the factor 1/ε in (Equation19b) will be explained later in (Equation24). For modelling and numerical simulations of data for scaling of potentials, interface and boundary conditions, we refer to [Citation25].

In (Equation19a), the functions (cˆ,ϕˆ)gi, R2n×R2R, i=1,,n, describing the boundary fluxes of species with respect to the traces cˆ:=(c|ωε+,c|ωε) and ϕˆ:=(ϕ|ωε+,ϕ|ωε) of the variables c=(c1,,cn) and ϕ, should satisfy (20a) balance of the mass:i=1ngi(cˆ,ϕˆ)=0;(20a) (20b) positive production rate at ωε+:gi(cˆ,ϕˆ)min(0,ci|ωε+)=0;(20b) (20c) uniform boundedness (Kg>0):|gi(cˆ,ϕˆ)|2Kg.(20c) The example of gi satisfying all assumptions (Equation20) can be found in [Citation9,Citation10], e.g. g1(cˆ,ϕˆ)=max(0,c1|ωε+)max(0,c2|ωε+)[k=1nmax(0,ck|ωε+)]2,g2(cˆ,ϕˆ)=g1(cˆ,ϕˆ), and gk(cˆ,ϕˆ)=0 for k3.

A weak formulation of the generalized PNP problem is the following one: Find (c1ε,,cnε) and ϕε such that for i=1,,n: (21a) ciεL(0,τ;L2(Qε)×L2(ωε))L2(0,τ;H1(Qε)×H1(ωε)),(21a) (21b) ϕεL(0,τ;H1(Qε)×H1(ωε)),ciεϕiεL2((0,τ)×(Qεωε)),(21b) which satisfy the Dirichlet boundary conditions (Equation17), the initial conditions (Equation15), the total mass balance and positivity conditions (Equation14), and fulfil the equations: (22a) 0τQεωεciεtc¯i+j=1nkBΘ∇cjε+εκ1QεΥj(cε)ϕεmi(Tε1Dij)c¯idxdt=0τωεε2gi(cˆε,ϕˆε)[[c¯i]]dSxdt,i=1,,n,(22a) (22b) Qεωε(ϕε)(Tε1A)ϕ¯1Qεk=1nzkckεϕ¯dx+αεωε[[ϕε]][[ϕ¯]]dSx=ωε(Tε1g)[[ϕ¯]]dSx,t(0,τ),(22b) for all test functions c¯iH1(0,τ;L2(Qε)×L2(ωε))L2(0,τ;H1(Qε)×H1(ωε)) and ϕ¯H1(Qε)×H1(ωε) such that c¯i=0 on (0,τ)×∂Ω and ϕ¯=0 on ∂Ω. In (Equation22a) the following notation was used for short: (23) Υj(c):=cjNAzj1Ck=1nzkck.(23) The time-derivative in (Equation22a) is understood in the weak sense such that 0τciεtc¯idt=0τciεc¯itdt+ciεc¯i|t=0τ. The factor 1/ε in the left-hand side of (Equation22b) comes from the discontinuous Poincaré inequality, see [Citation28, Lemma 3.3], that holds for fH1(Qε)×H1(ωε) with f=0 on ∂Ω: (24) fH1(Qε)×H1(ωε)2=Qεωε(f2+|f|2)dxkDPQεωε|f|2dx+1εωε[[f]]2dSx.(24) Under the assumptions made here, the following theorem is based on [Citation9,Citation10].

Theorem 3.1

Well-posedness

(i) There exists a solution (Equation21) of the generalized Poisson–Nernst–Planck problem (Equation22) satisfying the total mass balance (Equation14a). The positivity (Equation14b) is guaranteed locally at least for small τ(ε)τ0>0 for all ε0, where the uniform bound is provided by the local in time positivity ci0>0 of the limit solution of (Equation64). Moreover, if instead of (Equation11) the stronger assumption miDij=δijD~,i,j=1,,n, is imposed, then the non-negativity ciε0 is guaranteed globally for all τ>0.

(ii) The solution satisfies the following a-priori estimates, which are uniform in ε(0,ε0) for ε0>0 sufficiently small, with constants Kϕ,γc,Kc>0: (25a) |||cε|||2:=cεL(0,τ;L2(Qε)×L2(ωε))2+cεL2(0,τ;H1(Qε)×H1(ωε))2Kc+γcKϕ,(25a) (25b) ϕεL(0,τ;H1(Qε)×H1(ωε))2Kϕ.(25b)

4. Asymptotic analysis

We aim to homogenize the generalized PNP problem (Equation22) and to get residual error estimates. This task needs the asymptotic analysis as ε0+.

In the following, the Poincaré and trace inequalities will be used. For functions uH1(O) defined in a connected domain O=Y,Π,ω, there exists KP(O)>0 such that (26) uuOL2(O)2KP(O)uL2(O)2,uO:=1|O|Ou(y)dy.(26) In the particles ωε, applying to (Equation26) with O=ω the averaging operator Tε1 such that f=Tε1uH1(ωε) and using the integration rules (Equation4e) and (Equation4f) provides (27) 1ε2lIεffωεlL2(ωεl)2KP(ω)fL2(ωε)2.(27) In the pore phase, for fH1(Qε), f=0 on ∂Ω, the Poincaré inequality holds (28) fL2(Qε)2KP(Qε)fL2(Qε)2,KP(Qε)>0.(28) In the following, we write a unique Poincaré constant KP in (Equation26)–(Equation28) for short.

For a discontinuous across the interface ω function uH1(Π)×H1(ω), the trace theorem provides the following estimate with a constant K0>0: (29) [[u]]yL2(ω)2K0uL2(Π)×L2(ω)2+yuL2(Π)×L2(ω)2=K0uH1(Π)×H1(ω)2.(29) For fH1(Qε)×H1(ωε) in the two-phase domain such that f=Tε1u, applying the trace theorem and the integration rules (Equation4e), (Equation4f), and (Equation8d), from (Equation29) it follows (30) 1ε[[f]]L2(ωε)2K01ε2fL2(Ωε)×L2(ωε)2+fL2(Ωε)×L2(ωε)2K0ε2fH1(Ωε)×H1(ωε)2.(30) Based on [Citation13,Citation24], we formulate an auxiliary lemma for homogenization over the pore part Qε of the reference domain Ω.

Lemma 4.1

Asymptotic formula for restriction to pores

For given functions f,qH1(Ω), which are continuous over the interface ωε, the asymptotic representation in the pore space Qε with the porosity ϰ:=|Π|/|Y| holds as ε0+: (31) QεfqdxϰΩfqdx=O(ε).(31)

4.1. Cell problems

For homogenization of the periodic function g and periodic matrices A and D, three auxiliary problems below are formulated in the two-phase unit cell Πω.

First, for the interface data g we set the cell problem for Λ(y) as follows: (32a) divy(yΛ)A=0in Πω,(32a) (32b) [[(yΛ)A]]yν=0,(yΛ)Aν+α[[Λ]]y=gon ω,(32b) (32c) (yΛ)A(,k)|yk=0=(yΛ)A(,k)|yk=1,Λ|yk=0=Λ|yk=1for k=1,,d.(32c) Using the space of periodic functions H#1(Π):={uH1(Π):u|yk=0=u|yk=1, k=1,,d} we get the weak formulation of (Equation32): Find ΛH#1(Π)×H1(ω) such that (33) Πω(yΛ)Ayudy+ωα[[Λ]]y[[u]]ydSy=ωg[[u]]ydSy(33) for all test functions uH#1(Π)×H1(ω). Based on the standard elliptic theory, there exists a solution Λ defined up to a constant value in the cell Y .

Lemma 4.2

Asymptotic formula for periodic interface data

For a given function gL(0,τ;L2(ω)) and fixed ε>0, a periodic function (Tε1Λ)(x)=Λ({x/ε}) defined in (Equation33) satisfies the following asymptotic relation: (34) Qεωεε((Tε1Λ))(Tε1A)ϕ¯dx+ωεα[[Tε1Λ]][[ϕ¯]]dSx=ωε(Tε1g)[[ϕ¯]]dSx+O(ε),(34) for all test functions ϕ¯H1(Qε)×H1(ωε) such that ϕ¯=0 on ∂Ω.

Proof.

For ϕ¯H1(Qε)×H1(ωε) such that ϕ¯=0 on ∂Ω, we multiply (Equation32a) with Tεϕ¯(x,y) and integrate by parts for yΠω using (Equation32b) such that 0=Πωdivy(yΛ)A(Tεϕ¯)dy=Πω(yΛ)Ay(Tεϕ¯)dy+ω(α[[Λ]]yg)[[Tεϕ¯]]ydSyY(yΛ)Aν(Tεϕ¯)dSy. After integration of this relation over xΩε, using the periodicity in (Equation32c) for (yΛ)Aν on Y, we get (35) ΩεΠω(yΛ)Ay(Tεϕ¯)dydx+Ωεω(α[[Λ]]yg)[[Tεϕ¯]]ydSydx=ΩεYΩε(yΛ)Aν(Tεϕ¯)dSydx.(35) Adding to the first integral over Ωε in the left-hand side of (Equation35) the term in ΩΩε, which is of the order O(ε), we apply to (Equation35) the integration rules (Equation4f) and (Equation8c) from Section 2. The resulting integral in the right-hand side of (Equation35) is integrated by parts in ΩΩε using ϕ¯=0 on ∂Ω such that ε|Y|Ωεε(Tε1Λ)(Tε1A)νϕ¯dSx=ε2|Y|ΩΩεdiv[(Tε1Λ)(Tε1A)]ϕ¯(Tε1Λ)(Tε1A)ϕ¯dx=O(ε), where the factor ε2 is cancelled according to (Equation4f), and |ΩΩε|=O(ε). It follows (Equation34) and finishes the proof.

Based on Lemma 4.2, the corrector ε(Tε1Λ) will appear in expansion (Equation66b) of the solution ϕε of the inhomogeneous equation (Equation22b) after homogenization.

Second, for the permittivity matrix A(y) we formulate the following boundary value problem for a vector-function Φ=(Φ1,,Φd)(y) in the two-phase unit cell: (36a) divy(yΦ+I)A=0in Πω,(36a) (36b) [[(yΦ+I)A]]yν=0,(yΦ+I)Aν+α[[Φ]]y=0on ω,(36b) (36c) (yΦ+I)A(,k)|yk=0=(yΦ+I)A(,k)|yk=1,Φ|yk=0=Φ|yk=1for k=1,,d.(36c) In (Equation36), the divergence divy is taken for every Φi(y), the notation yΦ(y)Rd×d for yΠω stands for the matrix of derivatives with entries (yΦ)ij=Φi/yj for i,j=1,,d, and IRd×d is the identity matrix.

The weak form of (Equation36) implies: Find Φ(H#1(Π)×H1(ω))d such that (37) Πω(yΦ+I)Ayudy+ωα[[Φ]]y[[u]]ydSy=0(37) for all uH#1(Π)×H1(ω). A solution Φ exists up to a constant in the cell Y .

Based on Φ, another corrector will appear in the asymptotic expansion (Equation66b) as argued in the following lemma.

Lemma 4.3

Asymptotic formula for periodic permittivity matrix

  1. For the solution Φ of the cell problem (Equation37) the following representation holds: (38) (yΦ(y)+I)A(y)=A0+B1(y),yΠω,(38) where the constant d-by-d matrix A0 is given in the cell Y by the averaging A0:=(yΦ+I)AΠω, it is symmetric positive definite: (39) there exist a_00 such that ξA0ξa_0|ξ|2for ξRd.(39) The d-by-d matrix B1 in (Equation38) has the form in Πω: (40) (B1)kl=m=1dbklm,m(1),where bklm,m(1):=bklm(1)ym,(40) which components are skew-symmetric: (41) bklm(1)+bkml(1)=0,k,l,m=1,,d,(41) the average B1Πω=0, and the matrix B1 is divergence-free as follows (42) l,m=1dbklm,lm(1)=0,where bklm,lm(1):=2bklm(1)ylym,(42) and satisfies the following conditions at the interface: (43) [[B1]]yν=0,(A0+B1)ν=α[[Φ]]yon ω.(43)

  2. Assume that the solution of (Equation36) is such that Φ and yΦ are uniformly bounded in Πω. For given functions ϕ¯H1(Qε)×H1(ωε) and ϕ0H3(Ω), the following asymptotic formula holds with an arbitrary weight δ>0: (44) IA0Qεωε(ϕ1)(Tε1A)ϕ¯dx+ωεαε[[ϕ1]][[ϕ¯]]dSxωεδ[[ϕ¯]]2dSx+Kδε,with some K>0,for IA0:=Qεωε(ϕ0)A0ϕ¯dx+ωε(ϕ0)(A0ν)[[ϕ¯]]dSx,(44) where the notation ϕ1:=ϕ0+ε(ϕ0)(Tε1Φ)ηΩε, and ηΩε is a smooth cut-off function supported in Ωε and equals one outside an ϵ-neighbourhood of Ωε.

Proof.

(i) For the vector-valued solution Φ of (Equation37), the representation (Equation38) with properties (Equation39)–(Equation42) follows from the Helmholtz theorem, see [Citation17, Section 1.1]. The interface conditions (Equation43) are obtained after substitution of (Equation38) into (Equation36b) because of [[A0]]=0.

(ii) Let ϕ¯H1(Qε)×H1(ωε) and ϕ0H3(Ω) be given. To prove (Equation44), we rewrite IA0 in virtue of the integration rules (Equation4f) and (Equation8c) in the micro-variable y: (45) IA0=1ε2|Y|ΩΠω(y(Tεϕ0))(TεA0)y(Tεϕ¯)dy+ω(y(Tεϕ0))(TεA0)ν[[Tεϕ¯]]ydSydx.(45) For the constant matrix A0=TεA0 holds. Then, expressing A0 from (Equation38), using the product rule (y(Tεϕ0))yΦ=(y[(y(Tεϕ0))Φ])Φy(y(Tεϕ0)), the chain rule εTε(ϕ0)=y(Tεϕ0), and the notation ϕ~1:=ϕ0+ε(ϕ0)(Tε1Φ), we rearrange the following terms: (y(Tεϕ0))(TεA0)=(y(Tεϕ0))(A+(yΦ)AB1)=(y(Tεϕ~1))AΦy(y(Tεϕ0))A(y(Tεϕ0))B1. Taking into account this formula, IA0 in (Equation45) is equivalent to: (46) IA0=1ε2|Y|ΩΠω(y(Tεϕ~1))Ay(Tεϕ¯)Φy(y(Tεϕ0))Ay(Tεϕ¯)dy+ω(y(Tεϕ0))A0ν[[Tεϕ¯]]ydSy+IB1dx,(46) with the integral IB1 written component-wisely as follows: IB1:=Πω(y(Tεϕ0))B1y(Tεϕ¯)dy=Πωk,l,m=1d(Tεϕ0),kbklm,m(1)(Tεϕ¯),ldy. Recalling B1 from (Equation40), we integrate by parts IB1 and use the fact that B1 is divergence-free according to (Equation42) such that k,l,m=1d(Tεϕ0),kbklm,lm(1)=0 to get (47) IB1=Πωk,l,m=1d(Tεϕ0),klbklm,m(1)(Tεϕ¯)dy+ω(y(Tεϕ0))B1ν[[Tεϕ¯]]ydSyY(y(Tεϕ0))B1ν(Tεϕ¯)dSy.(47) After integration by parts the second time and rearranging the mixed derivatives (Tεϕ0),klm such that k,l,m=1d(Tεϕ0),klmbklm(1)=0 because bklm(1) is skew-symmetric as written in (Equation41), we proceed (Equation47): IB1=Πωk,l,m=1d(Tεϕ0),klbklm(1)(Tεϕ¯),mdy+ω(y(Tεϕ0))B1νk,l,m=1d(Tεϕ0),klbklm(1)νm[[Tεϕ¯]]dSy+IY, where IY:=Y(k,l,m=1d(Tεϕ0),klbklm(1)νm(y(Tεϕ0))B1ν)(Tεϕ¯)dSy.

Substituting the expression of IB1 into (Equation46) and using the formula at ω: α[[Tεϕ~1]]y=α[[Tεϕ0]]y+(y(Tεϕ0))α[[Φ]]y=(y(Tεϕ0))(A0+B1)ν following from (Equation43) and [[Tεϕ0]]y=0, with the help of the integration rules (Equation4f) and (Equation8c) we rewrite IA0 again with respect to the macro-variable x in the form: (48) IA0=Qεωε(ϕ~1)ε(Tε1Φ)x(ϕ0)(Tε1A)ϕ¯k,l,m=1dk,l,m=1dεϕ,kl0(Tε1bklm(1))ϕ¯,mdx+ωεαε[[ϕ~1]]k,l,m=1dεϕ,kl0(Tε1bklm(1))νm[[ϕ¯]]dSx+IΩε,(48) where the last two terms in the integral over Qεωε have the asymptotic order O(ε), and IY is transformed to the integral over Ωε such that IΩε:=Ωεk,l,m=1dεϕ,kl0(Tε1bklm(1))νm(ϕ0)ε(Tε1B1)νϕ¯dSx. Here, the factor ϵ appears due to the integration rule over the boundary Y analogously to (Equation8c), the chain rule gives (Tεϕ0),kl=ε2Tε(ϕ,kl0) and y(Tεϕ0)=εTε(ϕ0), while in the second term ϵ appears since (49) B1=m=1dymbklm(1)=m=1dεxm(Tε1bklm(1))=εTε1B1.(49) By this, the factor ε2 is cancelled by division by ε2 in (Equation46).

We estimate the interface term in the integral over ωε in the right-hand side of the Equation (Equation48) by Young's inequality with a weight δ>0 as follows: (50) ωεk,l,m=1dεϕ,kl0(Tε1bklm(1))νm[[ϕ¯]]dSxωεε24δk,l,m=1dϕ,kl0(Tε1bklm(1))νm2k,l,m=1dϕ,kl0(Tε1bklm(1))νm2+δ[[ϕ¯]]2dSxωεδ[[ϕ¯]]2dSx+Kδε,K>0,(50) since |ωε|=O(ε1). Applying Green's formula in the boundary layer ΩΩε and using ϕ¯=0 on ∂Ω leads to the asymptotic expansion of the boundary term: (51) IΩε=ΩΩεk,l,m=1dεϕ,kl0(Tε1bklm(1))ϕ¯,m+ϕ,kl0(εTε1bklm(1)),mϕ¯(ϕ0)(εTε1B1)ϕ¯div(ϕ0)(εTε1B1)ϕ¯dx=O(ε).(51) Here the ϵ-order is due to the fact that |ΩΩε|=O(ε), the uniform boundedness of εTε1B1 and the chain rule Tε1bklm,m(1)=(εTε1bklm(1)),m according to (Equation49).

Gathering in (Equation48) the asymptotic terms of the same order ϵ and accounting for formulas (Equation50) and (Equation51), the following estimate takes place with some K>0: (52) IA0Qεωε(ϕ~1)(Tε1A)ϕ¯dxωεαε[[ϕ~1]][[ϕ¯]]dSxωεδ[[ϕ¯]]2dSx+Kδε.(52) For a cut-off function ηΩε supported in Ωε we set ϕ1:=ϕ0+ε(ϕ0)(Tε1Φ)ηΩε such that ϕ1=0 in ΩΩε, the jump [[ϕ1]]=[[ϕ~1]] at ωε, and (53) ϕ1ϕ~1H1(Qε)×H1(ωε)=O(ε).(53) From (Equation52) and (Equation53) if follows (Equation44) and the assertion of Lemma 4.3.

Third, for a diffusivity matrix D corresponding to the assumption (Equation12) in Theorem 5.1 below, in analogy with (Equation36), we establish the cell problem for N=(N1,,Nd)(y): (54a) divy(yN+I)D=0in Πω,(54a) (54b) [[(yN+I)D]]yν=0,(yN+I)Dν=0on ω,(54b) (54c) (yN+I)D(,k)|yk=0=(yN+I)D(,k)|yk=1,N|yk=0=N|yk=1for k=1,,d.(54c) The system (Equation54) differs from (Equation36) by the interface condition and implies the following weak formulation: Find a vector-function N(H#1(Π)×H1(ω))d such that (55) Πω(yN+I)Dyudy=0(55) for all test functions uH#1(Π)×H1(ω). A solution of (Equation55) exists and is defined up to a piecewise constant in Πω. Moreover, since ω¯Y is assumed, this fact follows that N=−y and yN=I in ω. Based on N, the following lemma justifies the use of the corrector ε(ci0)(Tε1N) in the formula (Equation66a).

Lemma 4.4

Asymptotic formula for periodic diffusivity matrix

  1. For the solution N of the cell problem (Equation55) the following representation holds: (56) (yN(y)+I)D(y)=D0+B2(y),(56) where the d-by-d matrix D0 is constant in the cell Y and given by D0:=(yN+I)DΠω=(yN+I)DΠ, it is symmetric positive definite: (57) there exist d_00 such thatξD0ξd_0|ξ|2for ξRd.(57) The d-by-d matrix B2 has the following form in Πω: (58) (B2)kl=m=1dbklm,m(2),k,l=1,,d.(58) Its components bklm(2) are skew-symmetric, B2Πω=0, and B2 is divergence-free in the manner of (Equation41) and (Equation42). At the interface the conditions hold (59) [[B2]]yν=0,(D0+B2)ν=0on ω.(59)

  2. Assume N(W1,(Π)×W1,(ω))d. For ci¯L2(0,τ;H1(Qε)×H1(ωε)) such that c¯i=0 on ∂Ω and arbitrary ci0L2(0,τ;H3(Ω)), the following asymptotic formula with ci1:=ci0+ε(ci0)(Tε1N)ηΩε holds (60) 0τID0dt=0τQεωε(ci1)(Tε1D)c¯idxdt+O(ε),ID0:=Qεωε(ci0)D0c¯idx+ωε(ci0)D0ν[[c¯i]]dSx.(60)

Proof.

The proof is analogous to those from the previous Lemma 4.3 until (Equation47). Indeed, we derive similar to (Equation45) and (Equation46) formulas in micro-variables: (61) ID0=1ε2|Y|ΩΠω(y(Tεc~i1))Dy(Tεc¯i)Ny(y(Tεci0))Dy(Tεc¯i)dy+ω(y(Tεci0))D0ν[[Tεc¯i]]ydSy+IB2dx,(61) with c~i1:=ci0+ε(ci0)(Tε1N) and IB2:=Πω(y(Tεci0))B2y(Tεc¯i)dy. Likewise (Equation47), integration by parts of IB2 follows that (62) IB2=Πωk,l,m=1d(Tεci0),klbklm,m(2)(Tεc¯i)dy+ω(y(Tεci0))B2ν[[Tεc¯i]]ydSyY(y(Tεci0))B2ν(Tεc¯i)dSy.(62) After substitution of (Equation62) in (Equation61), the integral over ω disappears due to the interface condition (Equation59).

Returning to the micro-variables x with the help of the chain rule (/ym)Tε=εTε(/xm), the second term in the integral over Πω in (Equation61) has the asymptotic order O(ε). The integral over Y in (Equation62) divided by ε2 is transformed to the integral over Ωε with the factor 1/ε, and after integration by parts in the boundary layer ΩΩε, it is of the order O(ε), too.

The principal difference from Lemma 4.3 consists in estimation of the domain integral in IB2.

By adding and subtracting the averaged values, we rewrite equivalently Πωk,l,m=1d(Tεci0),klbklm,m(2)(Tεc¯i)dy=I1+I2, using the property B2Πω=0, and I1:=Πωk,l,m=1d(Tεci0),klbklm,m(2)(Tεc¯iTεc¯iΠω)dy,I2:=Tεc¯iΠωΠωk,l,m=1d[(Tεci0),kl(Tεci0),klΠω]bklm,m(2)dy. We rewrite I1 and I2 in the macro-variable x in all local cells using the integration rules (Equation4c) and (Equation8c), applying the chain rule (/ym)Tε=εTε(/xm) to ci0 and to B2 (see (Equation49)), then apply to the result the Cauchy–Schwarz inequality and the Poincaré inequality (Equation27). First, there are some constants 0K1K2 and K30 such that 1ε2|Y|ΩI1dx=jIεΠεjωεjk,l,m=1dci,kl0(εTε1bklm,m(2))(c¯ic¯iΠεjωεj)dxK1ci0H2(Πεωε)B2L(Πω)c¯iL2(Πεωε)K2ci0H2(Πεωε)(K3+yNL(Πω))εc¯iL2(Πεωε)=O(ε), where we have used the fact that the integral over the boundary layer ΩΩε of Tε1(Tεc¯iTεc¯iΠω) is zero due to the definition of the operator Tε1 in ΩΩε. Similarly, there exists K40 such that 1ε2|Y|ΩI2dxK4k,l=1dε(ci,kl0)L2(Πεωε)(K3+yNL(Πω))c¯iL2(Πεωε)=O(ε). Finally, we integrate the estimate of ID0 over the time t(0,τ) for further use.

The functions ci0 and ϕ0 will associate the averaged solution in the homogenization problem presented in the next section.

5. The main homogeneous result

In this section, we establish the averaged PNP equations for the functions (c0,ϕ0)(t,x) in the time-space domain (0,τ)×Ω as follows: (63a) ci0tdivkBΘ(ci0)D0=0for i=1,,n,(63a) (63b) div(ϕ0)A0=ϰk=1nzkck0,where ϰ=|Π||Y|,(63b) which are supported by the Dirichlet boundary and initial conditions: (63c) ci0=ciDandϕ0=ϕDon (0,τ)×∂Ω,ci0=ciininΩ.(63c) In (Equation63), the averaged matrices A0=(yΦ+I)AΠω and D0=(yN+I)DΠ are from Lemma 4.3 and Lemma 4.4, the matrix D is from (Equation12), the vectors N and Φ are the solutions of the two-phase cell problems (Equation55) and (Equation37), respectively.

From the standard existence theorems on elliptic and parabolic systems, the solution ϕ0L(0,τ;H1(Ω)) and ci0L(0,τ;L2(Ω))L2(0,τ;H1(Ω)) of the linear problem (Equation63) exists and fulfils the following variational equations: (64a) 0τΩci0tc¯i+kBΘ(ci0)D0c¯idxdt=0,for i=1,,n,(64a) (64b) Ω(ϕ0)A0ϕ¯ϰk=1nzkck0ϕ¯dx=0,(64b) for all test functions c¯iL2(0,τ;H01(Ω)) and ϕ¯H01(Ω).

The main result of this paper is the following theorem.

Theorem 5.1

Averaged problem and correctors

Let the solutions N, Φ of the two-phase cell problems (Equation55), (Equation37), and yN, yΦ be uniformly bounded in Πω, the averaged solutions ϕ0L(0,τ;H3(Ω)) and ci0L2(0,τ;H3(Ω)), i=1,,n. Then a solution (cε,ϕε) of the inhomogeneous PNP problem (Equation22) and the solution (c0,ϕ0) of the homogeneous PNP problem (Equation64) satisfy the residual error estimates: (65) |||cεc1|||2=O(ε),ϕεϕ2L(0,τ;H1(Qε)×H1(ωε))2=O(ε),(65) with the norm |||||| defined in (Equation25a), and the approximate functions are (66a) ci1:=ci0+ε(ci0)(Tε1N)ηΩε,(66a) (66b) ϕ2:=ϕ1+ε(Tε1Λ)ηΩε,ϕ1:=ϕ0+ε(ϕ0)(Tε1Φ)ηΩε.(66b) In (Equation66), the vector Λ is a solution of the two-phase cell problem (Equation33), and ηΩε is the cut-off function from Lemmas 4.3 and 4.4.

Proof.

Based on the asymptotic results of Section 3, we will prove the error estimates (Equation65). In particular, this will justify the averaged problem (Equation63).

Estimate of cεc1. We start with derivation of an asymptotic equation for ci1 as i=1,,n. We apply to div((ci0)D0) Green's formulas on the pore phase: (67a) Qεdiv(ci0)D0c¯i+(ci0)D0c¯idx=ωε+(ci0)D0νc¯idSx,(67a) for all c¯iH1(Qε) such that c¯i=0 on ∂Ω, and on the solid phase: (67b) ωεdiv(ci0)D0c¯i+(ci0)D0c¯idx=ωε(ci0)D0νc¯idSx,(67b) for all c¯iH1(ωε). Summing up the Equations (Equation67), using the diffusion equation (Equation63a) and the continuity of (ci0)D0ν across ωε, the variational problem (Equation64a) in Ω can be expressed equivalently over the two-phase domain as follows: (68) 0τQεωεci0tc¯i+kBΘ(ci0)D0c¯idxdt+0τωεkBΘ(ci0)D0ν[[c¯i]]dSxdt=0,(68) for all discontinuous over ωε test functions c¯iL2(0,τ;H1(Qε)×H1(ωε)) such that c¯i=0 on ∂Ω. Further, we employ the asymptotic arguments as ε0+.

We apply to the left-hand side of (Equation68) the asymptotic formula (Equation60) from Lemma 4.4, which implies: (69) 0=0τQεωεci0tc¯i+kBΘ(ci1)(Tε1D)c¯idxdt+O(ε),(69) where ci1 is defined in (Equation66a). In virtue of the relation ci1t=t[ci0+ε(ci0)(Tε1N)ηΩε]=ci0t+O(ε), then (Equation69) can be rewritten in terms of ci1 in the asymptotically equivalent form: (70) 0τQεωεci1tc¯i+kBΘ(ci1)(Tε1D)c¯idxdt=O(ε).(70) We continue with an asymptotic expansion of the perturbed problem (Equation22a). Due to the assumption (Equation12) on the diffusivity matrices and the uniform estimate |Υj(cε)|(|zj|+i=1n|zi|)C/NA, which follows that εκΥj(cε)=O(ε) for κ1, the Equation (Equation22a) is expressed in the asymptotic form: (71) 0τQεωεciεtc¯i+kBΘ(ciε)(Tε1D)c¯idxdt=0τωεε2gi(cˆε,ϕˆε)[[c¯i]]dSxdt+O(ε).(71) Since |ωε|=O(ε1), the interface integral over ωε in (Equation71) is estimated by Young's inequality due to the boundedness property (Equation20c) and the trace theorem (Equation30): (72) ωεε2gi(cˆε,ϕˆε)[[c¯i]]dSxε214ωε|gi(cˆε,ϕˆε)|2dSx+ωε|[[c¯i]]|2dSxε2Kg4|ωε|+K0εc¯iH1(Qε)×H1(ωε)2=O(ε).(72) Next, we subtract the Equation (Equation70) from (Equation71) and utilize (Equation72) to obtain that (73) 0τQεωε(ciεci1)tc¯i+kBΘ(ciεci1)(Tε1D)c¯idxdt=O(ε).(73) Integrating by parts over time in the first term in (Equation73) implies (74) Qεωε(ciεci1)22t=0τ+0τkBΘ(ciεci1)(Tε1D)(ciεci1)dtdx=O(ε).(74) The initial difference here (ciεci1)|t=0=ε(ciin)(Tε1N)ηΩε=O(ε). Using the uniform positive definiteness (Equation13) of D, after taking the supremum over τ(0,τ¯) and summing up (Equation74) over i=1,,n we arrive at the first estimate in (Equation65): (75) i=1nsupt(0,τ¯)Qεωε(ciεci1)2dx+0τ¯Qεωε(ciεci1)2dxdt=O(ε).(75) In particular, applying the triangle inequality for ci1 given by the sum in (Equation66a), due to the uniform boundedness of N, yN, and c0L2(0,τ;H1(Ω))d, from (Equation75) it follows the estimate which will be used further in (Equation82): (76) cεc0L(0,τ;L2(Qε)×L2(ωε))22ncεc1L(0,τ;L2(Qε)×L2(ωε))2+2nε2(c0)(Tε1N)ηΩεL(0,τ;L2(Qε)×L2(ωε))2=O(ε).(76)

Estimate of ϕεϕ2. Similarly to (Equation67), we apply to the term div((ϕ0)A0) the following Green's formulas on the both phases Qε and ωε: (77a) Qε(ϕ0)A0ϕ¯+div(ϕ0)A0ϕ¯dx=ωε+(ϕ0)A0νϕ¯dSx,(77a) (77b) ωε(ϕ0)A0ϕ¯+div(ϕ0)A0ϕ¯dx=ωε(ϕ0)A0νϕ¯dSx,(77b) for test functions ϕ¯H1(Qε) such that ϕ¯=0 at ∂Ω, and ϕ¯H1(ωε), respectively. We sum up the Equations (Equation77), use the Poisson equation (Equation63b) and the continuity of (ϕ0)A0ν across the interface ωε. Applying the asymptotic formula (Equation31) from Lemma 4.1 we rewrite (Equation64b) over the two-phase domain as follows: (78) Qεωε(ϕ0)A0ϕ¯1Qεk=1nzkck0ϕ¯dx+ωε(ϕ0)A0ν[[ϕ¯]]dSx=O(ε),(78) for all test functions ϕ¯H1(Qε)×H1(ωε) such that ϕ¯=0 at ∂Ω.

Applying the inequality (Equation44) from Lemma 4.3 with ϕ1:=ϕ0+ε(ϕ0)(Tε1Φ)ηΩε proceeds the expansion (Equation78) with some K>0 as (79) Qεωε(ϕ1)(Tε1A)ϕ¯1Qεk=1nzkck0ϕ¯dx+ωεαε[[ϕ1]][[ϕ¯]]dSxωεδ[[ϕ¯]]2dSx+Kε.(79) Next, we add to (Equation79) the Equation (Equation34) describing Λ from Lemma 4.2 and use the definition of ϕ2=ϕ1+ε(Tε1Λ)ηΩε to get (80) Qεωε(ϕ2)(Tε1A)ϕ¯1Qεk=1nzkck0ϕ¯dx+ωεαε[[ϕ2]]Tε1g[[ϕ¯]]dSxωεδ[[ϕ¯]]2dSx+Kε.(80) The subtraction of (Equation80) from the perturbed equation (Equation22b) implies that (81) Qεωε(ϕεϕ2)(Tε1A)ϕ¯dx+ωεαε[[ϕεϕ2]][[ϕ¯]]dSxQεk=1nzk(ckεck0)ϕ¯dxωεδ[[ϕ¯]]2dSx+Kε.(81) After substitution in (Equation81) the test function ϕ¯:=ϕεϕ2, which is zero at ∂Ω, using Young's inequality with a weight δ1>0 and applying the asymptotic bound (Equation76) of (ciεci0), we obtain the asymptotic inequality for δ<α/ε0 such that α/εδ>(αδε0)/ε>0 for 0<ε<ε0: (82) 0Qεωε(ϕεϕ2)(Tε1A)(ϕεϕ2)dx+ωεαεδ[[ϕεϕ2]]2dSxZ¯22δ1k=1nckεck0L2(Qε)2+δ12ϕεϕ2L2(Qε)2+Kε=δ12ϕεϕ2L2(Qε)2+O(ε),(82) where Z¯:=maxk{1,,n}|zk|. For δ1 chosen small enough, using the uniform positive definiteness of A in (Equation10) and the lower bound (Equation24), taking the supremum over t(0,τ) in (Equation82) follows the second estimate in (Equation65) and finishes the proof.

6. Discussion

Passing to the limit in (Equation14), we derive the total mass balance and the non-negativity for the averaged species concentrations c0.

According to the governing relations (Equation9c) and (Equation9d), we can introduce the entropy variables (μ10,,μn0), v0=(v10,,vd0), and p0 corresponding to the solution of the averaged problem (Equation63) as follows: μi0:=kBΘln(βici0);ηv0+p0=j=1nzjcj0ϕ0,div v0=0. We observe the following technical assumptions used for the homogenization:

  • the asymptotic factor εκ, κ1, in the electrochemical potentials μi in (Equation9c);

  • the asymptotic factor ε2 by the interface reactions gi(,) in (Equation19a);

  • asymptotic decoupling of the diffusivity matrices Dij in (Equation12).

Our future work is pointed towards possible relaxing these assumptions.

Acknowledgements

The authors thank two referees for the comments which helped to improve the manuscript.

Disclosure statement

No potential conflict of interest was reported by the authors.

Additional information

Funding

The authors are supported by the Austrian Science Fund (FWF) Project P26147-N26: ‘Object identification problems: numerical analysis’ (PION). V.A.K. thanks the Austrian Academy of Sciences (OeAW), [the RFBR and JSPS research project 19-51-50004], and A.V.Z. thanks IGDK1754 for partial support.

References

  • Dreyer W, Guhlke C, Müller R. Overcoming the shortcomings of the Nernst–Planck model. Phys Chem Chem Phys. 2013;15:7075–7086. doi: 10.1039/c3cp44390f
  • Dreyer W, Guhlke C, Müller R. Modeling of electrochemical double layers in thermodynamic non-equilibrium. Phys Chem Chem Phys. 2015;17:27176–27194. doi: 10.1039/C5CP03836G
  • Fuhrmann J. Comparison and numerical treatment of generalized Nernst–Planck models. Comput Phys Commun. 2015;196:166–178. doi: 10.1016/j.cpc.2015.06.004
  • Fuhrmann J, Guhlke C, Linke A, et al. Models and numerical methods for electrolyte flows. WIAS Preprint. 2018;2525. Available from: http://www.wias-berlin.de/preprint/2525/wias_preprints_2525.pdf.
  • Burger M, Schlake B, Wolfram MT. Nonlinear Poisson–Nernst–Planck equations for ion flux through confined geometries. Nonlinearity. 2012;25:961–990. doi: 10.1088/0951-7715/25/4/961
  • Herz M, Ray N, Knabner P. Existence and uniqueness of a global weak solution of a Darcy–Nernst–Planck–Poisson system. GAMM–Mitt. 2012;35:191–208. doi: 10.1002/gamm.201210013
  • Roubíček T. Incompressible ionized non-Newtonean fluid mixtures. SIAM J Math Anal. 2007;39:863–890. doi: 10.1137/060667335
  • Roubíček T. Incompressible ionized fluid mixtures: a non-Newtonian approach. IASME Trans. 2005;2:1190–1197.
  • Kovtunenko VA, Zubkova AV. Solvability and Lyapunov stability of a two-component system of generalized Poisson–Nernst–Planck equations. In: Maz'ya V, Natroshvili D, Shargorodsky E, Wendland WL, editors. Recent trends in operator theory and partial differential equations (The Roland Duduchava Anniversary Volume), Operator theory: advances and applications; Vol. 258. Basel: Birkhaeuser; 2017. p. 173–191.
  • Kovtunenko VA, Zubkova AV. On generalized Poisson–Nernst–Planck equations with inhomogeneous boundary conditions: a-priori estimates and stability. Math Meth Appl Sci. 2017;40:2284–2299.
  • Kovtunenko VA, Zubkova AV. Mathematical modeling of a discontinuous solution of the generalized Poisson–Nernst–Planck problem in a two-phase medium. Kinet Relat Mod. 2018;11(1):119–135. doi: 10.3934/krm.2018007
  • Allaire G, Brizzi R, Dufrêche JF, et al. Ion transport in porous media: derivation of the macroscopic equations using upscaling and properties of the effective coefficients. Comp Geosci. 2013;17:479–495. doi: 10.1007/s10596-013-9342-6
  • Belyaev AG, Pyatnitskii AL, Chechkin GA. Averaging in a perforated domain with an oscillating third boundary condition. Mat Sb. 2001;192:3–20. doi: 10.4213/sm576
  • Evendiev M, Zelik SV. Attractors of the reaction-diffusion systems with rapidly oscillating coefficients and their homogenization. Ann Instit H Poincare. 2002;19:961–989. doi: 10.1016/S0294-1449(02)00115-4
  • Mielke A, Reichelt S, Thomas M. Two-scale homogenization of nonlinear reaction-diffusion systems with slow diffusion. J Netw Heterog Media. 2014;9(2):353–382. doi: 10.3934/nhm.2014.9.353
  • Sazhenkov SA, Sazhenkova EV, Zubkova AV. Small perturbations of two-phase fluid in pores: Effective macroscopic monophasic viscoelastic behavior. Sib Èlektron Mat Izv. 2014;11:26–51.
  • Zhikov VV, Kozlov SM, Oleinik OA. Homogenization of differential operators and integral functionals. Berlin: Springer-Verlag; 1994.
  • Bunoiu R, Timofte C. Homogenization of a thermal problem with flux jump. Netw Heterog Media. 2016;11:545–562. doi: 10.3934/nhm.2016009
  • Gagneux G, Millet O. Homogenization of the Nernst–Planck–Poisson system by two-scale convergence. J Elast. 2014;114:69–84. doi: 10.1007/s10659-013-9427-4
  • Gahn M, Neuss-Radu M, Knabner P. Homogenization of reaction-diffusion processes in a two-component porous medium with nonlinear flux conditions at the interface. SIAM J App Math. 2016;76:1819–1843. doi: 10.1137/15M1018484
  • Khoa VA, Muntean A. Corrector homogenization estimates for a non-stationary Stokes–Nernst–Planck–Poisson system in perforated domains. arXiv:1710.09166v1 [math.NA]. 2017; Available from: https://arxiv.org/abs/1710.09166v1.
  • Ray N, Eck C, Muntean A, et al. Variable choices of scaling in the homogenization of a Nernst–Planck–Poisson problem. Vol. 344. Erlangen-Nürnberg: Inst. für Angewandte Mathematik; 2011.
  • Schmuck M, Bazant MZ. Homogenization of the Poisson–Nernst–Planck equations for ion transport in charged porous media. SIAM J Appl Math. 2015;75:1369–1401. doi: 10.1137/140968082
  • Fellner K, Kovtunenko VA. A discontinuous Poisson–Boltzmann equation with interfacial transfer: homogenisation and residual error estimate. Appl Anal. 2016;95:2661–2682. doi: 10.1080/00036811.2015.1105962
  • Efendiev Y, Iliev O, Taralova V. Upscaling of an isothermal li-ion battery model via the homogenization theory. Berichte des Fraunhofer ITWM. 2013;230. Available from: http://publica.fraunhofer.de/documents/N-256468.html.
  • Allaire G, Brizzi R, Dufrêche JF, et al. Role of nonideality for the ion transport in porous media: derivation of the macroscopic equations using upscaling. Phys D. 2014;282:39–60. doi: 10.1016/j.physd.2014.05.007
  • Allaire G, Mikelić A, Piatnitski A. Homogenization of the linearized ionic transport equations in rigid periodic porous media. J Math Phys. 2010;51:123103. doi: 10.1063/1.3521555
  • Hummel HK. Homogenization for heat transfer in polycrystals with interfacial resistances. Appl Anal. 2000;75:403–424. doi: 10.1080/00036810008840857
  • Cioranescu D, Damlamian A, Donato P, et al. The periodic unfolding method in domains with holes. SIAM J Math Anal. 2012;44(2):718–760. doi: 10.1137/100817942
  • Franců J. Modification of unfolding approach to two-scale convergence. Math Bohem. 2010;135:403–412.
  • Khludnev MA, Kovtunenko VA. Analysis of cracks in solids. Southampton–Boston: WIT Press; 2000.