From Elmer Wiki
Jump to: navigation, search

Until 2010, only expensive FEA suites like COMSOL, or specialized consultants, could do electrochemistry simulations. Full treatment of this topic requires coupled solutions of the Navier-Stokes equations with diffusion and electromigration (AdvectionDiffusion model) for convective mass transfer. With the addition of electromigration, Elmer now has this capability.


The Advection Diffusion equation is given in the Elmer Models Manual (using bold instead of arrows for vectors):

∂Ci/∂t = -∇⋅Ni + Si,

where Ni is the total convective and diffusive flux given by:

Ni = -v Ci - DiCi - Mi Ciφ.

The third term in the flux is the Nernst-Planck electromigration term using mobility Mi for species i and the gradient in electrical potential φ. The mobility Mi is in turn given by the Nernst-Einstein relation:

Mi = Di zi e/kT,

where zi is the ion charge for species i, e the electron charge, k Boltzmann's constant, and T the absolute temperature.


A new addition to the AdvectionDiffusion module implements this Nernst-Planck electromigration term. The implementation is based on the similarity between the electromigration flux and the convection flux v Ci, it adds the Nernst-Planck ion drift velocity to create an "effective velocity":

Ni = -DiCi + (v + Di zi eφ/kT) Ci.

The downside of this implementation is that because neither convection nor electromigration is in the variational form, flux boundary conditions ignore them, so it is tricky to set boundary conditions properly. Then again, this is generally true for electrochemical reactions... More on that another time.

This patch introduces only one new parameter or material property: the ion charge zi, the others are already in the AdvectionDiffusion module or Elmer's universal constants. If ion charge is zero, then this term vanishes. Ion charge is now in the material properties section of the advection-diffusion GUI XML file, and to the solver keywords.


This code does not create a velocity array, so one must use Navier-Stokes when modeling electromigration. If one runs with non-zero ion charge to induce electromigration but without Navier-Stokes, ElmerSolver fails with a segmentation fault.


[DISCLAIMER: This example is for purposes of illustration and mathematical consistency verification only. Material properties and results may differ from reality by orders of magnitude. Use at your own risk.]

Nernst-Planck example mesh
Gmsh-generated mesh for Nernst-Planck example
The following example describes electroplating a small hemisphere with copper. This problem has an axisymmetric geometry with a cylindrical container of aqueous copper sulfate (CuSO4) 1 cm in diameter and 1 cm deep, with the 0.5 cm radius hemisphere to be plated pointing downward into the center of the container (along the axis of symmetry). Only the aqueous solution has a mesh, which is densest near the hemisphere in order to resolve the diffusion/electromigration boundary layer there. To generate the mesh, run Gmsh and import the geometry from File:Np-simple.geo.

The model includes:

  • Electrical current using the StatCurrent module, with Nernst-Einstein ionic conductivity σ = σe + ∑ Mi ze NA Ci/MWi (not yet concentration-dependent) where σe is electronic conductivity (assumed zero), NA is Avogadro's number, and MWi is molecular weight (used because Ci is in kg/m3);
  • Copper ion transport using the AdvectionDiffusion module with the electromigration term above;
  • Heat transfer driven by Joule heating of the solution and an optional temperature difference across the cell;
  • (Incompressible) Navier-Stokes flow of the aqueous solution driven by thermal buoyancy, and possibly at some point by magneto-hydrodynamics.

It assumes that charge transfer at the cathode is very fast, so copper plating current is completely limited by mass transfer through the aqueous solution.

The initial condition has zero velocity, uniform temperature of 298 K, and copper ion concentration of 0.1 kg/m3. Boundary conditions are:

  • Cathode (hemisphere/quarter-circle): potential=0, concentration=0, temperature=303 K, no-slip wall.
  • Anode (right edge): potential=0.1, concentration=0.1 kg/m3 (copper anode dissolving due to anodic current), temperature=293 K, no-slip wall.
  • Free surface (top): vy=0, all other BCs natural.
  • Bottom wall: no-slip wall, all other BCs natural.
  • Symmetry axis (left): vx=0, all other BCs natural.

Material properties are mostly those of water:

  • Density: ρ = 1000 kg/m3
  • Heat capacity: cp = 4184 J/kg⋅K
  • Reference temperature: T0 = 298 K
  • Heat expansion coefficient: β = 10-4 K-1
  • Viscosity: μ = 0.001 Pa⋅s
  • Thermal conductivity: k = 0.6 W/m⋅K
  • Copper ion diffusivity: Di = 10-9 m2/s
  • Copper ion mobility: Mi = Di zi e/kT ≈ 8×10-8 m2/V⋅s
  • Electrical (ionic) conductivity: σi = Mi ze NA Ci/MWi, when C=0.1 kg/m3, σ=0.024 (Ωm)-1.

