Double layer capcacitance comparison

begin
  using Pkg
  Pkg.activate(joinpath(@__DIR__, "..", "docs"))
  using Revise
  using CairoMakie
  using GridVisualize
  using Colors
end

begin
    using PlutoUI
    using VoronoiFVM
    using ExtendableGrids
    using LinearAlgebra
    using NLsolve
    using Unitful
    using LessUnitful
    using LessUnitful.MoreUnitful
    using LiquidElectrolytes
    using Test
end

begin
    mylog=RLog()
    LiquidElectrolytes.rlog(x::Number)=mylog(x)
end

1. Intro

This code implements the model described in Müller, R., Fuhrmann, J., & Landstorfer, M. (2020). Modeling polycrystalline electrode-electrolyte interfaces: The differential capacitance. Journal of The Electrochemical Society, 167(10), 106512

The code is part of the LiquidElectrolytes.jl package.

Equation numbers refer to the paper.

Concentrations are given in number densities, and calculations are done in mole fractions.

If not stated otherwise, all calculations and calculation results are in coherent SI units.

2. Model data

EquilibriumData

begin
    """
    	EquilibriumData
    Data structure containg data for equilibrum calculations
    """
    Base.@kwdef mutable struct EquilibriumData
        N::Int64 = 2                     # number of ionic species
        T::Float64 = 298.15 * ufac"K"        # temperature
        kT::Float64 = ph"k_B" * T             # temperature
        p_ref::Float64 = 1.0e5 * ufac"Pa"        # referece pressure
        pscale::Float64 = 1.0 * ufac"GPa"         # pressure scaling nparameter
        E_ref::Float64 = 0.0 * ufac"V"           # reference voltage
        n0_ref::Float64 = 55.508 * ph"N_A" / ufac"dm^3"  # solvent molarity
        v0::Float64 = 1 / n0_ref              # solvent molecule volume
        χ::Float64 = 15                    # dielectric susceptibility
        z::Vector{Int} = [-1, 1]                # ion charge numbers
        κ::Vector{Int} = [10, 10]               # ion solvation numbers
        molarity::Float64 = 0.1 * ph"N_A" / ufac"dm^3"
        n_E::Vector{Float64} = [molarity, molarity]  # bulk ion number densities
        μ_e::Vector{Float64} = [0.0]             # grain facet electron chemical potential

        e::Float64 = ph"e"
        ε_0::Float64 = ph"ε_0"

        v::Vector{Float64} = derived(κ, v0, n_E, T).v   # ion volumes
        y_E::Vector{Float64} = derived(κ, v0, n_E, T).y_E # bulk ion mole fractions
        y0_E::Float64 = derived(κ, v0, n_E, T).y0_E       # bulk solvent mole fraction
        U_T::Float64 = derived(κ, v0, n_E, T).U_T     # Temperature voltage k_BT/e0
    end

    function EquilibriumData(electrolyte::AbstractElectrolyteData)
        return EquilibriumData(;
            N = electrolyte.nc,
            T = electrolyte.T,
            p_ref = electrolyte.p_bulk,
            pscale = electrolyte.pscale,
            E_ref = electrolyte.ϕ_bulk,
            n0_ref = ph"N_A" / electrolyte.v0,
            χ = electrolyte.ε - 1.0,
            z = electrolyte.z,
            κ = electrolyte.κ,
            molarity = ph"N_A" * electrolyte.c_bulk[1],
            n_E = ph"N_A" * electrolyte.c_bulk
        )
    end
end

EquilibriumData

debyelength(data)

$$L_{Debye}=\sqrt{ \frac{(1+χ)ε_0k_BT}{e^2n_E}}$$

LiquidElectrolytes.debyelength(data::EquilibriumData) = sqrt((1 + data.χ) * data.ε_0 * data.kT / (ph"e"^2 * data.n_E[1]))

debyelength(EquilibriumData(molarity=0.01ph"N_A"/ufac"dm^3"))|>u"nm"

1.9422608647967765 nm

set_molarity!(data,M)

function set_molarity!(data::EquilibriumData, M_E)
    n_E = M_E * ph"N_A" / ufac"dm^3"
    data.molarity = n_E
    return data.n_E = fill(n_E, data.N)
end

set_molarity! (generic function with 1 method)

dlcap0(data)

Double layer capacitance at $φ=0$

$$C_{dl,0}=\sqrt{\frac{2(1+χ) ε_0e^2 n_E}{k_BT}}$$

LiquidElectrolytes.dlcap0(data::EquilibriumData) = sqrt(2 * (1 + data.χ) * ph"ε_0" * ph"e"^2 * data.n_E[1] / (ph"k_B" * data.T));

Check with Bard/Faulkner: the value must be 22.8 μF cm^-2

let
    data=EquilibriumData()
    set_molarity!(data,0.01)
    data.χ=78.49-1
    cdl0=dlcap0(data)|>u"μF/cm^2"
    @assert cdl0 ≈ 22.84669184882525u"μF/cm^2"
end

Species indices

begin
    const iφ = 1
    const ip = 2
    const iA = 1
    const iC = 2
