Energetic Stable Discretization for Non-Isothermal Electrokinetics Model

We propose an edge averaged finite element(EAFE) discretization to solve the Heat-PNP (Poisson-Nernst-Planck) equations approximately. Our method enforces positivity of the computed charged density functions and temperature function. Also the thermodynamic consistent discrete energy estimate which resembles the thermodynamic second law of the Heat-PNP system is prescribed. Numerical examples are provided.

Authors

• 1 publication
• 27 publications
• 1 publication
• Helicity-conservative finite element discretization for MHD systems

We construct finite element methods for the magnetohydrodynamics (MHD) s...
07/15/2020 ∙ by Kaibo Hu, et al. ∙ 0

• Sea-ice dynamics on triangular grids

We present a stable discretization of sea-ice dynamics on triangular gri...
05/31/2020 ∙ by Carolin Mehlmann, et al. ∙ 0

• Energy-conserving time propagation for a geometric particle-in-cell Vlasov–Maxwell solver

This paper discusses energy-conserving time-discretizations for finite e...
10/09/2019 ∙ by Katharina Kormann, et al. ∙ 0

• A level-set based space-time finite element approach to the modelling of solidification and melting processes

We present a strategy for the numerical solution of convection-coupled p...
05/07/2021 ∙ by Leonardo Boledi, et al. ∙ 0

• New stable method to solve heat conduction problems in extremely large systems

We present a new explicit and stable numerical algorithm to solve the ho...
08/24/2019 ∙ by Endre Kovács, et al. ∙ 0

• Time discretization of an abstract problem applying to the linearized equations of coupled sound and heat flow

In this paper we deal with an abstract problem which includes the linear...
01/02/2020 ∙ by Shunsuke Kurima, et al. ∙ 0

• Heterogeneity Preserving Upscaling for Heat Transport in Fractured Geothermal Reservoirs

In simulation of fluid injection in fractured geothermal reservoirs, the...
02/20/2017 ∙ by Anna Nissen, et al. ∙ 0

This week in AI

Get the week's most popular data science and artificial intelligence research sent straight to your inbox every Saturday.

1 Background

The study of thermodynamic properties of complex fluid involving electrical solvent is of interest of many biological and physiological applications [21, 16, 7]. For example, electrokinetic microfluidic devices have been widely used in biomedical and biotechnological applications and chemical synthesis, which can be used for pumping and controlling the liquid flow in microfluidic systems or separating constituents suspended in a liquid [16, 12]. While low thermal conductivity devices are widely used, temperature effect should be includeded as it affects the electrical conductivity and can produce substantial heat driven flow. In addition, the interaction between heat flux and electric field could result in Joule heating effect as a byproduct of presense of electric field  [25, 13]. This effect is generated by the ohmic resistance of the electrolyte subjected to electric current and, clearly, the traditional isothermal models aer not adequate in such situation. More recent works account for thermal effects and show that these interactions can be modeled in a thermodynamic consistent way (see [19]).

Another important pool of applications stems from biology, namely, modeling ion channels. Such channels are formed by charged walls and play a fundamental role in controlling and regulating the nervous system. In addition to the three main types of ion channels: voltage gated, ligand gated and chemically gated ion channels, temperature controlled ion channels are also been found and studied in the literatude (see [4, 24]). For example, six members of the mammalian Transient receptor potential ion channels (TRP) respond to varied temperature thresholds. TRPV1, TRPV2, TRPV3, and TRPV4 are heat activated, whereas TRPM8 and TRPA1 are activated by cold [24].

Despite there is vast amount of works devoted on related experimental studies[13, 16], mathematical modeling[27, 22] and simulation[32], [26], thermodynamically consistent models are much less studied and developed.