All of these parameters and conditions are given in File:Np-case.sif

Analytical Calculations

Under these conditions, the maximum electric field |E| should be about 0.1 V/5 mm = 20 V/m. With copper mobility of 8x10-8 m2/V⋅s, this produces ion drift velocity Mi |E| around 1.5 μm/sec. If flow velocity is above this, it will disrupt the electromigration/diffusion boundary layer.

At short times (Dt << R2), the boundary layer will grow as the square root of time: δ ~ √Dt. At long times in an infinite domain, the boundary layer will stop growing around the radius size: δ ~ R, but in this finite domain, the concentration will slowly fall to zero.

With convection-diffusion, the boundary layer approaches a thickness given by the ratio Di/vN, where vN is the normal velocity. Using the ion drift velocity calculated above, the boundary layer thickness should be on the order of 0.6 mm.

[To do: add Grashof/Rayleigh number calculations to describe the nature of natural convection flow and estimate velocities.]

Numerical Results

Diffusion boundary layer at 7.5, 30, 120 and 480 seconds
Diffusion boundary layer at 7.5, 30, 120 and 480 seconds
With no temperature difference (cathode and anode at 298 K like initial condition) and no electromigration (zero ion charge), this results in simple growth of a diffusion boundary layer from the cathode hemisphere. This trend is shown in the series of figures here: the diffusion boundary layer is twice as thick at 480 seconds as at 120, and twice as thick then as at 30 (boundary layer thickness at 7.5 seconds is not accurate because it is comparable to element thickness). This corresponds roughly to "electroless" plating using a reducing agent in solution instead of electrons at the cathode.

To regenerate the diffusion-only result, remove the line "Concentration Ion Charge = 2" from Np-case.sif (above) and change the anode and cathode boundary condition temperatures to 298.

Electromigration boundary layer and flow velocities at 30, 120 and 480 seconds
Electromigration boundary layer and flow velocities at 30, 120 and 480 seconds
With no temperature difference at the boundaries but with ion charge of 2, the boundary layer grows as before until about 120 seconds, then stops growing at about 1 mm thickness by 480 seconds, as predicted above. Note that the boundary layer at 480 seconds is slightly thicker on the bottom than on the right, because of the higher electric field there. Joule heating produces a small temperature increase, leading to natural convection flow. But maximum vertical velocity is about 15 nm/s is much lower than the 1.5 μm/s ion drift velocity. And with the 0.5-1 cm domain size and the diffusivity here, this gives a Peclet number of 0.075-0.15, indicating that convection is not a significant means of mass transfer.

To regenerate this electromigration result, change the anode and cathode boundary condition temperatures to 298 in Np-case.sif (above).

Concentration boundary layer (contours), flow velocities (arrows) and temperatures (background shading) at 5, 30, 120 and 480 seconds
Concentration boundary layer (contours), flow velocities (arrows) and temperatures (background shading) at 5, 30, 120 and 480 seconds
The hot cathode and cold anode boundary conditions given above lead to a natural convection cell. The velocity reaches a maximum of around 1 mm/second at five seconds of simulation time, then declines as temperature gradients smooth out. At this velocity, convection is a very significant means of mass transfer, and it drags the low-copper region across the top then down in front of the anode, approaching steady state at 480 seconds of simulation time as shown here. This illustrates coupling of diffusion, Nernst-Planck electromigration, and flow driven by natural convection.

This simulation uses the full Np-case.sif input file as given above.

Future Work

It would be helpful for the AdvectionDiffusion module to export the Nernst-Einstein ionic conductivity (tensor field) so one could just provide the electronic conductivity into the StatCurrent module and have it add all of the ionic conductivities. This would let StatCurrent compute the full current accurately for Joule heating and magneto-hydrodynamic forcing.

Also, the example should have a better cathode boundary condition than zero concentration, full Butler-Volmer isn't necessary but at least a linearization would be nice. As it is, ionic conductivity is zero at the boundary, and there is only diffusive flux due to the concentration gradient, but this does not enter the StatCurrent calculations.

In view of this, it might be best for AdvectionDiffusion to export the actual ion current density (vector field), including electromigration and diffusion, for accurate coupling with StatCurrent and other modules.