end;

3. Model equations

Mole fractions

Equilibrium expression for mole fractions ($α≥0$) (16)

$$y_α(φ,p)=y_α^E\exp\left(\frac{-z_αe}{k_BT}(φ- φ^E)-\frac{v_α}{k_BT}(p-p^E)\right)$$

y_α(φ,p,α,data)

Ion molar fractions

function y_α(φ, p, α, data)
    η_φ = data.z[α] * data.e * (φ - data.E_ref)
    η_p = data.v[α] * (p * data.pscale - data.p_ref)
    return data.y_E[α] * exp(-(η_φ + η_p) / (data.kT))
end;

y0(p,data)

Solvent molar fraction

y0(p, data) = data.y0_E * exp(-data.v0 * (p * data.pscale - data.p_ref) / (data.kT));

Poisson equation

Poisson equation (32a)

$$-∇⋅(1+χ)ε_0∇φ = q(φ,p)$$

poisson_flux!(f,u,edge,data)

VoronoiFVM flux function for left hand side of Poisson equation

function poisson_flux!(f, u, edge, data)
    return f[iφ] = (1.0 + data.χ) * data.ε_0 * (u[iφ, 1] - u[iφ, 2])
end;

Space charge expression

Solvated ion volumes:

$$ v_α=(1+κ_α)v_0$$

Incompressibility condition (14)

$$\begin{aligned} 1&=∑\limits_α v_αn_α= n ∑\limits_α v_α y_α\\ n&=\frac1{∑\limits_α v_α y_α} \end{aligned}$$

Space charge

$$\begin{aligned} q(φ,p)&=e∑\limits_α z_αn_α = ne∑\limits_α z_αy_α\\ &=e\frac{∑\limits_α z_αy_α(\phi,p)}{∑\limits_α v_α y_α(\phi,p)} \end{aligned}$$

function spacecharge(φ, p, data)
    y = y0(p, data)
    sumyz = zero(eltype(p))
    sumyv = data.v0 * y
    for α in 1:(data.N)
        y = y_α(φ, p, α, data)
        sumyz += data.z[α] * y
        sumyv += data.v[α] * y
    end
    return data.e * sumyz / sumyv
end

spacecharge (generic function with 1 method)

Sum of mole fractions

function ysum(φ, p, data)
    sumy = y0(p, data)
    for α in 1:(data.N)
        sumy += y_α(φ, p, α, data)
    end
    return sumy
end

ysum (generic function with 1 method)

spacecharge_and_ysum!(f,u,node,data)

VoronoiFVM reaction function. This assumes that terms are on the left hand side. In addition to the space charge it calculates the residuum of the definition of $y_\alpha$ (32b):

$$∑_α y_α(φ,p)=1$$

However, direct usage of this equation leads to slow convergence of Newton's method. So we use

$$\log\left(∑_α y_α(φ,p)\right)=0$$

instead.

function spacecharge_and_ysum!(f, u, node, data)
    φ = u[iφ]
    p = u[ip]
    f[iφ] = -spacecharge(φ, p, data)
    return f[ip] = log(ysum(φ, p, data)) # this behaves much better with Newton's method
end;

update_derived!(data)

Update derived data in data record.

Calculate bulk mole fractions from incompressibiltiy:

$$\begin{aligned} ∑\limits_αv_αn_α^E&=1\\ n_0^E&=\frac1{v_0}\left(1-∑\limits_{α>0}v_αn_α^E\right)\\ n^E&=\frac1{v_0}\left(1-∑\limits_{α>0}v_αn_α^E\right)+ ∑\limits_{α>0}n_α^E\\ &=\frac1{v_0}\left(1-∑\limits_{α>0}(v_α-v_0)n_α^E\right)\\ &=\frac1{v_0}\left(1-∑\limits_{α>0}((1+ κ_α)v_0-v_0)n_α^E\right)\\ &=\frac1{v_0}\left(1-∑\limits_{α>0}κ_αv_0n_α^E\right)\\ &=\frac1{v_0}-∑\limits_{α>0}κ_αn_α^E\\ y_α^E&=\frac{n_α^E}{n^E} \end{aligned}$$

function derived(κ, v0, n_E, T)
    c0 = 1 / v0
    barc = 0.0
    v = (1.0 .+ κ) .* v0
    N = length(κ)
    for α in 1:N
        barc += n_E[α]
        c0 -= n_E[α] * (1 + κ[α])
    end
    barc += c0
    y_E = n_E / barc
    y0_E = c0 / barc
    U_T = ph"k_B" * T / ph"e"
    return (; v, y_E, y0_E, U_T)
end;

begin
    function update_derived!(data::EquilibriumData)
        (; κ, v0, n_E, T) = data
        return data.v, data.y_E, data.y0_E, data.U_T = derived(κ, v0, n_E, T)
    end
    update_derived!(::ElectrolyteData) = nothing
end

update_derived! (generic function with 2 methods)

