hp-FEM: Poisson-Nernst-Planck system

The PNP system is highly nonlinear and for a typical domain with two electrodes, largest differences in charge concentration occur in a very narrow region near the boundary. The computing power required for a full scale problem is significant.

We consider a rectangular 2D domain \(\Omega\subset\mathbb{R}^2 \ \) with boundaries \(\partial\Omega_{1\ldots 4}\subset\partial\Omega \ \).

Calculation domain

The Poisson-Nernst-Planck (PNP) system of equations (see Electromechanical transduction section).

As there is no flow through the domain’s boundary, the Nernst-Planck equation is equipped with a Neumann boundary condition

\[-D \frac{\partial C}{\partial n} - \mu F C \frac{\partial \phi} {\partial n} = 0.\]

Furthermore, we prescribe a positive constant voltage \(V_{pos} \ \) on \(\Omega_1 \ \) and zero voltage on \(\Omega_3 \ \):

\[\begin{split}\phi_{\partial\Omega_1}&=&V_{pos},\\ \phi_{\partial\Omega_3}&=&0.\end{split}\]

On the rest of the boundary, \(\phi \ \) has zero normal derivatives, and thus we prescribe a Neumann boundary condition

\[\frac{\partial \phi_{\Omega_2}}{\partial n}=\frac{\partial \phi_{\Omega_4}}{\partial n}=0.\]

Dimensionless weak form of the PNP system

To make our results easily reproducible, weak form of the PNP system of equations as well as the Jacobian matrix and residual vector that are used in actual computations are presented. To simplify notation, dimensionless formulation is derived. The following notations for the independent variables \(x \ \), \(y \ \), \(t \ \) and for the dependent variables \(C \ \) and \(\phi \ \) are used:

\[X=\frac xl,\ \ \ Y=\frac yl,\ \ \ \tau=\frac{tD}{\lambda_D l},\ \ \ \varphi=\frac{\phi F}{RT},\ \ \ c=\frac C C_0.\]

Here \(\lambda_D \ \) is the Debye screening length and it is expressed as follows:

\[\lambda_D=\sqrt{\frac{\varepsilon R T}{2 F^2 C_0}}.\]

After inserting the variables into the Poisson and Nernst-Planck equations, the system becomes:

\[\begin{split}\frac{DC_0}{\lambda_D l}\frac{\partial c}{\partial \tau} + \frac 1l \nabla_d\cdot \left(-\frac{DC_0}{l}\nabla_dc-c\frac{DC_0}{l}\nabla_d\varphi\right) &=& 0\label{eq:dim1}\\ - \frac{\varepsilon R T}{l^2F^2C_0}\nabla_d^2\varphi&=&c-1,\end{split}\]

where

\[\nabla_d=\left(\frac{\partial}{\partial X},\ \frac{\partial}{\partial Y}\right).\]

After simplifying the equations and denoting

\[\epsilon=\frac{\lambda_D}{l},\]

the dimensionless form of the PNP system of equations is

\[\begin{split}\frac{\partial c}{\partial \tau} -\epsilon\nabla_d^2 c - \epsilon \nabla_d \cdot \left( c\nabla_d \varphi\right) &=&0,\label{eq:dimensionless-nernst-planck}\\ -\nabla_d^2\varphi&=&\frac{c-1}{2\epsilon^2}.\end{split}\]

The Nernst-Planck equation boundary condition has the form

\[-\frac{\partial c}{\partial n}-c\frac{\partial\varphi}{\partial n}=0.\]

The weak form of the Nernst-Planck equation is

\[\int_{\Omega}\frac{\partial c}{\partial \tau}v^c d\mathbf{x}+ \epsilon\int_{\Omega}\nabla_d c\cdot\nabla_d v^c d\mathbf{x}+ \epsilon \int_{\Omega} c \left(\nabla_d\varphi\cdot\nabla_d v^c\right) d\mathbf{x}=0\]

and the weak form of the Poisson equation is

\[\int_{\Omega}\nabla_d\varphi\cdot\nabla_d v^\varphi d\mathbf{x}- \frac{1}{2\epsilon^2}\left[\int_{\Omega}cv^\varphi d\mathbf{x}- \int_{\Omega}v^\varphi d\mathbf{x}\right]=0.\]

The time discretized first part \(F^c \ \) of the residual vector \(F \ \) is

\[\begin{split}F_i^c\left(Y\right) & = & \int_{\Omega} \frac{c(Y)}{\delta\tau}v_i^c d\mathbf{x} - \int_{\Omega} \frac{c^{n}}{\delta\tau}v_i^c d\mathbf{x}\nonumber\\ &&+\frac 12 \epsilon \left[\int_{\Omega} \nabla_d c(Y) \cdot \nabla_d v_i^c d\mathbf{x}+ \int_{\Omega} \nabla_d c^{n} \cdot \nabla_d v_i^c d\mathbf{x}\right.\nonumber\\ &&+ \left.\int_{\Omega}c(Y) \left(\nabla_d \varphi(Y) \cdot \nabla_d v_i^c\right) d\mathbf{x}+ \int_{\Omega}c^{n} \left(\nabla_d \varphi^{n} \cdot \nabla_d v_i^c\right) d\mathbf{x}\right]\label{eq:Fc}\end{split}\]

