VoronoiFVMPoisson.jl

Realization of Poisson-Boltzmann and Poisson-Bikerman equations based on plain VoronoiFVM.jl (without LiquidElectrolytes). This notebook is meant for demonstation purposes and may be moved to VoronoiFVM.jl.

begin


    using PlutoUI, HypertextLiteral, UUIDs
    using LinearAlgebra
    using Interpolations
    using VoronoiFVM, GridVisualize, ExtendableGrids
    using LaTeXStrings
    using LessUnitful, Unitful

    if isdefined(Main, :PlutoRunner)
        using CairoMakie
        default_plotter!(CairoMakie)
        CairoMakie.activate!(type = "svg")
    end

end

Constants, Units & Parameters

See LessUnitful.jl ph and ufac and Unitful.jl for u

begin
    const ε = 78.49
    const ε_0 = ph"ε_0"
    const F = ph"N_A" * ph"e"
    const K = ufac"K"
    const nm = ufac"nm"
    const dm = ufac"dm"
    const V = ufac"V"
    const mol = ufac"mol"
end;

Parameters:

begin
    const molarity = 0.01
    const c_bulk = [molarity, molarity] * mol / dm^3 # bulk ion molar concentrations
    const z = [-1, 1] # species charge numbers
    const c̄ = 55.508mol / dm^3 # summary molar concentration
    const nref = 0 # grid refinement level
    const T = 298.15K # temperature
    const L = 20.0nm # computational domain size
    const vmax = 1.0 * V # maximum voltage
end;

Check for bulk electroneutrality:

@assert dot(c_bulk, z) == 0

Derived data:

begin
    const RT = ph"R" * T
    const c0_bulk = c̄ - sum(c_bulk) # solvent bulk molar concentration
    const hmin = 1.0e-1 * nm * 2.0^(-nref) # grid size at working electrode
    const hmax = 1.0 * nm * 2.0^(-nref) # grid size at bulk
    const l_debye = sqrt(ε * ε_0 * RT / (F^2 * c_bulk[1])) # Debye length
    const dlcap0 = sqrt(2 * ε * ε_0 * F^2 * c_bulk[1] / RT) # Double layer capacitance at point of zero charge (0V)
end;

Debye length= 4.3018 nm

Double layer capacitance at zero voltage for symmetric binary electrolyte = 22.847 μF cm^-2

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

52-element Vector{Float64}:
 0.0
 1.0000000000000002e-10
 2.0469806721888107e-10
 3.1431492001257436e-10
 4.290816462337372e-10
 5.492401903976795e-10
 6.750438637356726e-10
 ⋮
 1.5461201801124036e-8
 1.6287579454587685e-8
 1.715278088569287e-8
 1.8058630061610094e-8
 1.9007036640713593e-8
 2.0e-8

grid = simplexgrid(X)

ExtendableGrids.ExtendableGrid{Float64, Int32}
      dim =       1
   nnodes =      52
   ncells =      51
  nbfaces =       2

gridplot(grid, size = (600, 200))

Nonlinear Poisson equation

Let $\bar c$ be the summary concentration, and $y_i$ be the solute molar fractions, such that for the species molar concentrations one has $c_i=\bar c y_i$. Then the Poisson equation states:

$$\begin{aligned} -\nabla \cdot εε_0 \nabla \phi &= F\sum_{i=1}^n z_i \bar c y_i\\ \end{aligned}$$

Un VoronoiFVM, the left hand side is expressed via:

function pb_flux(y, u, edge, data)
    return y[1] = ε * ε_0 * (u[1, 1] - u[1, 2])
end

pb_flux (generic function with 1 method)

It remains to state the dependecy of the mole fractions on the potential.

We choose Dirichlet boundary conditions at the electrode (index 1) and the bulk (index 2):

const bcvals = [0.0, 0.0]

2-element Vector{Float64}:
 0.0
 0.0

function pb_bc(y, u, bnode, data)
    boundary_dirichlet!(y, u, bnode, species = 1, region = 2, value = bcvals[2])
    return boundary_dirichlet!(y, u, bnode, species = 1, region = 1, value = bcvals[1])
end

pb_bc (generic function with 1 method)

Classical Poisson-Boltzmann

$$\begin{aligned} y_i&= \frac{c_{i,bulk}}{\bar c}\exp\left(z_i\phi \frac{F}{RT}\right) \end{aligned}$$

pbo_molfrac(ϕ, i) = (c_bulk[i] / c̄) * exp(z[i] * ϕ * F / RT)

pbo_molfrac (generic function with 1 method)

function pbo_spacecharge(ϕ)
    sumyz = zero(ϕ)
    for i in 1:length(z)
        sumyz += z[i] * pbo_molfrac(ϕ, i)
    end
    return F * c̄ * sumyz
end

pbo_spacecharge (generic function with 1 method)