To begin the discussion of the non-isothermal model we first consider the isothermal case and the Poisson-Nernst-Planck (PNP) system. This model is widely used to describe the transport of charged particles  [23, 1, 11, 28]with low concentration and couples the well-known drift-diffusion equations for the concentration with the Poisson equation for the electrostatic potential with appropriate initial and boundary conditions [20]:

 −∇⋅(ϵ∇ϕ) =ecN∑i=1qiρi (Poisson Equation) ∂ρi∂t =−∇⋅→Ji (Drift Diffusion equation) →Ji =−Di∇ρi−qiμiρi∇ϕ, (Flux)

Here 0 is the diffusivity, is the valence number, and is the mobility of the -th ion species. The PNP system models the interaction of N ionic species through an electrostatic field, usually . is the electrostatic potential, stands for the charge density of the ith species.

Both ion channels and electroosmosis can be modeled by taking the background fluid into consideration and coupling the PNP system with the Navier-Stokes equation. The energetic variational formulation which deals with the isotropic case for the resulting system from a mechanical prospective is well elaborated in[31, 14]. The energy law for this system, as stated in [31] is:

 ddt{∫Ω12ρf|→u|2+N∑i=1ρi(ρi−1)+ϵ2|∇ϕ|2dx+∫ΓRκ2|ϕ|2ds} = −∫Ωμ2|∇u+∇uT2|2+N∑i=1Diρi|∇(logρi+qiϕ)|2dx,

Here, denotes the part of boundary where the Robin boundary condition: is imposed. The total energy consists of the sum of the kinetic energy of the background fluid, the electro potential energy and the free energy of each of the species. The dissipation of the total free energy is driven by the viscosity and diffusion. This coupled system can be derived through energetic variational approach (EVA, introduced in [6]) and is as follows:

 −∇⋅(ϵ∇ϕ) =ecN∑i=1qiρi (Poisson equation) ∂ρi∂t =−∇⋅→Ji (Drift diffusion equation) →Ji =−Di∇ρi−qiμiρi∇ϕ−ρi→u (Flux) ρf(→ut+→u⋅∇→u)+∇p =∇⋅∇u+∇uT2−N∑i=1qiρi∇ϕ (NS equation) ∇⋅→u =0 (Incompressibility)

As is evident from the equations above, this model does not take temperature effect into consideration.

A number of numerical solution approaches for the PNP system and relevant mathematical models models have been considered in the literature. For example, A. Flavell et al. (2014) [10] applied a conservative finite difference scheme which achieves second-order accuracy in both space and time and also conserves total concentration for each ion species. C. Liu et al. (2015) [18] used an energetically stable finite element method for the PNP-NS(Poisson-Nernst-Planck-Navier-Stokes) system. D. Xie et al. (2016) [29] used a non-local finite element method for the PB (Poisson-Boltzmann) equation to tackle the problem of solution singularity caused by point charge term.

In order to capture the temperature change [9], [3] considered a thermodynamic approach for incompressible Newtonian fluid coupled with temperature model. The resulting model is based on the first and second law of thermodynamics and in accordance with the Fourier’s Law assumes that the internal energy is proportional to the temperature and the heat flux is proportional to the temperature gradient. The equation for the temperature is derived through energy conservation and the whole system satisfies the inequality constraint for entropy production and this is in agreement with the second law of thermodynamics.

Some works in biology [27],[22], take both temperature and electrokinetics into consideration by a simple coupling which adds electrokinetic force into the Navier-Stokes equation and use PNP or Poisson-Boltzmann to model the electric field. The temperature effects are incorporated by including the Joule heat force term into the heat equation. One issue with this model is that it is unclear what the energy associated with this system, and, moreover, it is not evidennt whether the second thermodynamics law is violated by the model.

In our study, we are using a more consistent model to capture the interactions mentioned above. According to first law of thermodynamics, heat and work can convert to each other, so it is natural to consider this energy interchange in order to derive the equation for the temperature. By the second law of thermodynamics, the conversion between heat transfer and work follows the rule that keeps the entropy increasing. We follow the ideas in [19] and add this as an inequality constraint for the system. As it turns out, maintaining such inequality constraints is also crucial for the stability of the discretized model. This motivated the choice of discretization and we have applied the Edge Average Finite Element (EAFE, see [30]) discretization, which satisfies the discrete maximum principle for the temperature. The discrete maximum principle turns out to be the key to numerical stability and one of the key results that we present is the energy estimate satisfied by the numerical model.

