Navigation

Nernst-Planck Equation System

By David Pugal, University of Nevada, Reno.

Problem description

Equation reference: The first version of the following derivation was published in: IPMC: recent progress in modeling, manufacturing, and new applications D. Pugal, S. J. Kim, K. J. Kim, and K. K. Leang Proc. SPIE 7642, (2010). The following Bibtex entry can be used for the reference:

@conference{pugal:76420U,
        author = {D. Pugal and S. J. Kim and K. J. Kim and K. K. Leang},
        editor = {Yoseph Bar-Cohen},
        title = {IPMC: recent progress in modeling, manufacturing, and new applications},
        publisher = {SPIE},
        year = {2010},
        journal = {Electroactive Polymer Actuators and Devices (EAPAD) 2010},
        volume = {7642},
        number = {1},
        numpages = {10},
        pages = {76420U},
        location = {San Diego, CA, USA},
        url = {http://link.aip.org/link/?PSI/7642/76420U/1},
        doi = {10.1117/12.848281}
}

The example is concerned with the finite element solution of the Poisson and Nernst-Planck equation system. The Nernst-Planck equation is often used to describe the diffusion, convection, and migration of charged particles:

(1)\frac {\partial C} {\partial t} + \nabla \cdot (-D\nabla C - z \mu F C \nabla \phi) = - \vec {u} \cdot \nabla C.

The second term on the left side is diffusion and the third term is the migration that is directly related to the the local voltage (often externally applied) \phi. The term on the right side is convection. This is not considered in the current example. The variable C is the concentration of the particles at any point of a domain and this is the unknown of the equation.

One application for the equation is to calculate charge configuration in ionic polymer transducers. Ionic polymer-metal composite is for instance an electromechanical actuator which is basically a thin polymer sheet that is coated with precious metal electrodes on both sides. The polymer contains fixed anions and mobile cations such as H^{+}, Na^{+} along with some kind of solvent, most often water.

When an voltage V is applied to the electrodes, the mobile cations start to migrate whereas immobile anions remain attached to the polymer backbone. This creates spatial charges, especially near the electrodes. One way to describe this system is to solve Nernst-Planck equation for mobile cations and use Poisson equation to describe the electric field formation inside the polymer. The poisson equation is

(2)\nabla \cdot \vec{E} = \frac{F \cdot \rho}{\varepsilon},

where E could be written as \nabla \phi = - \vec{E} and \rho is charge density, F is the Faraday constant and \varepsilon is dielectric permittivity. The term \rho could be written as:

(3)\rho = C - C_{anion},

where C_{anion} is a constant and equals anion concentration. Apparently for IPMC, the initial spatial concentration of anions and cations are equal. The inital configuration is shown:

Initial configuration of IPMC.

The purploe dots are mobile cations. When a voltage is applied, the anions drift:

Bent IPMC

Images reference: IPMC: recent progress in modeling, manufacturing, and new applications D. Pugal, S. J. Kim, K. J. Kim, and K. K. Leang Proc. SPIE 7642, (2010) This eventually results in actuation (mostly bending) of the material (not considered in this section).

To solve equations (1) and (2) boundary conditions must be specified as well. When solving in 2D, just a cross section is considered. The boundaries are shown in:

IPMC boundaries

For Nernst-Planck equation (1), all the boundaries have the same, insulation boundary conditions:

(4)-D \frac{\partial C}{\partial n} - z \mu F C \frac{\partial \phi} {\partial n} = 0

For Poisson equation:

  1. (positive voltage): Dirichlet boundary \phi = 1V. For some cases it might be necessary to use electric field strength as the boundary condtition. Then the Neumann boundary \frac{\partial \phi}{\partial n} = E_{field} can be used.
  2. (ground): Dirichlet boundary \phi = 0.
  3. (insulation): Neumann boundary \frac{\partial \phi}{\partial n} = 0.

Weak form of the equations

To implement the (1) and (2) in Hermes2D, the weak form must be derived. First of all let’s denote:

  • K=z \mu F
  • L=\frac{F}{\varepsilon}

So equations (1) and (2) can be written:

(5)\frac{\partial C}{\partial t}-D\Delta C-K\nabla\cdot \left(C\nabla\phi\right)=0,

(6)-\Delta\phi=L\left(C-C_{0}\right),

Then the boundary condition (4) becomes

(7)-D\frac{\partial C}{\partial n}-KC\frac{\partial\phi}{\partial n}=0.

Weak form of equation (5) is:

(8)\int_{\Omega}\frac{\partial C}{\partial t}v d\mathbf{x} -\int_{\Omega}D\Delta Cv d\mathbf{x}-\int_{\Omega}K\nabla C\cdot \nabla\phi v d\mathbf{x} - \int_{\Omega}KC\Delta \phi v d\mathbf{x}=0,

where v is a test function \Omega\subset\mathbf{R}^{3}. When applying Green’s first identity to expand the terms that contain Laplacian and adding the boundary condition (7), the (8) becomes:

(9)\int_{\Omega}\frac{\partial C}{\partial t}v d\mathbf{x}+ D\int_{\Omega}\nabla C\cdot\nabla v d\mathbf{x}- K\int_{\Omega}\nabla C \cdot \nabla \phi v d\mathbf{x}+ K\int_{\Omega}\nabla\left(Cv\right)\cdot \nabla \phi d\mathbf{x}- D\int_{\Gamma}\frac{\partial C}{\partial n}v d\mathbf{S}- \int_{\Gamma}K\frac{\partial\phi}{\partial n}Cv d\mathbf{S}=0,

where the terms 5 and 6 equal 0 due to the boundary condition. By expanding the nonlinear 4th term, the weak form becomes:

(10)\int_{\Omega}\frac{\partial C}{\partial t}v d\mathbf{x}+ D\int_{\Omega}\nabla C \cdot \nabla v d\mathbf{x}- K\int_{\Omega}\nabla C \cdot \nabla \phi v d\mathbf{x}+ K\int_{\Omega}\nabla \phi \cdot \nabla C v d\mathbf{x}+ K\int_{\Omega} C \left(\nabla\phi\cdot\nabla v\right) d\mathbf{x}=0

As the terms 3 and 4 are equal and cancel out, the final weak form of equation (5) is

(11)\int_{\Omega}\frac{\partial C}{\partial t}v d\mathbf{x}+ D\int_{\Omega}\nabla C \cdot \nabla v d\mathbf{x}+ K\int_{\Omega} C \left(\nabla\phi\cdot\nabla v\right) d\mathbf{x}=0

The weak form of equation (6) with test function u is:

(12)-\int_{\Omega}\Delta\phi u d\mathbf{x}-\int_{\Omega}LCu d\mathbf{x}+ \int_{\Omega}LC_{0}u d\mathbf{x}=0.

After expanding the Laplace’ terms, the equation becomes:

(13)\int_{\Omega}\nabla\phi\cdot\nabla u d\mathbf{x}-\int_{\Omega}LCu d\mathbf{x}+ \int_{\Omega}LC_{0}u d\mathbf{x}=0.

Notice, when electric field strength is used as a boundary condition, then the contribution of the corresponding surface integral must be added:

(14)\int_{\Omega}\nabla\phi\cdot\nabla u d\mathbf{x}-\int_{\Omega}LCu d\mathbf{x}+ \int_{\Omega}LC_{0}u d\mathbf{x}+\int_{\Gamma}\frac{\partial \phi}{\partial n}u d\mathbf{S}=0.

However, for the most cases we use only Poisson boundary conditions to set the voltage. Therefore the last term of (14) is omitted and (13) is used instead in the following sections.

Jacobian matrix

Equation (10) is time dependent, thus some time stepping method must be chosen. For simplicity we start with first order Euler implicit method

(15)\frac{\partial C}{\partial t} \approx \frac{C^{n+1} - C^n}{\tau}

where \tau is the time step. We will use the following notation:

(16)C^{n+1} = \sum_{k=1}^{N^C} y_k^{C} v_k^{C}, \ \ \ \phi^{n+1} = \sum_{k=1}^{N^{\phi}} y_k^{\phi} v_k^{\phi}.

In the new notation, time-discretized equation (11) becomes:

(17)F_i^C(Y) = \int_{\Omega} \frac{C^{n+1}}{\tau}v_i^C d\mathbf{x} - \int_{\Omega} \frac{C^{n}}{\tau}v_i^C d\mathbf{x} + D\int_{\Omega} \nabla C^{n+1} \cdot \nabla v_i^C d\mathbf{x} + K \int_{\Omega}C^{n+1} (\nabla \phi^{n+1} \cdot \nabla v_i^C) d\mathbf{x},

and equation (13) becomes:

(18)F_i^{\phi}(Y) = \int_{\Omega} \nabla \phi^{n+1} \cdot \nabla v_i^{\phi} d\mathbf{x} - \int_{\Omega} LC^{n+1}v_i^{\phi} d\mathbf{x} + \int_{\Omega} LC_0 v_i^{\phi} d\mathbf{x}.

The Jacobian matrix DF/DY has 2\times 2 block structure, with blocks corresponding to

(19)\frac{\partial F_i^C}{\partial y_j^C}, \ \ \ \frac{\partial F_i^C}{\partial y_j^{\phi}}, \ \ \ \frac{\partial F_i^{\phi}}{\partial y_j^C}, \ \ \ \frac{\partial F_i^{\phi}}{\partial y_j^{\phi}}.

Taking the derivatives of F^C_i with respect to y_j^C and y_j^{\phi}, we get

(20)\frac{\partial F_i^C}{\partial y_j^C} = \int_{\Omega} \frac{1}{\tau} v_j^C v_i^C d\mathbf{x} + D\int_{\Omega} \nabla v_j^C \cdot \nabla v_i^C d\mathbf{x} + K\int_{\Omega} v_j^C (\nabla \phi^{n+1} \cdot \nabla v_i^C) d\mathbf{x},

(21)\frac{\partial F_i^C}{\partial y_j^{\phi}} = K \int_{\Omega} C^{n+1} (\nabla v_j^{\phi} \cdot \nabla v_i^C) d\mathbf{x}.

Taking the derivatives of F^{\phi}_i with respect to y_j^C and y_j^{\phi}, we get

(22)\frac{\partial F_i^{\phi}}{\partial y_j^C} = - \int_{\Omega} L v_j^C v_i^{\phi} d\mathbf{x},

(23)\frac{\partial F_i^{\phi}}{\partial y_j^{\phi}} = \int_{\Omega} \nabla v_j^{\phi} \cdot \nabla v_i^{\phi} d\mathbf{x}.

In Hermes, equations (17) and (18) are used to define the residuum F, and equations (20) - (23) to define the Jacobian matrix J. It must be noted that in addition to the implicit Euler iteration Crank-Nicolson iteration is implemented in the code (see the next section for the references of the source files).

Meshing and Multimesh

When it comes to meshing, it should be considered that the gradient of C near the boundaries will be higher than gradients of \phi. This allows us to create different meshes for those variables. In main.cpp. the following code in the main() function enables multimeshing

// Spaces for concentration and the voltage.
H1Space C(&Cmesh, C_bc_types, C_essential_bc_values, P_INIT);
H1Space phi(MULTIMESH ? &phimesh : &Cmesh, phi_bc_types, phi_essential_bc_values, P_INIT);

When MULTIMESH is defined in main.cpp. then different H1Spaces for phi and C are created. It must be noted that when adaptivity is not used, the multimeshing in this example does not have any advantage, however, when adaptivity is turned on, then mesh for H1Space C is refined much more than for phi.

Sample results (non-adaptive)

The following figure shows the calculated concentration C inside the IPMC.

Calculated concentration

As it can be seen, the concentration is rather uniform in the middle of domain. In fact, most of the concentration gradient is near the electrodes, within 5...10% of the total thickness. That is why the refinement

The voltage inside the IPMC forms as follows:

Calculated voltage inside the IPMC

Here we see that the voltage gradient is smaller and more uniform near the boundaries than it is for C. That is where the adaptive multimeshing can become useful.

Sample results (adaptive)

To be added soon.

Table Of Contents

Previous topic

Nernst-Planck

Next topic

Nernst-Planck-IPMC

Navigation

© Copyright 2009-2011, hp-FEM group at UNR. Last updated on Sep 15, 2014. Created using Sphinx 1.2.3.