let
    data = EquilibriumData()
    set_molarity!(data, 0.01)
    update_derived!(data)
    sumyz = 0.0
    sumyv = data.y0_E * data.v0
    sumy = data.y0_E
    for α in 1:data.N
        v = (1.0 + data.κ[α]) * data.v0
        sumyz += data.y_E[α] * data.z[α]
        sumyv += data.y_E[α] * v
        sumy += data.y_E[α]
    end
    @assert sumy ≈ 1.0
end

Bulk boundary condition

From (32d):

$$∇ φ\to 0$$

we choose homogeneous Neumann boundary conditions

$$\partial_n φ=0$$

which do not need any implementation.

Electrode boundary condition

See (32c) !!! bug in paper

$$φ|_{Σ_i}=\frac{1}{e} \mu_{e,i}- (E-E^{ref})$$

φ_Σ(ifacet,data,E)

Calculate potential boundary value for each facet from applied voltage E.

φ_Σ(ifacet, data, E) = data.μ_e[ifacet] / data.e - (E - data.E_ref);

4. System setup and solution

create_equilibrium_system(grid,data)

Create equlibrium system, enable species and apply zero voltage

function create_equilibrium_system(
        grid,
        data::EquilibriumData = EquilibriumData();
        Γ_bulk = 0
    )
    update_derived!(data)
    sys = VoronoiFVM.System(
        grid;
        data = data,
        flux = poisson_flux!,
        reaction = spacecharge_and_ysum!,
        species = [iφ, ip]
    )
    if Γ_bulk > 0
        boundary_dirichlet!(sys, iφ, Γ_bulk, 0.0)
    end
    return apply_voltage!(sys, 0)
end;

apply_voltage!(sys,E)

Apply voltage E to system.

function apply_voltage!(sys, E)
    data = sys.physics.data
    nbc = num_bfaceregions(sys.grid)
    nfacets = length(data.μ_e)
    @assert nbc > nfacets
    for ifacet in 1:nfacets
        boundary_dirichlet!(sys, iφ, ifacet, φ_Σ(ifacet, data, E))
    end
    return sys
end;

5. Postprocessing

calc_φ(sol,sys)

Obtain electrostatic potential from solution

calc_φ(sol, sys) = sol[iφ, :];

calc_p(sol,sys)

Obtain pressure from solution

calc_p(sol, sys) = sol[ip, :] * sys.physics.data.pscale;

c_num!(c,φ,p, data)

Calculate number concentration at discretization node

$$ n_α=ny_α$$

function c_num!(c, φ, p, data)
    y = y0(p, data)
    sumyv = data.v0 * y
    for α in 1:(data.N)
        c[α] = y_α(φ, p, α, data)
        sumyv += c[α] * data.v[α]
    end
    return c ./= sumyv
end;

function c0_num!(c, φ, p, data)
    y = y0(p, data)
    sumyv = data.v0 * y
    for α in 1:(data.N)
        c[α] = y_α(φ, p, α, data)
        sumyv += c[α] * data.v[α]
    end
    return y / sumyv
end;

calc_cnum(sol,sys)

Obtain ion number densities from system

function calc_cnum(sol, sys)
    data = sys.physics.data
    grid = sys.grid
    nnodes = num_nodes(grid)
    conc = zeros(data.N, nnodes)
    for i in 1:nnodes
        @views c_num!(conc[:, i], sol[iφ, i], sol[ip, i], data)
    end
    return conc
end;

function calc_c0num(sol, sys)
    data = sys.physics.data
    grid = sys.grid
    nnodes = num_nodes(grid)
    c0 = zeros(nnodes)
    conc = zeros(data.N)
    for i in 1:nnodes
        @views c0[i] = c0_num!(conc, sol[iφ, i], sol[ip, i], data)
    end
    return c0
end;

calc_cmol(sol,sys)

Obtain ion molarities (molar densities in mol/L) from system

calc_cmol(sol, sys) = calc_cnum(sol, sys) / (ph"N_A" * ufac"mol/dm^3");

calc_c0mol(sol, sys) = calc_c0num(sol, sys) / (ph"N_A" * ufac"mol/dm^3");

calc_QBL(sol,sys)

Obtain boundary layer charge (35)

$$Q^{BL}(E)=-\frac{1}{Σ} ∫_{Ω_E} q dx$$

calc_QBL(sol, sys) = VoronoiFVM.integrate(sys, spacecharge_and_ysum!, sol)[iφ, 1]

calc_QBL (generic function with 1 method)

Poly surface charge (34)

$$Q_s(E)= ∑_i s_i \hat Q_s(E-E^{ref} + \frac1{e}\mu_{e,i})$$

Single surface charge (26)

$$\hat Q_s = \frac{∑_\alpha z_αey_{s,α} +\sum_α \sum_β ν_{αβ}z_αey_{s,β}}{a_V^{ref}+ ∑_{α} a_α^{ref}z_αey_{s,α}} $$

We assume that there are no surface reactions, so we assume $\hat Q_s=0$ and $Q_s=0$.

6. Solution scenarios

dlcapsweep_equi(sys)

Calculate double layer capacitance. Return vector of voltages V and vector of double layer capacitances C.

