API

abstract type AbstractElectrolyteData

Abstract super type for electrolytes.

mutable struct ElectrolyteData{Tγ, Tcache, Texp, Tlog, Tflux} <: AbstractElectrolyteData

Data for electrolyte. It is defined using Base.@kwdef allowing for keyword constructors like

    ElectrolyteData(nc=3,z=[-1,2,1])

The struct has three groups of fields:

• Problem parameters: these are meant to be set by user in order to provide simulation parameters.
• Derived values: these are derived from problem parameters by default, or via update_derived! after setting some of the problems pararmeters
• Fixed values: these are mostly physical constants and values derived from them, they should not be changed.
• Reserved fields: they are used to pass information during simulations and should not be set by the user.

Fields (reserved fields are modified by some algorithms):

• nc::Int64: Number of charged species $N$.

  While in the default constructor, this value is primary to the value of cspecies, update_derived! treats it as a derived datum.

• cspecies::Vector{Int64}: Charged species list. Default: [1...nc]

• na::Int64: Number of surface species

• sspecies::Vector{Int64}: Surface species list. Default: [(nc + 1)...(nc + na)]

• iϕ::Int64: Index of electrostatic potential $ϕ$ in species list.

• ip::Int64: Index of pressure p in species list

• D::Vector{Float64}: Mobility coefficients $D_i\; (i=1…N)$

• z::Vector{Int64}: Charge numbers of ions $z_i\; (i=1…N)$

• M0::Float64: Molar weight of solvent $M_0$

• M::Vector{Float64}: Molar weights of ions $M_i\; (i=1…N)$

• v0::Float64: Molar volume of solvent $v_0$

• v::Vector{Float64}: Molar volumes of ions $v_i\; (i=1…N)$

• κ::Vector{Float64}: Solvation numbers of ions $κ_i\; (i=1…N)$

• actcoeff!::Any: actcoeff!(γ, c, p, ::ElectrolyteData)

  Activity coefficient function. Write activity coefficients $γ_i\; (i=1…N)$ into vector γ. Default: DGML_gamma!.

• c_bulk::Vector{Float64}: Bulk ion concentrations $c_i^b\; (i=1…N)$

• Γ_we::Int64: Working electrode boundary number

• Γ_bulk::Int64: Bulk boundary number

• ϕ_bulk::Float64: Bulk voltage $ϕ^b$

• p_bulk::Float64: Bulk pressure $p^b$

• T::Float64: Temperature $T$

• ε::Float64: Dielectric permittivity of solvent $ε$

• rexp::Any: Regularized exponential, default: exp (unregularized)

• rlog::Any: Regularized logarithm, default: log (unregularized)

• pscale::Float64: Pressure scaling factor. Default: 1.0e9

• eneutral::Bool: Local electroneutrality switch. Default: false

• upwindflux!::Any: Upwind flux caculation method for ionic species. This allows to choose between

  • μex_flux! (default, strongly preferrable): excess chemical potential (SEDAN) scheme
  • act_flux!: scheme based on reciprocal activity coefficients
  • cent_flux!: central scheme

• weights::Vector{Float64}: Species weights for norms in solver control.

• solvepressure::Bool: Solve for pressure.

  This is true by default. Setting this to false can serve two purposes:

  • Use the pressure from the solution of a flow equation
  • Ignore the pressure contribution to the excess chemical potential

• F::Float64: Faraday constant $F$ (fixed)

• ε_0::Float64: Dielectric permittivity of vacuum $ε_0$ (fixed)

• RT::Float64: Molar gas constant scaled with temperature $RT$ (derived)

• Mrel::Vector{Float64}: Solvated molar mass ratio $m_i = \frac{M_i}{M_0} + κ_i \; (i=1… N)$ (derived)

• vrel::Vector{Float64}: Solvated molar volume ratio $\hat v_i = \frac{v_i}{v_0} + κ_i \; (i=1… N)$ (derived)

• barv::Vector{Float64}: Solvated molar volume $\bar v_i = v_i + κ_i v_0 \; (i=1… N)$ (derived)

• tildev::Vector{Float64}: Pressure relevant volume $\tilde v_i = \bar v_i + m_i v_0 \; (i=1… N)$ (derived)

• γ_bulk::Vector{Float64}: Activity coefficients at bulk interface $γ_i^b= γ_i(c_1^b … c_N^b, p^b) \; (i=1… N)$ (derived)

• γk_cache::Any: Cache for activity coefficient calculation (reserved)

• γl_cache::Any: Cache for activity coefficient calculation (reserved)

• ϕ_we::Float64: Working electrode voltage $ϕ_{we}$ (reserved) Used by sweep algorithms to pass boundary data.

