VoronoiFVMPoisson2.jl

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

It demonstrates a way to calculate the double layer charge as a unknown of the system.

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))
const nv = ones(num_nodes(grid))
52-element Vector{Float64}:
 1.0
 1.0
 1.0
 1.0
 1.0
 1.0
 1.0
 ⋮
 1.0
 1.0
 1.0
 1.0
 1.0
 1.0

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 - q= 0\\ q&=F\sum_{i=1}^n z_i \bar c y_i\\ Q&=\int_\Omega q d\omega \end{aligned}$$

In this version of the implementation, we introduce three species: Electrostatic potential \(\phi\), charge density \(q\) and double layer charge \(Q\).

begin
    const iφ = 1
    const iq = 2
    const iQ = 3
end
3

In VoronoiFVM, the left hand side is expressed via:

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

It remains to state the dependency 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 = iφ, region = 2, value = bcvals[2])
    boundary_dirichlet!(y, u, bnode, species = iφ, region = 1, value = bcvals[1])
    return
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_spacechargedensity(ϕ)
    sumyz = zero(ϕ)
    for i in 1:length(z)
        sumyz += z[i] * pbo_molfrac(ϕ, i)
    end
    return F * c̄ * sumyz
end
pbo_spacechargedensity (generic function with 1 method)
function pbo_reaction(y, u, node, data)
    y[iφ] = u[iq] # implements the  `-q` part of of the lhs of the poission equation
    y[iq] = u[iq] - pbo_spacechargedensity(u[iφ]) # implements q=F\sum z_i c_i
    return
end
pbo_reaction (generic function with 1 method)

Implement the calculation of the double layer charge. This just adds one additional equation to calculate the value of the integral. VoronoiFVM uses sparsity detection in order to obtain add the sparsity pattern of this operator to the overall one.

function pb_generic(y0, u0, sys, data)
    y = reshape(y0, sys)
    u = reshape(u0, sys)
    y[iQ, 1] = u[iQ, 1] - u[iq, :] ⋅ nv # implements Q=\int_\Omega q
    return
end
pb_generic (generic function with 1 method)
begin
    pbo_system = VoronoiFVM.System(
        grid;
        reaction = pbo_reaction,
        flux = pb_flux,
        bcondition = pb_bc,
        generic = pb_generic,
        unknown_storage = :sparse # use sparse storage of the unknown vector
    )

    enable_species!(pbo_system, species = [iφ, iq], regions = [1])
    # add the double layer charge as additional species at boundary 1:
    enable_boundary_species!(pbo_system, iQ, [1])
    nv .= nodevolumes(pbo_system)
end;
VoronoiFVM.has_generic_operator(pbo_system)
true

Calculate the double layer capacitance vs voltage curve

