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):
where Ni is the total convective and diffusive flux given by:
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:
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":
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.]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
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.]
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.
To regenerate this electromigration result, change the anode and cathode boundary condition temperatures to 298 in Np-case.sif (above).
This simulation uses the full Np-case.sif input file as given above.
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.