• edgevelocity::Union{Float64, Vector{Float64}}: Edge velocity projection (reserved).

• scheme::Symbol: Scheme parameter. Deprecated and disabled. Use upwindflux! instead to change the discretization scheme.

update_derived!(electrolyte::ElectrolyteData)

Update derived electrolyte data. This needs to be called in order to update fields Mrel, vrel, barv, tildev, RT, γ_bulk of ElectrolyteData after changing some of the other parameters.

It also correctes electrolyte.nc in order to be consistent with electrolyte.cspecies.

Called on the passed electrolyte data by PNPSystem, PBSystem, dlcapsweep, ivsweep, cvsweep.

dlcap0(electrolyte)

Return double layer capacitance at zero voltage for symmetric binary electrolyte. The molarity is defined by the bulk concentration value $c_1^b$.

\[ c_{dl}=2 \frac{εε_0F^2c_1^b}{RT}\]

Example

using LessUnitful
ely = ElectrolyteData(c_bulk=fill(0.01ufac"mol/dm^3",2))
round(dlcap0(ely),sigdigits=5) |> u"μF/cm^2"
# output

22.847 μF cm^-2

debyelength(electrolyte)

Return debye length for symmetric binary electrolyte. The molarity is defined by the bulk concentration value $c_1^b$.

\[ l_{debye}=\sqrt{RT\frac{εε_0}{F^2c_1^b}}\]

using LessUnitful
ely = ElectrolyteData(c_bulk=fill(0.01ufac"mol/dm^3",2))
round(debyelength(ely),sigdigits=5) |> u"nm"
# output

4.3018 nm

chargedensity(c,electrolyte)

Calculate charge density from vector of concentrations (in one grid point).

chargedensity(sol,electrolyte)

Calculate charge density vector from solution (on whole grid)

chargedensity(tsol,electrolyte)

Calculate time dependent charge densities from from time/voltage dependent solution solution (on whole grid). Returns a VoronoiFVM.TransientSolution.

 chemical_potential(c, barc, p, barv, electrolyte)

Calculate chemical potential of species with concentration c

\[ μ = \bar v(p-p_{ref}) + RT\log \frac{c}{\bar c}\]

chemical_potentials!(μ,u,electrolyte)

Calculate chemical potentials from concentrations.

Input:

• μ: memory for result (will be filled)
• u: local solution vector (concentrations, potential, pressure)

Returns μ0, μ: chemical potential of solvent and chemical potentials of ions.

using LessUnitful
ely = ElectrolyteData(c_bulk = fill(0.01ufac"mol/dm^3", 2))
μ0, μ = chemical_potentials!([0.0, 0.0], vcat(ely.c_bulk, [0, 0]), ely)
round(μ0, sigdigits = 5), round.(μ, sigdigits = 5)
# output

(-0.89834, [-21359.0, -21359.0])

rrate(R0,β,A)
rrate(R0,β,A, electrolyte)

Reaction rate expression

rrate(R0,β,A)=R0*(exp(-β*A) - exp((1-β)*A))

iselectroneutral(cx::Vector,electrolyte)

Check for electroneutrality of concentration vector

iselectroneutral(tsol::TransientSolution,electrolyte)

Check for electroneutrality of transient solution

isincompressible(cx::Vector,electrolyte)

Check if the concentration vector fulfills incompressibility constraint, including the fact that the solvent concentration is nonnegative

isincompressible(tsol::TransientSolution,electrolyte)

Check for incompressibility of transient solution

c0_barc(u,electrolyte)

Calculate solvent concentration $c_0$ and summary concentration $\bar c$ from vector of concentrations c using the incompressibility constraint (assuming $κ_0=0$):

\[ \sum_{i=0}^N c_i (v_i + κ_iv_0) =1\]

This gives

\[ c_0v_0=1-\sum_{i=1}^N c_i (v_i+ κ_iv_0)\]

\[c_0= 1/v_0 - \sum_{i=1}^N c_i(\frac{v_i}{v_0}+κ_i)\]

Then we can calculate

\[ \bar c= \sum_{i=0}^N c_i\]

AbstractCellData

Abstract type for cell data. Users may implement their own subtypes.

electrolytes(::AbstractCellData)

Return a vector of electrolytes, one electrolyte for each cell region.

Subtypes of AbstractCellData can implement their own method for this.

working_electrode(d::AbstractCellData)

Return boundary index of working electrode. Subtypes of AbstractCellData can implement their own method for this.

bulk_electrode(d::AbstractCellData)

Return boundary index of bulk electrode. Subtypes of AbstractCellData can implement their own method for this.

norm_weights(d::AbstractCellData)

Weights of components in norm calculation. Subtypes of AbstractCellData can implement their own method for this.