The paper is organized as follows. We introduce the Heat-PNP equations in section 2 and the corresponding energy law in section 3. In section 4 we propose our discretization and prove an energy estimate for the discretized Heat-PNP system. In section 4 some numerical experiments are provided to validate the stability of our numerical scheme. We also provide numerical tests which show consistency with qualitative phenomena observed experimentally [21].

2 The Heat-PNP equations and their discretization

As it is well known [15] the PNP system is as follows:

 (2.1)

We treat the evolution equations for ions’ concentration as a system of continuity equations. After introducing the flux variable for each ion species, this sytem can be written as:

 ⎧⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪⎨⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪⎩∂ρi∂t=−∇⋅→Ji,→Ji=−Di∇ρi−ziekBTρi∇ϕ,∇⋅(ϵ∇ϕ)=−(ρ0+N∑i=1zieρi) (2.2)

We use the same fashion to rewrite the PNP system with temperature throughout this paper.

We consider the the temperature PNP system in without background fluid. According to [19], this system has the form:

 ⎧⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪⎨⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪⎩∂∂tρi+∇⋅(ρi→ui)=0,i=1,...,N,νiρi→ui=−kB∇(ρiT)−zieρi∇ϕ,i=1,...,N,−∇⋅ϵ∇ϕ=∑iρizie+ρf.(N∑i=1kBCiρi)∂T∂t+(N∑i=1kBCiρi→ui)⋅∇T+(N∑i=1kBρi∇⋅→ui)T=∇⋅k∇T+∑Ni=1νiρi|→ui|2+q,

Here, is the electrostatic potential , stands for the charge density of the ith species, T is the temperature function and is the macroscopic velocity of the ith species. All these variables are space time functions defined on . 0 , , and are correspondingly the diffusivity constant, the valence number and the constant mobility of the ith ion species.

2.1 Initial and Boundary Conditions

For the system (2), the initial conditions at are

 T(x,0)=T0(x),ρi(x,0)=ρi,0(x),i=1,…,N,x∈Ω,

The boundary conditions are of different types: We have homogeneous no-flux boundary conditions for each of the concentrations of ion species,

 ρi→ui⋅→n=0 on ∂Ω.

For the Poisson equation, we partition the boundary into disjoint parts:

 ϕ =δV on ΓD, ϵ∇ϕ⋅→n =E on ΓN, ϵ∇ϕ⋅→n+κϕ =C on ΓR,

For the temperature equation we have Dirichlet boundary conditions which, as we show later help us in deriving the discrete energy law.

 T=δT on ∂Ω.

3 Energy of the Heat-PNP system

According to [19], the Heat-PNP system satisfies simultaneously first and second energy laws of thermodynamics, which are:

 ⎧⎪ ⎪ ⎪ ⎪⎨⎪ ⎪ ⎪ ⎪⎩ddt(U+K)=∫Ω(q+N∑i=0ρi∂ψ∂t)dr+∫∂Ωk∇Tdr,ddtS=Δ+∫ΩqTdr+∫∂Ωk∇TTdr. (3.1)

These laws involve the internal energy functional , the kinetic energy functional of the system, and the entropy functional of the whole system. In addition, is the applied electric field, which, for simplicity, we assume to be independent of time in our system here.

More precisely, if we substitute the predefined entropy functional and entropy production functional for Heat PNP system we obtain:

 ⎧⎪ ⎪ ⎪ ⎪ ⎪ ⎪⎨⎪ ⎪ ⎪ ⎪ ⎪ ⎪⎩S=N∑i=1∫ΩkBρi(logρi−CilogT−Ci)dxΔ=N∑i=1∫Ωνiρi|→ui|2T+1k∣∣∣k∇TT∣∣∣2dx (3.2)

