*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 \ \).

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

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

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

## 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:

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

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

where

After simplifying the equations and denoting

the dimensionless form of the PNP system of equations is

The Nernst-Planck equation boundary condition has the form

The weak form of the Nernst-Planck equation is

and the weak form of the Poisson equation is

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

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

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,

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

## 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:

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.

## 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

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:

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:

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.

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.

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`