function dlcapsweep_equi(
        sys; vmax = 2 * ufac"V", nsteps = 21, δV = 1.0e-3 * ufac"V",
        molarity = nothing,
        verbose = false
    )

    if !isnothing(molarity)
        error("The molarity kwarg of dlcapsweep_equie has been removed. Pass the molarity information with set_molarity!.")
    end


    data = sys.physics.data
    update_derived!(data)
    apply_voltage!(sys, 0)

    c = VoronoiFVM.NewtonControl()
    #	c.damp_growth=1.1
    c.verbose = verbose
    c.tol_round = 1.0e-10
    c.max_round = 3
    c.damp_initial = 0.01
    c.damp_growth = 2

    inival = solve(sys; inival = 0, control = c)
    vstep = vmax / (nsteps - 1)

    c.damp_initial = 1

    function rundlcap(dir)
        volts = zeros(0)
        caps = zeros(0)
        volt = 0.0
        sol = inival
        for iv in 1:nsteps
            apply_voltage!(sys, volt)
            c.damp_initial = 1
            sol = solve(sys; inival = sol, control = c)

            Q = calc_QBL(sol, sys)
            apply_voltage!(sys, volt + dir * δV)
            c.damp_initial = 1
            sol = solve(sys; inival = sol, control = c)
            Qδ = calc_QBL(sol, sys)
            push!(caps, (Q - Qδ) / (dir * δV))
            push!(volts, volt)
            volt += dir * vstep
        end
        return volts, caps
    end
    Vf, Cf = rundlcap(1)
    Vr, Cr = rundlcap(-1)
    return vcat(reverse(Vr), Vf), vcat(reverse(Cr), Cf)
end

dlcapsweep_equi (generic function with 1 method)

7. The pressure Poisson equation

An alternative possibility to handle the pressure has been introduced in

J. Fuhrmann, “Comparison and numerical treatment of generalised Nernst–Planck models,” Computer Physics Communications, vol. 196, pp. 166–178, 2015..

Starting with the momentum balance in mechanical equilibrium

$$ \nabla p = -q\nabla \varphi$$

by taking the divergence on both sides of the equation, one derives the pressure Poisson problem

$$\begin{aligned} -\Delta p &= \nabla\cdot q\nabla \varphi & \text{in}\; \Omega\\ p&=p_{bulk} & \text{on}\; \Gamma_{bulk}\\ (\nabla p + q\nabla \varphi)\cdot \vec n &=0 & \text{on}\; \partial\Omega\setminus\Gamma_{bulk}\\ \end{aligned}$$

The bulk Dirichlet boundary condition for the pressure is necessary to make the solution unique. It is reasonable to set the $\varphi$ to a bulk value at $\Gamma_{bulk}$ as well, and to calculate $p_{bulk}$ from the molar fraction sum constraint.

function spacecharge!(f, u, node, data)
    φ = u[iφ]
    p = u[ip]
    return f[iφ] = -spacecharge(u[iφ], u[ip], data)
end;

function poisson_and_p_flux!(f, u, edge, data)
    f[iφ] = (1.0 + data.χ) * data.ε_0 * (u[iφ, 1] - u[iφ, 2])
    q1 = spacecharge(u[iφ, 1], u[ip, 1], data)
    q2 = spacecharge(u[iφ, 2], u[ip, 2], data)
    return f[ip] = (u[ip, 1] - u[ip, 2]) + (u[iφ, 1] - u[iφ, 2]) * (q1 + q2) / (2 * data.pscale)
end;

function create_equilibrium_pp_system(
        grid,
        data::EquilibriumData = EquilibriumData();
        Γ_bulk = 0
    )
    update_derived!(data)

    sys = VoronoiFVM.System(
        grid; data = data, flux = poisson_and_p_flux!,
        reaction = spacecharge!, species = [iφ, ip]
    )
    if Γ_bulk > 0
        logysum!(y, p) = y[1] = log(ysum(0.0, p[1], data))
        res = nlsolve(logysum!, [0.0]; autodiff = :forward, method = :newton, xtol = 1.0e-10, ftol = 1.0e-20)
        boundary_dirichlet!(sys, iφ, Γ_bulk, 0.0)
        boundary_dirichlet!(sys, ip, Γ_bulk, res.zero[1])
    end
    return apply_voltage!(sys, 0)
end;

function ysum(sys, sol)
    data = sys.physics.data
    n = size(sol, 2)
    sumy = zeros(n)
    for i in 1:n
        sumy[i] = ysum(sol[iφ, i], sol[ip, i], data)
    end
    return sumy
end

ysum (generic function with 2 methods)

8. Tests

begin
    SI(x) = Float64(Unitful.ustrip(Unitful.upreferred(1 * x)))
    const V = SI(Unitful.V)
    const eV = SI(Unitful.eV)
    const nm = SI(Unitful.nm)
    const cm = SI(Unitful.cm)
    const μF = SI(Unitful.μF)
end

1.0e-6

elydata = ElectrolyteData(
    c_bulk = fill(0.01 * ufac"mol / dm^3", 2),
    κ = zeros(2)
)