function pbo_reaction(y, u, node, data)
    return y[1] = pbo_spacecharge(u[1])
end

pbo_reaction (generic function with 1 method)

pbo_system = VoronoiFVM.System(
    grid;
    reaction = pbo_reaction,
    flux = pb_flux,
    bcondition = pb_bc,
    species = [1]
)

VoronoiFVM.System{Float64, Float64, Int32, Int64, Matrix{Int32}}(
  grid = ExtendableGrids.ExtendableGrid{Float64, Int32}(dim=1, nnodes=52, ncells=51,
  nbfaces=2),  
  physics = Physics(flux=pb_flux, storage=default_storage, reaction=pbo_reaction,
  breaction=pb_bc, ),  
  num_species = 1)

function sweep(system; nsweep = 100, δ = 1.0e-4)
    voltages = range(0, vmax, length = nsweep + 1)
    solutions = []
    dlcaps = zeros(0)
    for v in voltages
        bcvals[1] = v
        if length(solutions) == 0
            inival = 0.0
        else
            inival = solutions[end]
        end
        push!(solutions, solve(system; inival))
        q = integrate(system, system.physics.reaction, solutions[end])
        bcvals[1] = v + δ
        solδ = solve(system; inival = solutions[end])
        qδ = integrate(system, system.physics.reaction, solδ)
        push!(dlcaps, (qδ[1] - q[1]) / δ)
    end
    return voltages, dlcaps, TransientSolution(solutions, voltages)
end

sweep (generic function with 1 method)

pbo_volts, pbo_caps, pbo_sols = sweep(pbo_system)

(0.0:0.01:1.0, [0.2279211284512289, 0.23232155094885917, 0.24560434836521683, 0.26828027521705844, 0.301224127664014, 0.34571356466221725, 0.40348767573306077, 0.47683067960557074, 0.5686896505781303, 0.6828408873541664  …  4.535790741564179e12, 6.694035309081421e12, 9.879227520047607e12, 1.458001517569336e13, 2.1517557125893555e13, 3.175615793874756e13, 4.686654535775879e13, 6.916683932619873e13, 1.0207818020028808e14, 1.506495739066797e14], TransientSolution{Any, 3, Vector{Any}, StepRangeLen{Float64, Base.TwicePrecision{Float64}, Base.TwicePrecision{Float64}, Int64}}(Any[[0.0 0.0 … 0.0 0.0], [0.01 0.009674563602456738 … 9.237370354343381e-6 6.465143843346336e-36], [0.02 0.01933737655569516 … 1.8302825687845453e-5 1.2809966069615975e-35], [0.03 0.028976332547113244 … 2.7037039748681027e-5 1.8922955816241113e-35], [0.04 0.03857861313360881 … 3.530395636764442e-5 2.470888872055713e-35], [0.05 0.04813033116640495 … 4.299813647931687e-5 3.0093969027009016e-35], [0.06 0.05761617795125789 … 5.0048127083344423e-5 3.502818748971762e-35], [0.07 0.06701908305595497 … 5.641602729517804e-5 3.948501765569903e-35], [0.08 0.07631990299437341 … 6.209404779787784e-5 4.3459007859324827e-35], [0.09 0.08549716482188524 … 6.709916609103969e-5 4.696204048408479e-35]  …  [0.91 0.2723334589126947 … 9.830385041770582e-5 6.880188938255974e-35], [0.92 0.2727500282846376 … 9.831524078776162e-5 6.880986138952869e-35], [0.93 0.2731603776598562 … 9.83264020040468e-5 6.881767301404764e-35], [0.94 0.27356468421654406 … 9.83373419242993e-5 6.882532975566934e-35], [0.95 0.27396311783185523 … 9.834806803279557e-5 6.883283685256463e-35], [0.96 0.2743558414708359 … 9.835858746279237e-5 6.884019929722911e-35], [0.97 0.27474301155016845 … 9.836890701734921e-5 6.884742185105666e-35], [0.98 0.27512477827863524 … 9.837903318867109e-5 6.885450905787713e-35], [0.99 0.2755012859760444 … 9.838897217608784e-5 6.886146525653989e-35], [1.0 0.2758726733722127 … 9.839872990278725e-5 6.886829459262535e-35]], 0.0:0.01:1.0, NewtonSolverHistory[]))

plotcdl(pbo_volts, pbo_caps, pb_v)

plotsols(pbo_sols, pb_v, pbo_molfrac)

"Poisson-Bikerman"

Here we use the expression

$$\begin{aligned} y_i&= \frac{\frac{c_{i,bulk}}{c_{0,bulk}}\exp\left(z_i\phi \frac{F}{RT}\right)}{1+\sum_{j=1}^n \frac{c_{j,bulk}}{c_{0,bulk}}\exp\left(z_i\phi \frac{F}{RT}\right)} \end{aligned}$$