where \(i = 1, 2, \ldots, N^c \ \). Analogously, the second part \(F^{\varphi} \ \) of the residual vector \(F \ \) can be written

\[F_i^{\varphi}\left(Y\right) = \int_{\Omega} \nabla_d \varphi(Y) \cdot \nabla_d v_i^{\varphi} d\mathbf{x} -\frac{1}{2\epsilon^2}\left[ \int_{\Omega} c(Y)v_i^{\varphi} d\mathbf{x} - \int_{\Omega} v_i^{\varphi} d\mathbf{x}\right]\]

where \(i = N^c + 1, N^c + 2, \ldots, N^c + N^{\varphi} \ \). The nonlinear discrete problem that needs to be solved at the end of each time step thus has the form \(F(Y) = 0 \ \).

The Jacobian matrix \(J(Y) = DF/DY \ \) has a \(2\times 2 \ \) block structure,

\[\begin{split}J(Y) = \left( \begin{array}{cc} \displaystyle \frac{\partial F_i^c}{\partial y_j^c} & \displaystyle \frac{\partial F_i^c}{\partial y_j^{\varphi}} \\ \displaystyle \frac{\partial F_i^{\varphi}}{\partial y_j^c} & \displaystyle \frac{\partial F_i^{\varphi}}{\partial y_j^{\varphi}} \end{array}\right),\end{split}\]

and its entries are obtained by calculating the partial derivatives of \(F \ \) with respect to the components of the coefficient vector \(Y \ \):

\[\begin{split}\frac{\partial F_i^c}{\partial y_j^c}(Y) &=& \int_{\Omega} \frac{1}{\delta\tau} v_j^c v_i^c d\mathbf{x} + \frac 12 \epsilon\left[\int_{\Omega} \nabla_d v_j^c \cdot \nabla_d v_i^c d\mathbf{x} + \int_{\Omega} v_j^c \left(\nabla_d \varphi(Y) \cdot \nabla_d v_i^c\right) d\mathbf{x}\right],\label{eq:dFcdyc}\\ \frac{\partial F_i^c}{\partial y_j^{\varphi}}(Y) &=& \frac 12 \epsilon \int_{\Omega} c(Y) \left(\nabla_d v_j^{\varphi} \cdot \nabla_d v_i^c\right) d\mathbf{x},\\ \frac{\partial F_i^{\varphi}}{\partial y_j^c}(Y) &=& - \frac{1}{2\epsilon^2}\int_{\Omega} v_j^c v_i^{\varphi} d\mathbf{x},\\ \frac{\partial F_i^{\varphi}}{\partial y_j^{\varphi}}(Y) &=& \int_{\Omega} \nabla_d v_j^{\varphi} \cdot \nabla_d v_i^{\varphi} d\mathbf{x}.\end{split}\]

Selected calculation results

Running the simulation with different adaptivity modes and meshes showed that the multi-mesh hp-FEM configuration resulted in the smallest problems and similar error convergence compared to any single-mesh configuration. However, multi-mesh problems generally resulted in somewhat longer computing times. The problem size is illustrated for HP_ANISO adaptivity mode:

Number of DOFs as a function of physical time for single-mesh and multi-mesh configurations with HP ANISO adaptivity mode.

Number of DOFs as a function of physical time for single-mesh and multi-mesh configurations with HP ANISO adaptivity mode. Copyright 2012 Global Science Press.

The role of anisotropic mesh refinements is illustrated in the following figure. Typical results for the HP_ISO, HP_ANISO_H, HP_ANISO adaptivity modes in terms of DOF are shown. It can be seen that HP ISO is notoriously inefficient.

Number of DOF as function of physical time for multi-mesh configuration with HP\_ANISO, HP\_ANISO\_H, and HP\_ISO adaptivity modes (log y scale).

Number of DOF as function of physical time for multi-mesh configuration with HP_ANISO, HP_ANISO_H, and HP_ISO adaptivity modes (log y scale). Copyright 2012 Global Science Press.

Calculations with realistic BCs

In real physics calculations, the applied voltage on boundary \(\partial\Omega_1 \ \) is not constant. This can be, for instance, due to the high resistance of the electrodes. To see how the i HP_ANISO adaptivity works for such situations, the voltage on the boundary was applied as follows

\[\phi_{\Omega_1}\left( x \right)=0.5\left[V \right] \frac{x\left[ m \right ]}{\text{width}_{\Omega_1}\left[ m \right]}+0.5\left[ V \right],\]