ElectrolyteData(2, 0, 3, 4, [2.0e-9, 2.0e-9], [1, -1], 0.0180153, [0.0180153, 0.0180153], 1.805054151624549e-5, [1.805054151624549e-5, 1.805054151624549e-5], [0.0, 0.0], [10.0, 10.0], 0.0, 0.0, 2, 0.0, 1, 298.15, 2478.957029602388, 96485.33212331001, 78.49, 8.8541878128e-12, 1.0e9, false, :μex, true, [1.805054151624549e-5, 1.805054151624549e-5, 1.0, 0.0], 0.0, :DGL)

dlcap_exact = 0.22846691848825248

0.22846691848825249

@test dlcap0(elydata) ≈ dlcap_exact

Test Passed

begin
    equidata = EquilibriumData()
    set_molarity!(equidata, 0.01)
    equidata.χ = 78.49 - 1
end

77.49

@test dlcap0(equidata) |> unitfactor ≈ dlcap_exact

Test Passed

@test dlcap0(EquilibriumData(elydata)) ≈ dlcap_exact

Test Passed

begin
    Vmax = 2 * V

    L = 20nm

    hmin = 0.05 * nm

    hmax = 0.5 * nm

    X = ExtendableGrids.geomspace(0, L, hmin, hmax)

    grid = ExtendableGrids.simplexgrid(X)

    data = EquilibriumData(elydata)
end

EquilibriumData(2, 298.15, 4.1164049935e-21, 0.0, 1.0e9, 0.0, 3.3362659810399997e28, 2.9973629371369076e-29, 77.49, [1, -1], [0, 0], 6.02214076e24, [6.02214076e24, 6.02214076e24], [0.0], 1.602176634e-19, 8.8541878128e-12, [2.9973629371369076e-29, 2.9973629371369076e-29], [0.0001805054151624549, 0.0001805054151624549], 0.9996389891696752, 0.02569257912108585)

sys_sy = create_equilibrium_system(grid, data)

VoronoiFVM.System{Float64, Float64, Int32, Int64, Matrix{Int32}}(
  grid = ExtendableGrids.ExtendableGrid{Float64, Int32}(dim=1, nnodes=103, ncells=102,
  nbfaces=2),  
  physics = Physics(data=Main.var"workspace#3".EquilibriumData, flux=poisson_flux!,
  storage=default_storage, reaction=spacecharge_and_ysum!, ),  
  num_species = 2)

sys_pp = create_equilibrium_pp_system(grid, data, Γ_bulk = 2)

VoronoiFVM.System{Float64, Float64, Int32, Int64, Matrix{Int32}}(
  grid = ExtendableGrids.ExtendableGrid{Float64, Int32}(dim=1, nnodes=103, ncells=102,
  nbfaces=2),  
  physics = Physics(data=Main.var"workspace#3".EquilibriumData,
  flux=poisson_and_p_flux!, storage=default_storage, reaction=spacecharge!, ),  
  num_species = 2)

inival = unknowns(sys_sy, inival = 0);

molarities = [0.001, 0.01, 0.1]

3-element Vector{Float64}:
 0.001
 0.01
 0.1

function capscalc(sys, molarities)
    result = []
    for imol in 1:length(molarities)
        if isa(sys.physics.data, EquilibriumData)
            set_molarity!(sys.physics.data, molarities[imol])
            t = @elapsed volts, caps = dlcapsweep_equi(sys, vmax = 1V, nsteps = 101)
        else
            sys.physics.data.c_bulk .= molarities[imol] * ufac"mol/dm^3"
            t = @elapsed r = dlcapsweep(
                sys,
                voltages = range(-1, 1, length = 201)
            )
            volts = voltages(r)
            caps = r.dlcaps
        end
        cdl0 = dlcap0(sys.physics.data)
        @info "elapsed=$(t)"
        push!(
            result,
            (
                voltages = volts,
                dlcaps = caps,
                cdl0 = cdl0,
                molarity = molarities[imol],
            )
        )
    end
    return result
end

capscalc (generic function with 1 method)

Algebraic pressure equation

result_sy = capscalc(sys_sy, molarities)