function dlcapsweep(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 = unknowns(system, inival = 0)
        else
            inival = solutions[end]
        end
        push!(solutions, solve(system; inival))
        # Here, we can directly access the double layer charge as
        # solution component
        Q = solutions[end][iQ, 1]
        bcvals[1] = v + δ
        solδ = solve(system; inival = solutions[end])
        Qδ = solδ[iQ, 1]
        push!(dlcaps, (Qδ[1] - Q[1]) / δ)
    end
    return voltages, dlcaps, TransientSolution(solutions, voltages)
end
dlcapsweep (generic function with 1 method)
pbo_volts, pbo_caps, pbo_sols = dlcapsweep(pbo_system)
(0.0:0.01:1.0, [0.2279211284512252, 0.2323215509488505, 0.24560434836520817, 0.2682802752170324, 0.3012241276640487, 0.34571356466214787, 0.40348767573309546, 0.4768306796052932, 0.5686896505780261, 0.6828408873538888  …  4.535790741564331e12, 6.694035309081421e12, 9.879227520046996e12, 1.4580015175693969e13, 2.1517557125893555e13, 3.1756157938746336e13, 4.686654535776123e13, 6.916683932620117e13, 1.020781802002832e14, 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.0 0.0 … 0.0 0.0; 0.0 NaN … NaN NaN], [0.01 0.009674563602456738 … 9.237370354343367e-6 6.465143843346327e-36; 770183.0785876448 743926.521956174 … 693.796260457962 -5.572267160095641e-42; 0.002293713725788274 NaN … NaN NaN], [0.02 0.01933737655569516 … 1.830282568784546e-5 1.280996606961598e-35; 1.658521748094714e6 1.5934423500260697e6 … 1374.680491459436 1.0980080326360017e-41; 0.004675118475748602 NaN … NaN NaN], [0.03 0.028976332547113244 … 2.703703974868103e-5 1.8922955816241115e-35; 2.801298125527285e6 2.6679554281862536e6 … 2030.6861595805808 -1.873828355746101e-40; 0.007235274526353839 NaN … NaN NaN], [0.04 0.03857861313360881 … 3.530395636764442e-5 2.4708888720557127e-35; 4.3738282129090745e6 4.115911824820652e6 … 2651.594463344157 -1.503801473677301e-39; 0.0100721266209927 NaN … NaN NaN], [0.05 0.04813033116640495 … 4.299813647931688e-5 3.009396902700902e-35; 6.617357535984294e6 6.132836264874609e6 … 3229.4861453607864 -6.394687358043636e-39; 0.013294321683177773 NaN … NaN NaN], [0.06 0.05761617795125789 … 5.004812708334442e-5 3.502818748971761e-35; 9.87607118565252e6 8.983616422027325e6 … 3758.9945290832356 -1.8422140289333155e-38; 0.017025509074486918 NaN … NaN NaN], [0.07 0.06701908305595497 … 5.6416027295178016e-5 3.9485017655699003e-35; 1.4649896059320828e7 1.3030361033735603e7 … 4237.272933352258 -4.0190527334113987e-38; 0.021409345229854687 NaN … NaN NaN], [0.08 0.07631990299437343 … 6.209404779787786e-5 4.345900785932484e-35; 2.1671195817754738e7 1.87668447237299e7 … 4663.7362709194185 -7.141991323295859e-38; 0.02661549435087813 NaN … NaN NaN], [0.09 0.08549716482188524 … 6.709916609103969e-5 4.696204048408479e-35; 3.2017124510484297e7 2.685982648686085e7 … 5039.659626573801 -1.0858635591392667e-37; 0.03284703175053636 NaN … NaN NaN]  …  [0.91 0.2723334591068045 … 9.830385042594758e-5 6.880188938832809e-35; 2.3261904016877727e21 3.8712356463144646e10 … 7383.37926678981 -6.389658871551468e-33; 1.1630952008882016e11 NaN … NaN NaN], [0.92 0.2727500284677846 … 9.831524079548552e-5 6.880986139493458e-35; 3.433050943435032e21 3.934513937035178e10 … 7384.234775789538 -6.165399150565574e-33; 1.7165254717624973e11 NaN … NaN NaN], [0.93 0.27316037783279323 … 9.832640201129167e-5 6.881767301911823e-35; 5.066583875365027e21 3.9978585882881004e10 … 7385.073073489546 -5.952038688677802e-33; 2.533291937728161e11 NaN … NaN NaN], [0.94 0.2735646843799637 … 9.83373419311011e-5 6.882532976042983e-35; 7.477393312557031e21 4.0612677614489006e10 … 7385.894750069348 -5.748908232557507e-33; 3.73869665632483e11 NaN … NaN NaN], [0.95 0.27396311798639467 … 9.834806803918683e-5 6.88328368570378e-35; 1.1035327180218686e22 4.124739691918747e10 … 7386.7003676582835 -5.555387825881395e-33; 5.517663590156324e11 NaN … NaN NaN], [0.96 0.2743558416170826 … 9.835858746880248e-5 6.884019930143553e-35; 1.6286216450586584e22 4.188272685263071e10 … 7387.490462021435 -5.370902597815067e-33; 8.143108225340942e11 NaN … NaN NaN], [0.97 0.2747430116886644 … 9.836890702300541e-5 6.884742185501539e-35; 2.4035612351468206e22 4.2518651135954315e10 … 7388.265544123866 -5.194918918836546e-33; 1.201780617578242e12 NaN … NaN NaN], [0.98 0.2751247784098808 … 9.83790331939987e-5 6.885450906160586e-35; 3.5472367867813983e22 4.3155154121883064e10 … 7389.026101581772 -5.026940945425737e-33; 1.7736183933955977e12 NaN … NaN NaN], [0.99 0.2755012861005022 … 9.838897218110978e-5 6.88614652600547e-35; 5.235102246407546e22 4.37922207629412e10 … 7389.772600015129 -4.866507654782349e-33; 2.6175511232087393e12 NaN … NaN NaN], [1.0 0.2758726734903106 … 9.839872990752456e-5 6.886829459594095e-35; 7.726097009500335e22 4.4429836581612366e10 … 7390.50548430307 -4.7131898214291714e-33; 3.8630485047552017e12 NaN … NaN NaN]], 0.0:0.01:1.0, NewtonSolverHistory[]))
pbo_sols[end]
3×52 VoronoiFVM.SparseSolutionArray{Float64, 2, Int32}:
 1.0           0.275873      0.202794     0.169966   …     9.83987e-5    6.88683e-35
 7.7261e22     4.44298e10    2.58457e9    7.20254e8     7390.51         -4.71319e-33
 3.86305e12  NaN           NaN          NaN              NaN           NaN
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_spacechargedensity(ϕ)
    sumyz = zero(ϕ)
    for i in 1:length(z)
        sumyz += z[i] * pbi_molfrac(ϕ, i)
    end
    return F * c̄ * sumyz