Next, we use the divergence theorem to derive the second law for our system in terms of the unknon variables:

 ddtN∑i=1∫ΩkBρi(logρi−CilogT−Ci)dx=−N∑i=1∫Ωνiρi|→ui|2T+1k|k∇TT|2−∇⋅(k∇TT)dx (3.3)

In fact, to derive the second energy law of thermodynamics from the equation set is straightforward. We multiply by both sides of the heat equation (2) and integrate over the domain. For the left hand side we then we have

 ∫ΩN∑i=1kBCiρi∂T∂t⋅1T+N∑i=1kBCiρi∇⋅(uiT)⋅1T+N∑i=1kB(1−Ci)(ρi)∇⋅uidx = ∫ΩN∑i=1kBCiρi∂logT∂t+∫ΩN∑i=1kBCi(∇(ρiuiT)T−ui∇ρi)+N∑i=1kB(1−Ci)(∇(ρiui−ui∇ρi))dx = ∫ΩN∑i=1kBCiρi∂logT∂t+∫ΩN∑i=1kBCi(∇(ρiui)−ρiui∇TT−ui∇ρi)+N∑i=1kB(1−Ci)(ρi∂t−ρi∂tlogρi)dx = ∫ΩN∑i=1kBCiρi∂logT∂t+∫ΩN∑i=1kBCi(−ρi∂t+ρi∂tlogT−−ρi∂tlogρi)+N∑i=1kB(1−Ci)(ρi∂t−ρi∂tlogρi)dx = ddtN∑i=1∫ΩkBρi(logρi−CilogT−Ci)dx.

Next, for the right hand side we have:

 ∫Ω∇⋅k∇TT+N∑i=1νiρi|ui|2T+qTdx = −N∑i=1∫Ωνiρi|ui|2T+1k|k∇TT|2−∇⋅(k∇TT)

Notice that in the derivation we used the continuity equations. Also, since this energy law is for closed systems, zero flux boundary condition is used when integrating by parts and finally we used that the heat source vanishes, namely, we have .

Clearly, the energy law (3.3) makes sense when ion concentrations and temperature are positive and the solution satisfies the following regularity assumptions:

 ϕ∈H1ΓD≡{ν∈H1(Ω)|ν|ΓD=δV} ρi∈W≡H1∩L∞(Ω)

3.1 Log-density formulation and its energy

We introduce a change of variables (also known as a transformation for the ion concentrations) which is as follows

 ηi(x,t) =logρi(x,t), i=1,...,N ξ(x,t) =logT(x,t)

We note that such change of variables requires that the concentrations are positive and we have

 ηi=logρi∈W≡H1∩L∞(Ω)↪H2

The embedding above holds for space dimensions . The Heat-PNP system then is written in a variational (weak) form: Find , with and such that:

 ⎧⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪⎨⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪⎩(ϵ∇ϕ,∇v)+⟨κϕ,v⟩ΓR−N∑i=1qi(eηi,v)=⟨C,v⟩ΓR+⟨S,v⟩ΓN(∂∂teηi,w)+(∇(eηieξ)+eηizie∇ϕνi,∇w)=0(N∑i=1kBCieηi∂∂teξ,w)−(N∑i=1kBCieξui,∇(eηiw))+(k∇eξ,∇w)=(N∑i=1νieηi|ui|2,w) (3.4)

Here we denote by the macroscopic velocity of ith species defined as , and for the energy law (3.3) we have

 ddtN∑i=1∫ΩkBeηi(ηi−Ciξ−Ci)dx=−N∑i=1∫Ωνieηi|ui|2T+1k|k∇TT|2−∇⋅(k∇TT) (3.5)