working_electrode_voltage(d::AbstractCellData)

Return voltage applied at working electrode. Subtypes of AbstractCellData can implement their own method for this.

working_electrode_voltage!(d::AbstractCellData, v)

Set voltage applied at working electrode. Subtypes of AbstractCellData can implement their own method for this.

pressure_index(d::AbstractCellData)

Index of pressure in species list. Subtypes of AbstractCellData can implement their own method for this.

voltage_index(d::AbstractCellData)

Index of voltage in species list. Subtypes of AbstractCellData can implement their own method for this.

check_celldata(d::AbstractCellData)

Check if subtypes of AbstractCellData implement all necessary methods.

DGML_gamma!(γ, c, p, electrolyte)

Activity coefficients according to Dreyer, Guhlke, Müller, Landstorfer.

\[γ_i = \exp\left(\frac{\tilde v_i p}{RT}\right) \left(\frac{\bar c}{c_0}\right)^{m_i} \frac{1}{v_0\bar c}\]

Input:

• c: vector of concentrations
• p: pressure
• electrolyte: instance of ElectrolyteData

Output:

• γ (mutated) activity coefficients

abstract type AbstractElectrochemicalSystem

Abstract super type for electrochemical systems

unknowns(sys)

Return vector of unknowns initialized with bulk data.

PNPSystem(grid;
         celldata=ElectrolyteData(),
         bcondition=(f, u, n, e)-> nothing,
         reaction=(f, u, n, e)-> nothing,
         kwargs...)

Create VoronoiFVM system for generalized Poisson-Nernst-Planck. Input:

• grid: discretization grid
• celldata: instance of ElectrolyteData or of subtype of AbstractCellData
• bcondition: boundary condition
• reaction : reactions of the bulk species
• kwargs: Keyword arguments of VoronoiFVM.System

celldata(sys)

Extract celldata from system.

electrolytedata(sys)

Extract electrolyte data from system.

   solventconcentration(U::Array, electrolyte)

Calculate vector of solvent concentrations from solution array.

chemical_potentials(solution, electrolyte)

Calculate chemical potentials from solution. Returns μ0, μ, where μ0 is the vector of solvent chemical potentials, and μ is the nc×nnodes matrix of solute chemical potentials.

electrochemical_potentials(solution, electrolyte)

Calculate electrochemical potentials from solution. and return nc×nnodes matrix of solute electrochemical potentials measured in Volts.

μex_flux!(f,dϕ,ck,cl,γk,γl,electrolyte; evelo=0.0)

Excess chemical potential (Sedan) upweind flux, based on modified Scharfetter-Gummel scheme, see Gaudeul/Fuhrmann 2022.

Appearantly first described by Yu, Zhiping and Dutton, Robert, SEDAN III, www-tcad.stanford.edu/tcad/programs/sedan3.html

see also the 198? Fortran code available via https://web.archive.org/web/20210518233152/http://www-tcad.stanford.edu/tcad/programs/oldftpable.html

Verification calculation is in the paper.

act_flux!(ic,dϕ,ck,cl,γk,γl,electrolyte; evelo=0)

Scharfetter-Gummel upwind flux expression based on activities, see Fuhrmann, CPC 2015. As shown in Chainais-Hilliaret, Cances, Fuhrmann, Gaudeul, 2019, this is consistent to thermodynamic equilibrium but shows inferior convergence behavior.

cent_flux!(ic,dϕ,ck,cl,γk,γl,electrolyte; evelo = 0)

Flux expression based on central differences, see Gaudeul/Fuhrmann 2022. As shown in the paper this is covergent, consistent to thermodynamic equilibrium but may show inferior convergence behavior.

PBSystem(grid;
             celldata=ElectrolyteData(),
         bcondition=(f, u, n, e)-> nothing,
         kwargs...)

Create VoronoiFVM system generalized Poisson-Boltzmann. Input:

• grid: discretization grid
• celldata: composite struct containing electrolyte data
• bcondition: boundary condition
• kwargs: Keyword arguments of VoronoiFVM.System

struct PNPStokesSolver{Tflow, Tpnp}

Stokes solver struct.

• flowsolver::Any

• pnpsolver::Any

RLog(eps)

Callable struct for regularized logarithm. Smooth linear continuation for x<eps. This means we can calculate a "logarithm" of a small negative number. Objects mylog::RLog of this type are meant to replace the logarithm function and allow to invoke mylog(x).

RExp(trunc)

Callable struct for regularized exponential. Linear continuation for x>trunc, returns 1/rexp(-x) for x<-trunc. Objects myexp::RExp of this type are meant to replace the exponential function and allow to invoke myexp(x).