RasmahRasmah

Electrochemical (Nernst–Planck) Topology Optimization

println("Julia ", VERSION)
Julia 1.12.6

Overview

In an electrolyte, ions move under both diffusion and the electric field. At thermal equilibrium the ionic double layer obeys the nonlinear Poisson–Boltzmann (Nernst–Planck) equation (Newman & Thomas-Alyea 2004; Bard & Faulkner 2001)

\[-\nabla\cdot\bigl(\varepsilon(\rho)\nabla\varphi\bigr) + 2 q n_i \sinh\!\left(\frac{\varphi}{V_T}\right) = 0, \qquad \varepsilon(\rho) = \varepsilon_{min} + \rho^{p}(\varepsilon_0 - \varepsilon_{min})\]

where φ is the potential, V_T = RT/F the thermal voltage, and n_i the bulk ion density. Rasmah solves this by a dense Newton and maximizes the ionic conduction current I = Σ_anode [K(ρ)φ]ᵢ through the cell — so the design concentrates high-conductivity electrolyte along the current path between electrodes, subject to a volume budget.

using Rasmah
using LinearAlgebra

gg = LevelSetGrid(14, 14, 1 / 13, 1 / 13)
m  = grid_triangle_mesh(gg)
N  = size(m.vertices, 2)

# anode (bottom-left, φ = 0) and cathode (top-right, φ = V)
fixed = Int[]; gvals = Float64[]
for i in 1:N
    x, y = m.vertices[1, i], m.vertices[2, i]
    (x < 1e-12 && y < 0.4) && (push!(fixed, i); push!(gvals, 0.0))
    (x > 1 - 1e-12 && y > 0.6) && (push!(fixed, i); push!(gvals, 1.0))
end

nernst_planck_topopt distributes the high-permittivity electrolyte ρ to maximize the ionic current (objective = −I, so the reported history decreases as the current grows).

x, hist = nernst_planck_topopt(m, 1.0, fixed, gvals, 0.5; q = 1.0, ni = 0.0, VT = 1.0, max_iter = 40)

println("objective (−I): ", round(hist[1]; digits = 3), " → ", round(hist[end]; digits = 3))
println("final volume frac = ", round(sum(x) / length(x); digits = 3))

vtk(m; field = x, fieldname = "electrolyte density", edges = true)
objective (−I): -0.074 → -0.384
final volume frac = 0.498

Verification

The physics is anchored by the Debye length λ_D = √(ε RT / (2 F² c₀)), the screening length of the linearized double layer (a closed form from the same Nernst–Planck model the optimizer drives):

ε, F, RT, c0 = 8.854e-12, 96485.0, 2478.9, 1.0
λD = debye_length(ε, RT, F, c0)
println("Debye length = ", round(λD; sigdigits = 4), " m   (nanometre screening length)")
Debye length = 1.086e-9 m   (nanometre screening length)

Next steps