Furthermore, if we denote , , and for simplicity we ignore the constant , and under the appropriate zero flux boundary condition we arrive at the following simplified form of (3.3):

 ddtN∑i=1∫Ωeηi(ηi−ξ−1)dx=−N∑i=1∫Ωeηi|ui|2eξ+k|∇ξ2|. (3.6)

The derivation of the the energy law (3.6) from the variational form (3.4) is as follows: We choose the test function to be in the temperature equatio and we obtain

 (N∑i=1eηi∂∂teξ,e−ξ)−(N∑i=1eξui,∇(eηie−ξ))+(k∇eξ,∇e−ξ) =(N∑i=1eηi|ui|2,e−ξ) ∫ΩN∑i=1eηi∂ξ∂t−(eξ(∇ηi+∇ξ)+zieϕ,e−ξeηi(∇ηi−∇ξ))−k|∇ξ2| =N∑i=1∫Ωeηi|ui|2eξ ∫ΩN∑i=1eηi∂ξ∂t−(eξeηi(∇ηi+∇ξ)+ezieηiϕ,(∇ηi−∇ξ))−k|∇ξ2| =N∑i=1∫Ωeηi|ui|2eξ

Next, we test the continuity equation with :

 (∂∂teηi,ηi−ξ)+(∇(eηieξ)+eηizieϕ,∇(ηi−ξ)) =0 (eηi∂∂tηi,ηi−ξ) =−(eηieξ∇(ηi+ξ)+eηizieϕ,∇(ηi−ξ))

Combining the temperature and the continuity equations together shows that:

 ∫ΩN∑i=1eηi∂ξ∂t+(eηi∂∂tηi,ηi−ξ) =−N∑i=1∫Ωeηi|ui|2eξ+k|∇ξ2| ∫ΩN∑i=1ddteηi(ηi−ξ−1) =−N∑i=1∫Ωeηi|ui|2eξ+k|∇ξ2|

3.2 The discrete formulation

We consider a mesh of simplices covering our computational domain (triangles in2D or tetrahedra in 3D). As approximating space we take the space for piecewise linear, with respect to , continuous polynomials [2],

 Wh≡{wh∈H1|wh|τ∈P1 for % all τ in Th}⊂H1
 Vh,ΓD≡{vh∈Wh|v|ΓD=h|ΓD}

The finite element solution to the Heat-PNP system then is defined using these finite element spaces. As a time marching scheme we choose the backward Euler scheme which is implicit and therefore stable. The discrete variational form then is: Find and

 ⎧⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪⎨⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪⎩(ϵ∇ϕjh,∇vh)+⟨κϕjh,vh⟩ΓR−N∑i=1qi(eηji,h,vh)=⟨C,vh⟩ΓR+⟨S,vh⟩ΓN1Δtj(eηji,h,wh)+(∇(eηji,heξjh+qiϕjh)νi,∇wh)=1Δtj(eηj−1i,h,wh)1Δtj(N∑i=1kBCieηji,heξjhξjh,wh)−(N∑i=1kBCieξjhuji,h,∇(ηji,h,wh))+(k∇eξjh,∇wh)=(N∑i=1νieηji,h|uji,h|2,wh)+1Δtj(∑Ni=1kBCieηj−1i,heξjhξj−1h,wh) (3.7)

The approximation to the initial conditions uses the standard interpolation operator

and is as follows:

 ηi(x,0)=Ih(log(ρi,0(x))) for x ∈Ω,i = 1,...,N ξ(x,0)=Ih(log(T0(x))) for x ∈Ω,

3.3 A discrete energy estimate

Next result shows a discrete energy estimate which holds in case when the mesh is quasi-uniform mesh and the stiffness matrix on each Picard iteration is an -matrix. In this case, we can show that the corresponding discrete solution satisfies maximum principle and with this condition satisfied, we have discrete energy estimate and mass conservation.

Theorem 1.

Suppose and satisfy equations (3.7) for and the assumptions for the discrete maximum principle of temperature equation are satisfied. Then the mass is conserved for each ion species

 ∫Ωeηji,h(x,t)dx=∫Ωeη0i,h(x,t)dx, for i=1,..,N,j=1,...m. (3.8)