3-element Vector{Any}:
 (voltages = [-1.0000000000000007, -0.9900000000000007, -0.9800000000000006, -0.9700000000000006, -0.9600000000000006, -0.9500000000000006, -0.9400000000000006, -0.9300000000000006, -0.9200000000000006, -0.9100000000000006  …  0.9100000000000006, 0.9200000000000006, 0.9300000000000006, 0.9400000000000006, 0.9500000000000006, 0.9600000000000006, 0.9700000000000006, 0.9800000000000006, 0.9900000000000007, 1.0000000000000007], dlcaps = [1.6054456294094521, 1.6167384247496308, 1.6282615466951533, 1.640020213455884, 1.6520214461541727, 1.664274423724521, 1.6767906669992172, 1.6895840173747345, 1.7026703965759005, 1.7160673594509568  …  1.7160673594500686, 1.7026703965754564, 1.6895840173742904, 1.6767906669992172, 1.664274423724077, 1.6520214461546168, 1.640020213455884, 1.6282615466960415, 1.6167384247496308, 1.6054456294098962], cdl0 = 0.07224758324229108, molarity = 0.001)
 (voltages = [-1.0000000000000007, -0.9900000000000007, -0.9800000000000006, -0.9700000000000006, -0.9600000000000006, -0.9500000000000006, -0.9400000000000006, -0.9300000000000006, -0.9200000000000006, -0.9100000000000006  …  0.9100000000000006, 0.9200000000000006, 0.9300000000000006, 0.9400000000000006, 0.9500000000000006, 0.9600000000000006, 0.9700000000000006, 0.9800000000000006, 0.9900000000000007, 1.0000000000000007], dlcaps = [1.5431651974693317, 1.5531651801965474, 1.5633745223961881, 1.5737978182785994, 1.5844387601986476, 1.5953004223328016, 1.6063856678210087, 1.6176976474433502, 1.629240344612537, 1.6410191116280615  …  1.6410191116280615, 1.629240344612537, 1.6176976474433502, 1.6063856678205646, 1.5953004223332456, 1.5844387601990917, 1.5737978182785994, 1.5633745223961881, 1.5531651801969915, 1.5431651974688876], cdl0 = 0.22846691848825246, molarity = 0.01)
 (voltages = [-1.0000000000000007, -0.9900000000000007, -0.9800000000000006, -0.9700000000000006, -0.9600000000000006, -0.9500000000000006, -0.9400000000000006, -0.9300000000000006, -0.9200000000000006, -0.9100000000000006  …  0.9100000000000006, 0.9200000000000006, 0.9300000000000006, 0.9400000000000006, 0.9500000000000006, 0.9600000000000006, 0.9700000000000006, 0.9800000000000006, 0.9900000000000007, 1.0000000000000007], dlcaps = [1.487924343455127, 1.4968306876359527, 1.50590131135786, 1.5151449781125592, 1.5245703702571056, 1.534185754234585, 1.5439986856011778, 1.554015788434171, 1.564242639701785, 1.5746837819499149  …  1.5746837819494708, 1.564242639701785, 1.554015788434171, 1.543998685602066, 1.534185754235029, 1.5245703702579938, 1.5151449781125592, 1.50590131135786, 1.4968306876355086, 1.4879243434542389], cdl0 = 0.7224758324229108, molarity = 0.1)

Pressure poisson problem

result_pp = capscalc(sys_pp, molarities)

3-element Vector{Any}:
 (voltages = [-1.0000000000000007, -0.9900000000000007, -0.9800000000000006, -0.9700000000000006, -0.9600000000000006, -0.9500000000000006, -0.9400000000000006, -0.9300000000000006, -0.9200000000000006, -0.9100000000000006  …  0.9100000000000006, 0.9200000000000006, 0.9300000000000006, 0.9400000000000006, 0.9500000000000006, 0.9600000000000006, 0.9700000000000006, 0.9800000000000006, 0.9900000000000007, 1.0000000000000007], dlcaps = [1.6053696728359768, 1.616661957323906, 1.6281845565306874, 1.6399426876914092, 1.6519433710473486, 1.6641957847953748, 1.6767114492113322, 1.689504205306669, 1.7025899746241357, 1.7159863120124186  …  1.7159863120110863, 1.7025899746236917, 1.689504205306669, 1.6767114492113322, 1.6641957847953748, 1.6519433710464604, 1.6399426876914092, 1.6281845565306874, 1.6166619573243501, 1.6053696728359768], cdl0 = 0.07224758324229108, molarity = 0.001)
 (voltages = [-1.0000000000000007, -0.9900000000000007, -0.9800000000000006, -0.9700000000000006, -0.9600000000000006, -0.9500000000000006, -0.9400000000000006, -0.9300000000000006, -0.9200000000000006, -0.9100000000000006  …  0.9100000000000006, 0.9200000000000006, 0.9300000000000006, 0.9400000000000006, 0.9500000000000006, 0.9600000000000006, 0.9700000000000006, 0.9800000000000006, 0.9900000000000007, 1.0000000000000007], dlcaps = [1.5431581779439263, 1.553158121668119, 1.5633674242052642, 1.5737906797639312, 1.584431580668344, 1.5952932010696586, 1.6063784040536433, 1.6176903403488652, 1.6292329932814376, 1.641011715070917  …  1.6410117150704728, 1.6292329932809935, 1.6176903403493093, 1.6063784040540874, 1.5952932010692145, 1.5844315806679, 1.5737906797639312, 1.56336742420482, 1.553158121668119, 1.5431581779443704], cdl0 = 0.22846691848825246, molarity = 0.01)
 (voltages = [-1.0000000000000007, -0.9900000000000007, -0.9800000000000006, -0.9700000000000006, -0.9600000000000006, -0.9500000000000006, -0.9400000000000006, -0.9300000000000006, -0.9200000000000006, -0.9100000000000006  …  0.9100000000000006, 0.9200000000000006, 0.9300000000000006, 0.9400000000000006, 0.9500000000000006, 0.9600000000000006, 0.9700000000000006, 0.9800000000000006, 0.9900000000000007, 1.0000000000000007], dlcaps = [1.4879243433880696, 1.4968306875720039, 1.505901311292579, 1.51514497804639, 1.524570370194489, 1.5341857541719683, 1.5439986855363408, 1.554015788369334, 1.564242639637392, 1.5746837818828574  …  1.5746837818828574, 1.564242639636948, 1.554015788369334, 1.5439986855354526, 1.5341857541719683, 1.5245703701936009, 1.51514497804639, 1.505901311292135, 1.4968306875720039, 1.4879243433885136], cdl0 = 0.7224758324229108, molarity = 0.1)

