Ferri/Ferrocyanide Reaction with diffusion in Julia – Part 2
In the post Simulating Ferri/Ferrocyanide Reaction in Julia, I described my first attempt to simulate a ferry/ferrocyanide reaction in Julia by solving a so-called split-ODE, where one ODE was generated from a PDE containing the diffusion equation. The other ODE was generated from a reaction network, describing the reaction at the electrode.
One problem is that the PDE is constructed using DiffEqOperators.jl, which is deprecated. The library is replaced by MethodOfLines.jl. My old code also did not run well in a Pluto notebook. This post describes my first attempt to simulate the reaction using MethodOfLines.jl. I based my code on the example on solving the heat equation.
The following packages are used:
using OrdinaryDiffEq using ModelingToolkit using MethodOfLines using DomainSets using Gadfly
Most of the constants used in the original post are used again:
const width = 0.01 const Δx = 0.0002 # The distance between points in the cell const Nₓ = floor(Int, width/Δx) # The amount of discrete points in the cell const A = 1e-6 const F = 96485.3329 # Faraday constant const R = 8.31446261815324 # Universal gas constant const E₀ = 0.37 # standard electrode potential const k₀ = 2.5e-4 # reaction rate constant const T = 298.15 # Temperture const α = 0.5 const λₐ = (-α*F)/(T*R) const λᵦ = ((1-α)*F)/(T*R) const Dₒ = 0.6e-9 # Diffusion coefficient for ferrocyanide const Dᵣ = 0.7e-9 # Diffusion coefficient for ferricyanide const cₘₐₓ = 1.0 # The highest initial concentration
Next I’ve created a single function, which defines the PDE, discretises it into an ODE and solves it:
function solve_reaction_diffusion(t_max, E; savelog=false, c0ₒ=0.0, c0ᵣ=cₘₐₓ) @parameters x t @variables cₒ(..) cᵣ(..) Dt = Differential(t) Dx = Differential(x) Dxx = Differential(x)^2 kb(E) = k₀ * exp(λₐ * (E - E₀)) kf(E) = k₀ * exp(λᵦ * (E - E₀)) r(x, cₒ, cᵣ, E) = (x <= Δx) * (kf(E) * cᵣ - kb(E) * cₒ) eqs = [ Dt(cₒ(x, t)) ~ Dₒ * Dxx(cₒ(x, t)) + r(x, cₒ(x, t), cᵣ(x, t), E(t)), Dt(cᵣ(x, t)) ~ Dᵣ * Dxx(cᵣ(x, t)) - r(x, cₒ(x, t), cᵣ(x, t), E(t)), ] bcs = [ cₒ(x, 0) ~ c0ₒ, cᵣ(x, 0) ~ c0ᵣ, Dx(cₒ(0, t)) ~ 0, Dx(cᵣ(0, t)) ~ 0, Dx(cₒ(width, t)) ~ 0, Dx(cᵣ(width, t)) ~ 0 ] domains = [ t ∈ Interval(0.0, t_max), x ∈ Interval(0.0, width) ] @named pdesys = PDESystem(eqs, bcs, domains,[x, t],[cₒ(x, t), cᵣ(x, t)]) discretization = MOLFiniteDifference([x => Δx], t) prob = discretize(pdesys,discretization) s = savelog ? [10^i for i ∈ 0:0.2:log10(t_max)] : (0:(t_max/1000):t_max) sol = solve(prob , BS3() , saveat=s , save_start=true , abstol=1e-10 , reltol=1e-10 , maxiters=1e7 ) sol, x, t, cₒ(x,t), cᵣ(x,t) end
This function has a few oddities. My initial idea, was to incorporate the reaction part of the equation in the boundary conditions like this:
eqs = [ Dt(cₒ(x, t)) ~ Dₒ * Dxx(cₒ(x, t)), Dt(cᵣ(x, t)) ~ Dᵣ * Dxx(cᵣ(x, t)), ] bcs = [ cₒ(x, 0) ~ c0ₒ, cᵣ(x, 0) ~ c0ᵣ, Dt(cₒ(0, t)) ~ kf(E(t)) * cᵣ(x, t) - kb(E(t)) * cₒ(x, t), Dt(cᵣ(0, t)) ~ kf(E(t)) * cₒ(x, t) - kb(E(t)) * cᵣ(x, t), Dx(cₒ(width, t)) ~ 0, Dx(cᵣ(width, t)) ~ 0 ]
But for some reason, this does not work. All the concentrations stay at their initial value. No clue why this exactly happens. But the trick to move to move the reaction part one step away from the boundary works.
The output of the solver is slightly different than the first post. In order to retrieve the data I need the symbols, so that is why I return those as well. because the output is different, all the helper functions need to be updated as well:
function current(Mₒ, Mᵣ, t, E) function redox_rate(cₒ, cᵣ, E) η = E - E₀ kₐ = k₀ * exp(λₐ * η) kᵦ = k₀ * exp(λᵦ * η) kₐ * cₒ - kᵦ * cᵣ end mₒ = @view Mₒ[2, :] mᵣ = @view Mᵣ[2, :] r = redox_rate.(mₒ, mᵣ, E.(t)) b = F * A b .* r end function plot_I_vs_E(I, E) plot(x=E, y=I , Guide.xlabel("E (V)") , Guide.ylabel("I (A)") , Guide.title("I vs E") , Geom.line(preserve_order=true) , Theme(default_color="purple") ) end function cv(t) s = 0.01 t₁ = inv(s) i = floor((t - t₁)/(2t₁)) (t < t₁) * (s * t + E₀) + (t >= t₁) * ( (i % 2 == 0) * (-s * (t - (2i+1)*t₁) + E₀ + 1.0) + (i % 2 != 0) * (s * (t - (2i+1)*t₁) + E₀ - 1.0) ) end
I’ve rewritten the function which generates the cyclic voltammetry signal, such that it does not use and if-then-else statements anymore. Note sure if this makes a performance difference. But at least it can directly be used in the PDE this way.
Solving the equations results in the same plot as before:
sol, x, t, cₒ, cᵣ = solve_reaction_diffusion(800, cv, savelog=false, c0ₒ=0.5, c0ᵣ=0.5) I = current(sol[cₒ], sol[cᵣ], sol.t, cv) E = cv.(sol.t) plot_I_vs_E(I, E)