function pbi_molfrac(ϕ, j)
    N = length(z)
    denom = one(ϕ)
    for i in 1:length(z)
        denom += exp(z[i] * ϕ * F / RT) * c_bulk[i] / c0_bulk
    end
    return exp(z[j] * ϕ * F / RT) * c_bulk[j] / (c0_bulk * denom)
end

pbi_molfrac (generic function with 1 method)

function pbi_spacecharge(ϕ)
    sumyz = zero(ϕ)
    for i in 1:length(z)
        sumyz += z[i] * pbi_molfrac(ϕ, i)
    end
    return F * c̄ * sumyz
end

pbi_spacecharge (generic function with 1 method)

function pbi_reaction(y, u, node, data)
    return y[1] = pbi_spacecharge(u[1])
end

pbi_reaction (generic function with 1 method)

pbi_system = VoronoiFVM.System(grid; reaction = pbi_reaction, flux = pb_flux, bcondition = pb_bc, species = [1])

VoronoiFVM.System{Float64, Float64, Int32, Int64, Matrix{Int32}}(
  grid = ExtendableGrids.ExtendableGrid{Float64, Int32}(dim=1, nnodes=52, ncells=51,
  nbfaces=2),  
  physics = Physics(flux=pb_flux, storage=default_storage, reaction=pbi_reaction,
  breaction=pb_bc, ),  
  num_species = 1)

pbi_volts, pbi_caps, pbi_sols = sweep(pbi_system)

(0.0:0.01:1.0, [0.227921128294752, 0.23231665625685072, 0.24558296238934468, 0.26822443355533837, 0.3011027537085965, 0.34547100866769154, 0.4030233712739098, 0.47596242949837775, 0.5670867874531943, 0.6799016647912604  …  1.6615077566939362, 1.6484924477255802, 1.6351761185573466, 1.621610192410472, 1.6078596289670344, 1.5940011669801635, 1.5801209693666607, 1.5663116858721793, 1.5526690268252707, 1.5392880156506195], TransientSolution{Any, 3, Vector{Any}, StepRangeLen{Float64, Base.TwicePrecision{Float64}, Base.TwicePrecision{Float64}, Int64}}(Any[[0.0 0.0 … 0.0 0.0], [0.01 0.009674565740806013 … 9.237401913675033e-6 6.46516593140938e-36], [0.02 0.01933739458710759 … 1.8303078174837718e-5 1.281014278275629e-35], [0.03 0.02897639891080578 … 2.7037892148765896e-5 1.892355240260166e-35], [0.04 0.03857879032931354 … 3.530597861058703e-5 2.4710304068324803e-35], [0.05 0.04813073300819722 … 4.3002093324967484e-5 3.009673838401386e-35], [0.06 0.05761700690912229 … 5.0054986396716e-5 3.5032988255077535e-35], [0.07 0.06702069411926119 … 5.642697455746507e-5 3.9492679535937835e-35], [0.08 0.07632291087241937 … 6.211050841066457e-5 4.3470528478866273e-35], [0.09 0.08550262057379068 … 6.712283647877935e-5 4.697860715356911e-35]  …  [0.91 0.6216348682552673 … 0.00010614269340449226 7.428821779973225e-35], [0.92 0.6292533108610063 … 0.00010624974683503027 7.436314343341397e-35], [0.93 0.6368907011963554 … 0.00010635660638928038 7.443793337500572e-35], [0.94 0.6445474374833199 … 0.00010646329044417423 7.451260048633965e-35], [0.95 0.6522238341336632 … 0.00010656980411580769 7.458714834812039e-35], [0.96 0.6599201034526659 … 0.00010667613931128036 7.466157129620734e-35], [0.97 0.6676363402950177 … 0.00010678227517388479 7.473585473319788e-35], [0.98 0.6753725105302449 … 0.00010688817895392738 7.48099757379065e-35], [0.99 0.683128444098655 … 0.00010699380731248986 7.488390397785434e-35], [1.0 0.6909038332509377 … 0.0001070991080341135 7.495760290796002e-35]], 0.0:0.01:1.0, NewtonSolverHistory[]))

plotcdl (generic function with 1 method)

function plotsols(pbsols, pbv, pb_molfrac)
    vis = GridVisualizer(size = (700, 200), layout = (1, 2), xlabel = "x/nm")
    ϕ = pbsols(pbv)[1, :]
    cp = c̄ * pb_molfrac.(ϕ, 1) / ufac"mol/dm^3"
    cm = c̄ * pb_molfrac.(ϕ, 2) / ufac"mol/dm^3"
    scalarplot!(vis[1, 1], X / nm, ϕ, color = :green, ylabel = "ϕ/V")
    scalarplot!(vis[1, 2], X / nm, cp, color = :red, ylabel = "c/(mol/L)")
    scalarplot!(vis[1, 2], X / nm, cm, color = :blue, clear = false)
    return reveal(vis)
end

plotsols (generic function with 1 method)