end
pbi_spacechargedensity (generic function with 1 method)
function pbi_reaction(y, u, node, data)
    y[iφ] = u[iq]
    y[iq] = u[iq] - pbi_spacechargedensity(u[iφ])
    return
end
pbi_reaction (generic function with 1 method)
begin

    pbi_system = VoronoiFVM.System(
        grid;
        reaction = pbi_reaction,
        flux = pb_flux,
        bcondition = pb_bc,
        generic = pb_generic,
        unkown_storage = :sparse
    )

    enable_species!(pbi_system, species = [iφ, iq], regions = [1])
    enable_boundary_species!(pbi_system, iQ, [1])

end
pbi_volts, pbi_caps, pbi_sols = dlcapsweep(pbi_system)
(0.0:0.01:1.0, [0.22792112829474656, 0.2323166562568724, 0.24558296238939673, 0.2682244335553471, 0.3011027537086486, 0.34547100866758745, 0.40302337127411797, 0.4759624294985512, 0.5670867874531943, 0.6799016647915379  …  1.6615077567161407, 1.648492447738903, 1.6351761185662284, 1.6216101924282356, 1.6078596289714753, 1.5940011670112497, 1.5801209693444562, 1.5663116858499748, 1.5526690268430343, 1.5392880156372968], TransientSolution{Any, 3, Vector{Any}, StepRangeLen{Float64, Base.TwicePrecision{Float64}, Base.TwicePrecision{Float64}, Int64}}(Any[[0.0 0.0 … 0.0 0.0; 0.0 0.0 … 0.0 0.0; 0.0 0.0 … 0.0 0.0], [0.01 0.009674565740806013 … 9.237401913675065e-6 6.4651659314094e-36; 770161.7929485585 743907.4659817506 … 693.7986307859829 4.887886501479563e-42; 0.002293697778634902 0.0 … 0.0 0.0], [0.02 0.01933739458710759 … 1.8303078174837674e-5 1.281014278275626e-35; 1.6583313857287976e6 1.5932736862664495e6 … 1374.6994550178708 -2.4669743573360575e-41; 0.004674983468875458 0.0 … 0.0 0.0], [0.03 0.02897639891080578 … 2.703789214876588e-5 1.8923552402601653e-35; 2.8005284525481714e6 2.6672851317663076e6 … 2030.7501808917302 -1.7603939703404386e-40; 0.007234774241549017 0.0 … 0.0 0.0], [0.04 0.03857879032931354 … 3.530597861058704e-5 2.471030406832481e-35; 4.371501221533081e6 4.1139336116121896e6 … 2651.7463482960966 -1.4740458217579362e-39; 0.010070777407589814 0.0 … 0.0 0.0], [0.05 0.04813073300819722 … 4.3002093324967484e-5 3.009673838401386e-35; 6.611230735607479e6 6.127788495814615e6 … 3229.783333125046 -6.309165037159845e-39; 0.013291219913180824 0.0 … 0.0 0.0], [0.06 0.05761700690912229 … 5.005498639671611e-5 3.503298825507761e-35; 9.861096225129828e6 8.971752594172677e6 … 3759.5097137046223 -1.816785039015882e-38; 0.017018994557003328 0.0 … 0.0 0.0], [0.07 0.06702069411926119 … 5.642697455746506e-5 3.9492679535937835e-35; 1.4614839448078213e7 1.3003887141167352e7 … 4238.095153689524 -3.950104810928035e-38; 0.021396388408061188 0.0 … 0.0 0.0], [0.08 0.07632291087241937 … 6.211050841066459e-5 4.3470528478866273e-35; 2.159126333078234e7 1.870988050535633e7 … 4664.972585899112 -6.993374179358427e-38; 0.02659058250018928 0.0 … 0.0 0.0], [0.09 0.08550262057379068 … 6.712283647877935e-5 4.697860715356911e-35; 3.1837919333442215e7 2.6740696820269518e7 … 5041.43745124439 -1.0564943530408122e-37; 0.0328001393477958 0.0 … 0.0 0.0]  …  [0.91 0.6216348682552673 … 0.00010614269340449226 7.428821779973225e-35; 5.355707815488365e9 5.355706892479685e9 … 7972.140186865812 -7.390556397265991e-42; 2.271748421230688 0.0 … 0.0 0.0], [0.92 0.6292533108610063 … 0.00010624974683503027 7.436314343341398e-35; 5.355707815492339e9 5.355707129328256e9 … 7980.1807762440185 4.9776824259402936e-43; 2.2882993414473707 0.0 … 0.0 0.0], [0.93 0.6368907011963554 … 0.0001063566063892804 7.443793337500575e-35; 5.355707815495032e9 5.355707305777183e9 … 7988.206804063241 -2.564120636256701e-41; 2.3047185816724305 0.0 … 0.0 0.0], [0.94 0.6445474374833199 … 0.00010646329044417426 7.451260048633966e-35; 5.355707815496857e9 5.35570743713736e9 … 7996.219650580693 1.1930529491956295e-41; 2.321003374388665 0.0 … 0.0 0.0], [0.95 0.6522238341336632 … 0.00010656980411580759 7.458714834812033e-35; 5.355707815498093e9 5.355707534859713e9 … 8004.219700056957 -8.41145412002479e-42; 2.3371515345401632 0.0 … 0.0 0.0], [0.96 0.659920103452666 … 0.00010667613931128037 7.466157129620735e-35; 5.355707815498931e9 5.355707607503623e9 … 8012.20634464818 -3.818329108636591e-42; 2.353161586684778 0.0 … 0.0 0.0], [0.97 0.6676363402950177 … 0.00010678227517388473 7.473585473319784e-35; 5.355707815499498e9 5.35570766146342e9 … 8020.178017844658 -1.5372162151885555e-41; 2.369032871630028 0.0 … 0.0 0.0], [0.98 0.6753725105302449 … 0.00010688817895392736 7.480997573790647e-35; 5.355707815499884e9 5.355707701513126e9 … 8028.132259874226 -4.617711785315587e-42; 2.384765626593753 0.0 … 0.0 0.0], [0.99 0.683128444098655 … 0.00010699380731248987 7.488390397785434e-35; 5.355707815500143e9 5.355707731214766e9 … 8036.065815628824 2.73877951142391e-42; 2.4003610334619734 0.0 … 0.0 0.0], [1.0 0.6909038332509377 … 0.00010709910803411353 7.495760290796003e-35; 5.355707815500321e9 5.355707753224251e9 … 8043.974763298014 9.625616705354513e-42; 2.41582123102185 0.0 … 0.0 0.0]], 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)[iφ, :]
    cp = c̄ * pb_molfrac.(ϕ, 1) / ufac"mol/dm^3"
    cm = c̄ * pb_molfrac.(ϕ, 2) / ufac"mol/dm^3"
    q = pbsols(pbv)[iq, :] / (F * 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)
    scalarplot!(vis[1, 2], X / nm, q, color = :green, clear = false)
    return reveal(vis)
end
plotsols (generic function with 1 method)