COMSOL 6.4 - Diffuse Double Layer with Charge Transfer

Diffuse Double Layer with Charge Transfer

In the direct vicinity of an electrode surface (in the range of up to a few nanometers), in the diffuse double layer, the assumption of electroneutrality is not valid due to charge separation. Typically the diffuse double layer may be of interest when modeling very thin layers of electrolyte, for instance in electrochemical capacitors and in atmospheric corrosion problems.

To model the behavior of the diffuse double layer, one needs to solve for the Nernst–Planck equations for all of the ions, in combination with the Poisson’s equation for the potential. The combination of these equations is frequently referred to as the Poisson-Nernst–Planck (PNP) equations.

A problem that arises when modeling the PNP equations is that of how to handle the boundary condition for the potential equation. In this example an assumption of a Stern layer with a constant capacity is used to derive surface charge boundary conditions for Poisson’s equation.

The model reproduces the results of Bazant and Chu (see Ref. 1 and Ref. 2).

See also the Diffuse Double Layer tutorial for an introduction to double layer modeling.

Model Definition

The model geometry is in 1D (a single interval between 0 and L) and consists of one single domain, representing the electrolyte phase, including the diffuse double layer.

Domain equations

The concentrations, ci (SI unit: mol/m3, i=+,-), of two ions of opposite charge
(zi = +1/−1) are solved for in the electrolyte phase. The fluxes (Ji, SI unit: mol/(m2·s)) of these are described by the Nernst–Planck equation

with Di (SI unit: m2/s) being the diffusion coefficient, um,i (SI unit: s·mol/kg) the mobility, F (SI unit: C/mol) Faraday’s constant, and ϕ (SI unit: V) the potential.

Assuming no homogeneous reactions in the electrolyte, the governing equations for the two species become:

For the potential, Poisson’s equation states

where ε is the permittivity (SI unit: F/m) and ρ the charge density (SI unit: C/m3), depending on the ion concentrations according to:

Boundary conditions

The boundaries reside in the reaction plane of the electrodes on each side. The same electrode reaction, in which the positive ion, S+,participates, takes place on both electrodes.

The reaction rate r (SI unit: mol/(m2·s)) is

where Ka and Kc (SI unit: m/s) are the anodic and cathodic rate constants, cM the metal species activity (SI unit: mol/m3, constant) and αa and αc the anodic and cathodic transfer coefficients. ϕΔ (SI unit: V) is the difference in potential between the metal phase, ϕM (SI unit: V), and the reaction plane:

The electrode reaction renders an inward flux for the positive ion according to

on both boundaries. For the negative ion, a zero flux condition is used.

Assuming the reaction plane to be placed at the boundary between the inner (compact) and diffuse double layer, and with the assumption of a Stern compact layer of a constant thickness, λS (SI unit: m), one can derive the following Robin type of boundary condition for the potential:

This condition reduces to a Dirichlet voltage condition for λS = 0, that is, in the absence of a Stern layer. For the case of a nonzero stern layer thickness, the condition can be reformulated as a surface charge condition

Cell potential equation

The problem formulated above can now be solved for given voltages of ϕM in the metal electrode phase for each side. Typically one grounds one electrode and specifies the cell voltage as V so that

However, to solve for a given cell current density, icell (SI unit: A/m2), with V not known a priori, an additional global equation, solving for V, is used, fulfilling the condition:

Global concentration constraint for the negative ion

When solving this system for a stationary solution, the negative ion concentration needs an additional “boot-strap” to render a stable, unique, solution. This is done by adding the following global constraint to the equation system:

where c0 is the initial ion concentration (SI unit: mol/m3), equal for both ions.

The constraint assures that the total number of negative ions is preserved during the iterative solver process. (For time-dependent simulations the constraint can be omitted.)

Dimensionless numbers and parameter values

A number of dimensionless numbers can be derived that govern the behavior of the cell. The problem is solved using a parametric study for a dimensionless parameter εD = (0.001, 0.01, 0.1), defined as

where λD is the Debye length.

The current of the cell is defined via the dimensionless number j = 0.9,

where iD is the Nernstian limiting current density.

The cathodic reaction rate constant is defined using the dimensionless number kc = 10,

The rate of the anodic reaction term is governed by the dimensionless number kr = 10,

and the Stern layer thickness is set using the dimensionless number δ = 0.1,

Results and Discussion

The following dimensionless variables are used when presenting the results:

Figure 1 shows the dimensionless concentration,

. The concentration gradients are steepest close to the electrodes.

Figure 1: Dimensionless concentration profile.

Figure 2 shows the dimensionless charge density profile. Charge separation occurs close to the electrodes. For higher values of εD, the region of charge separation, the diffuse double layer, stretches further into the domain. This is expected since higher εD values effectively mean a shorter domain length.

Figure 2: Dimensionless charge density profile.

Figure 3 shows the potential profile. For higher values of εD the voltage over the cell decreases. This is an expected result since a shorter domain length shortens the potential losses due to ion transport.