Poisson Nernst-Planck from LiquidElectrolytes

function pnp_bcondition(f, u, bnode, data::ElectrolyteData)
    (; iϕ, Γ_we, ϕ_we) = data
    boundary_dirichlet!(f, u, bnode, species = iϕ, region = Γ_we, value = ϕ_we)
    return bulkbcondition(f, u, bnode, data)
end

pnp_bcondition (generic function with 1 method)

sys_pnp = PNPSystem(grid; bcondition = pnp_bcondition, celldata = elydata)

VoronoiFVM.System{Float64, Float64, Int32, Int64, Matrix{Int32}}(
  grid = ExtendableGrids.ExtendableGrid{Float64, Int32}(dim=1, nnodes=103, ncells=102,
  nbfaces=2),  
  physics = Physics(data=LiquidElectrolytes.ElectrolyteData, flux=pnpflux,
  storage=pnpstorage, reaction=_pnpreaction, breaction=pnp_bcondition, ),  
  num_species = 4)

result_pnp = capscalc(sys_pnp, molarities)

3-element Vector{Any}:
 (voltages = [-1.0, -0.99, -0.98, -0.97, -0.96, -0.95, -0.94, -0.93, -0.92, -0.91  …  0.91, 0.92, 0.93, 0.94, 0.95, 0.96, 0.97, 0.98, 0.99, 1.0], dlcaps = [1.6059847397498217, 1.6172895790589337, 1.6288249802798305, 1.6405962480314784, 1.6526105132363966, 1.664877079909921, 1.6774075984216097, 1.6902160325571458, 1.703318407355603, 1.7167323510358656  …  1.7165966087828721, 1.703185871697599, 1.6900865214530825, 1.6772809426024438, 1.6647531290603368, 1.6524891391211938, 1.6404773460942934, 1.6287084689103892, 1.6171753969107527, 1.605872841508038], cdl0 = 0.0722475832422911, molarity = 0.001)
 (voltages = [-1.0, -0.99, -0.98, -0.97, -0.96, -0.95, -0.94, -0.93, -0.92, -0.91  …  0.91, 0.92, 0.93, 0.94, 0.95, 0.96, 0.97, 0.98, 0.99, 1.0], dlcaps = [1.5437027401388193, 1.553714059565614, 1.5639350174501132, 1.5743701532011656, 1.585023116450479, 1.5958969580598392, 1.6069945430263743, 1.6183190533602954, 1.6298745338794163, 1.6416664263640968  …  1.6415473149899995, 1.6297578193480433, 1.6182046727397648, 1.606882449771696, 1.5957871169547388, 1.5849154982872093, 1.574264729065078, 1.5638317544208746, 1.5536129170889978, 1.5436036676863552], cdl0 = 0.22846691848825249, molarity = 0.01)
 (voltages = [-1.0, -0.99, -0.98, -0.97, -0.96, -0.95, -0.94, -0.93, -0.92, -0.91  …  0.91, 0.92, 0.93, 0.94, 0.95, 0.96, 0.97, 0.98, 0.99, 1.0], dlcaps = [1.488410003780416, 1.4973251551175792, 1.5064050640933857, 1.5156584988273991, 1.5250941276301333, 1.5347201849191805, 1.5445441801809068, 1.554572681734534, 1.5648112081212417, 1.5752642480171275  …  1.57515864327884, 1.564707767669482, 1.5544713653525832, 1.5444449380064285, 1.5346229564272562, 1.5249988419174798, 1.5155650769083095, 1.5063134206405593, 1.4972352024900104, 1.4883216552963319], cdl0 = 0.7224758324229109, molarity = 0.1)

Poisson-Boltzmann from LiquidElectrolytes

function pb_bcondition(f, u, bnode, data)
    (; Γ_we, Γ_bulk, ϕ_we, iϕ, ip) = data
    ## Dirichlet ϕ=ϕ_we at Γ_we
    boundary_dirichlet!(f, u, bnode, species = iϕ, region = Γ_we, value = ϕ_we)
    boundary_dirichlet!(f, u, bnode, species = iϕ, region = Γ_bulk, value = data.ϕ_bulk)
    return boundary_dirichlet!(f, u, bnode, species = ip, region = Γ_bulk, value = data.p_bulk)

end

pb_bcondition (generic function with 1 method)