Moreover, the discrete analogue of the energy estimate (second law of thermodynamics) holds

 N∑i=1∫Ωeηji,h(ηji,h−ξjh−1)dx +m∑j=1ΔtjN∑i=1∫Ω−eηji,h|ui|2eξjh+|k∇eξjh|2⋅(1+Mh12) ≤ N∑i=1∫Ωeη0i,h(η0i,h−ξ0h−1)dx

where and h is the mesh size.

Proof.

For simplicity, let , , , , and , and the corresponding discrete energy estimate is done as follows. First, we choose in the first equation of (3.7) and this gives us:

 1Δtj∫Ωeηji,h−eηj−1i,h=0

For the energy estimate, we test the last equation in (3.7) with . The latter is a valid test function because and is the interpolation operator mapping to to the piecewise linear, continuous .This gives us:

 1Δtj(N∑i=1eηj+1i,heξjhξj+1h,Ih(e−ξjh))−1Δtj(N∑i=1eηji,heξjhξjh,Ih(e−ξjh))I −(N∑i=1eξjhuji,h,∇(eηji,hIh(e−ξjh)))II +(k∇eξjh,∇Ih(e−ξjh))III = (N∑i=1νieηji,h|uji,h|2,1eξjh)IV

We analyze each term in the above equation sum and we begin by rewriting I on the left hand side. We have

 I = 1Δtj(N∑i=1eηj+1i,heξjhξj+1h,Ih(e−ξjh))−1Δtj(N∑i=1eηji,heξjhξjh,Ih(e−ξjh)) = 1Δtj((N∑i=1eηj+1i,h,ξj+1h)−(N∑i=1eηji,h,ξjh))(Ih(e−ξjh)eξjh)

Next, we rewrite also II using integration by parts and the zero flux boundary conditions:

 II = −(N∑i=1eξjhuji,h,∇(eηji,hIh(e−ξjh))) = (N∑i=1∇(eξjhuji,h)eηji,hIh(e−ξjh)) = (Ih(e−ξjh)eξjh)⋅((N∑i=1uji,heηji,h,∇ξjh)+(N∑i=1∇uji,heηji,h,1)) = (Ih(e−ξjh)eξjh)⋅(1Δtj(eηj+1i,h,ηj+1i,h−ξjh−1)−1Δtj(eηji,h,ηji,h−ξjh−1))

In the last step we used the discretized continuity equation to estimate term III on the left hand side as follows:

 III = (k∇eξjh,∇Ih(e−ξjh)) = k|∇eξjh|2−(k∇eξjh,∇Ih(e−ξjh)−∇e−ξjh) ≤ k|∇eξjh|2+k∥∇eξjh∥2⋅∥∇Ih(e−ξjh)−∇e−ξjh∥2

Here we have used the assumption that , for any and that these constants are uniformly bounded by the constant . As shown in [5, Theorem 3.1.6] we have the following interpolation estimate

 ∥∇Ih(e−ξjh)−∇e−ξjh∥L2 ≤Ch|e−ξjh|2,Ω=C(Ω)hM.

term4 on the right hand side:

 (N∑i=1eηji,h|uji,h|2,1eξjh)=⎛⎜ ⎜⎝∑Ni=1eηji,h|uji,h|2eξjh⎞⎟ ⎟⎠(Ih(e−ξjh)eξjh)

Combining all the terms above and by a discrete Grnwall argument of telescopic sum from time step t=0 to m and, as a result, we obtain the following energy estimate:

 m∑j=1Δtj(I+II) ≥ m∑j=1Δtj(IV−(Ch|Ih(e−ξjh)−e−ξjh|2,Ω⋅∥∇eξjh∥2))

only keep the term in I to the right side, we have:

 N∑i=1∫Ωeηji,h(ηji,h−ξjh−1)dx + m∑j=1ΔtjN