where \(\text{width}_{\Omega_1} \ \) is the width of the boundary. The given boundary is effectively a linear increase of the voltage from \(\phi_{\Omega_1}\left(x = 0 \right)=0.5\ \mathrm{V} \ \) to \(\phi_{\Omega_1}\left(x=\text{width}_{\Omega_1}\right) = 1.0\ \mathrm{V} \ \). Now the concentration gradient \(\nabla c \ \) and the voltage gradient \(\nabla \varphi \ \) are no longer effectively in 1D. The calculated scaled values \(c \ \) and \(\varphi \ \) in \(\Omega \ \) and corresponding meshes and polynomial degrees of the elements at \(t=0.1\ \mathrm{s} \ \) are shown:

Solutions :math:`c \ ` and :math:`\varphi \ ` and corresponding polynomial degrees of the elements at t = 0.1s. HP\_ANISO refinement mode was used. The height in the solution graphs indicates the value.

Solutions \(c \ \) and \(\varphi \ \) and corresponding polynomial degrees of the elements at t = 0.1s. HP_ANISO refinement mode was used. The height in the solution graphs indicates the value. Copyright 2012 Global Science Press.

The HP_ANISO adaptivity algorithm has particularly increased the polynomial degree and refined the mesh near \(\Omega_1 \ \) where a sharp concentration peak exists. This example clearly illustrates how the solution of PNP with non-uniform boundary conditions is very dynamic in time and how the HP_ANISO time dependent adaptivity finds an optimal mesh and polynomial space to adapt to the dynamics of the problem.

Length scale analysis

The Debye length \(\lambda_D \ \) is the screening length in the electrolyte solutions. Its numerical value shows the thickness of the charged layer in the vicinity of the boundaries \(\partial\Omega_1 \ \) and \(\partial\Omega_2 \ \). In all the previous simulations, the Debye screening length was determined by the constants and was calculated to be \(\lambda_D=1.7\ \mu m \ \). It is known that computation gets increasingly difficult when reducing the value of \(\lambda_D \ \). It was our interest to see how small screening lengths can Hermes HP_ANISO automatic adaptivity handle. The parameter \(\varepsilon \ \) was varied as follows:

\[\varepsilon_n=\varepsilon\times 0.5^n,\]

The simulations were run for each \(\varepsilon_n \ \) and corresponding \(\lambda_D \ \) value and maximum number of degrees of freedom and cumulative CPU time were recorded. The simulation time \(t \ \) for each \(\lambda_D \ \) was chosen to be \(\tau \ \) — the characteristic time scale — and each simulation was divided equally into fifteen time steps. The following figure shows the maximum and average number of degrees of freedom during calculation as a function of the Debye length.

Problem size depending on the Debye length. Initially the coarse mesh was used. It was necessary to use the fine mesh for smaller :math:`\lambda_D \ ` values.

Problem size depending on the Debye length. Initially the coarse mesh was used. It was necessary to use the fine mesh for smaller \(\lambda_D \ \) values. Copyright 2012 Global Science Press.

The simulations up to \(0.52\ \mathrm{nm} \ \) screening length were carried out on the initial coarse mesh. However, from \(\lambda_D > 0.52\ \mathrm{nm} \ \), the finer initial mesh had to be used so the existence of the large gradients of the physical fields \(c \ \) and \(\varphi \ \) near the boundaries could be captured in the first place.

The fine mesh allowed simulations with the Debye length down to \(0.40\ \mathrm{nm} \ \). The calculated \(c \ \) and \(\varphi \ \) at \(t=\tau \ \) for \(\lambda_D=0.40\ \mathrm{nm} \ \) are shown in the following figure.

Calculated fields :math:`c \ ` and :math:`\varphi \ ` at :math:`t=\tau=0.81\ ms \ ` for :math:`\lambda_D=0.4\ nm. \ `

Calculated fields \(c \ \) and \(\varphi \ \) at \(t=\tau=0.81\ \mathrm{ms} \ \) for \(\lambda_D=0.4\ \mathrm{nm.} \ \) Copyright 2012 Global Science Press.

It appears that when using even finer initial mesh and higher initial polynomial degrees, even smaller Debye lengths could be used when necessary. The polynomial space of \(c \ \) had consistently higher maximum polynomial degree than that of \(\varphi \ \), however, the difference was less noticeable for smaller Debye lengths.

See also

Copyright 2012 Global Science Press. This article may be downloaded for personal use only. Any other use requires prior permission of the author and Global Science Press.

The following article appeared in Communications in Computational Physics and may be found at Vol 11, No 1, page 249-270

PDF file of the full article Modeling Ionic Polymer-Metal Composites with Space-Time Adaptive Multimesh hp-FEM can be downloaded here