sys_pb = PBSystem(grid; celldata = deepcopy(elydata), bcondition = pb_bcondition)

VoronoiFVM.System{Float64, Float64, Int32, Int64, Matrix{Int32}}(
  grid = ExtendableGrids.ExtendableGrid{Float64, Int32}(dim=1, nnodes=103, ncells=102,
  nbfaces=2),  
  physics = Physics(data=LiquidElectrolytes.ElectrolyteData, flux=pbflux,
  storage=default_storage, reaction=pbreaction, breaction=pb_bcondition, ),  
  num_species = 4)

result_pb = capscalc(sys_pb, molarities)

3-element Vector{Any}:
 (voltages = [-1.0, -0.99, -0.98, -0.97, -0.96, -0.95, -0.94, -0.93, -0.92, -0.91  …  0.91, 0.92, 0.93, 0.94, 0.95, 0.96, 0.97, 0.98, 0.99, 1.0], dlcaps = [1.6059847397675853, 1.617289579085579, 1.6288249802798305, 1.6405962480359193, 1.6526105132408375, 1.6648770798965984, 1.6774075984216097, 1.6902160325615867, 1.7033184073644847, 1.7167323510358656  …  1.7165966087784312, 1.7031858717020398, 1.6900865214530825, 1.6772809426024438, 1.6647531290470141, 1.6524891391256347, 1.6404773461031752, 1.6287084689015074, 1.61717539689743, 1.6058728415213608], cdl0 = 0.0722475832422911, molarity = 0.001)
 (voltages = [-1.0, -0.99, -0.98, -0.97, -0.96, -0.95, -0.94, -0.93, -0.92, -0.91  …  0.91, 0.92, 0.93, 0.94, 0.95, 0.96, 0.97, 0.98, 0.99, 1.0], dlcaps = [1.5437027401965508, 1.5537140596144638, 1.5639350174456723, 1.5743701531567567, 1.5850231164460382, 1.5958969580509574, 1.6069945430308152, 1.6183190533691771, 1.6298745338749754, 1.6416664263640968  …  1.6415473149855586, 1.6297578193391615, 1.6182046727397648, 1.606882449771696, 1.5957871169547388, 1.5849154982650049, 1.5742647290473144, 1.563831754407552, 1.5536129170712343, 1.5436036676597098], cdl0 = 0.22846691848825249, molarity = 0.01)
 (voltages = [-1.0, -0.99, -0.98, -0.97, -0.96, -0.95, -0.94, -0.93, -0.92, -0.91  …  0.91, 0.92, 0.93, 0.94, 0.95, 0.96, 0.97, 0.98, 0.99, 1.0], dlcaps = [1.4884100039180836, 1.4973251551486655, 1.5064050640933857, 1.5156584988496036, 1.5250941276301333, 1.534720184941385, 1.544544180167584, 1.554572681730093, 1.5648112081212417, 1.5752642480082457  …  1.5751586432743991, 1.564707767669482, 1.5544713653437015, 1.5444449380019876, 1.5346229564006109, 1.524998841908598, 1.5155650769171913, 1.5063134206805273, 1.4972352025965918, 1.488321655291891], cdl0 = 0.7224758324229109, molarity = 0.1)

Result comparison

function resultcompare(r1, r2; tol = 1.0e-3)
    for i in 1:length(r1)
        for f in fieldnames(typeof(r1[i]))
            if !isapprox(r1[i][f], r2[i][f]; rtol = tol)
                return false
            end
        end
    end
    return true
end

resultcompare (generic function with 1 method)

@test resultcompare(result_pp, result_sy; tol = 5.0e-3)

Test Passed

@test resultcompare(result_pb, result_pp; tol = 5.0e-3)

Test Passed

@test resultcompare(result_pp, result_pnp; tol = 5.0e-3)

Test Passed

@test resultcompare(result_pb, result_pnp; tol = 1.0e-10)

Test Passed

Result plot

Compare with Fig 4.2 of Fuhrmann (2015)

function capsplot(vis, result, title)
    hmol = 1 / length(result)
    for imol in 1:length(result)
        c = RGB(1 - imol * hmol, 0, imol * hmol)
        scalarplot!(
            vis, result[imol].voltages, result[imol].dlcaps / (μF / cm^2),
            color = c, clear = false, label = "$(result[imol].molarity)M", markershape = :none, title = title
        )
        scalarplot!(
            vis, [0], [result[imol].cdl0] / (μF / cm^2),
            clear = false, markershape = :circle, markersize = 8, label = ""
        )
    end
    return vis
end

capsplot (generic function with 1 method)

let
    vis = GridVisualizer(Plotter = CairoMakie, legend = :lt, layout = (2, 2), size = (650, 650))

    capsplot(vis[1, 1], result_sy, "Algebraic pressure")
    capsplot(vis[1, 2], result_pp, "Pressure Poisson")
    capsplot(vis[2, 1], result_pb, "Poisson-Boltzmann")
    capsplot(vis[2, 2], result_pnp, "Poisson-Nernst-Planck")

    reveal(vis)
end