# Convert this cell to markdown in order to enable Pluto's inbuilt package managerif isdefined(Main, :PlutoRunner) using Pkg docsdir = joinpath(@__DIR__, "..", "docs") if isdir(docsdir) Pkg.activate(docsdir) end using Reviseendbegin using SciMLBase: ODEProblem, solve using OrdinaryDiffEqLowOrderRK: DP5 using OrdinaryDiffEqRosenbrock: Rosenbrock23 using OrdinaryDiffEqTsit5: Tsit5 using Catalyst using VoronoiFVM: VoronoiFVM, enable_species!, enable_boundary_species! using VoronoiFVM: ramp, boundary_dirichlet! using ExtendableGrids: simplexgrid using GridVisualize: GridVisualize, GridVisualizer, reveal, scalarplot!, gridplot, available_kwargs doplots = isdefined(Main, :PlutoRunner) if doplots using Plots: Plots, plot, theme using PlotThemes Plots.gr() Plots.theme(:dark) GridVisualize.default_plotter!(Plots) end import PlutoUI import Latexify using Test PlutoUI.TableOfContents(; depth = 4)endTowards Heterogeneous Catalysis
How to model and simulate heterogeneous catalysis with Catalyst.jl and VoronoiFVM.jl.
Mass action kinetics
Sources:
General notation: Horn/Jackson 1972
Textbook: Érdi/Tóth 1989
Assume \(j=1 … M\) reversible reactions of educt (substrate) species to product species
$$ α_1^j S_1 + α_2^j S_2 + … + α_n^jS_n \underset{k_j^+}{\stackrel{k_j^-}{\longrightleftharpoons}} β_1^jS_1 + β_2^j S_2 + … + β_n^j S_n$$
or equivalently,
$$∑_{i=1}^n α_{i}^j S_i \underset{k_j^+}{\stackrel{k_j^-}{\longrightleftharpoons}} ∑_{i=1}^n β_i^j S_i $$
The rate of these reactions depend on the concentrations \([S_i]\) of the species. Within \(Catalyst.jl\), due to consistency with the derivation from stochastic approaches, the default "combinatoric rate law" is
$$ r_j=k_j^+ ∏_{i=1}^n \frac{[S_i]^{α_i^j}}{α_i^j!} - k_j^- ∏_{i=1}^n \frac{[S_i]^{β_i^j}}{β_i^j!} $$
while it appears that in most textbooks, the rate law is
$$ r_j=k_j^+ ∏_{i=1}^n [S_i]^{α_i^j} - k_j^- ∏_{i=1}^n [S_i]^{β_i^j}.$$
See the corresponding remark in the Catalyst.jl docs and the github issue. We will stick to the secobd version which can be achieved by setting combinatoric_ratelaws to false at appropriate places.
Later in the project we will see that instead of the concentrations, we need to work with so called activities.
The numbers \(σ_i^j=α_i^j-β_i^j\) are the net stoichiometric coefficients of the system.
The rate differential equations then are (TODO: check this)
$$ ∂_t [S_i] + \sum_{j=1}^M \sigma_i^jr_j = 0$$
These assume that the reactions take place in a continuously stirred tank reactor (CSTR) which means that we have complete mixing, and species concentrations are spatially constant, and we have just one concentration value for each species at a given point of time.
Example 1: A \(\longrightleftharpoons\) B
$$\begin{aligned} A& \underset{k_1^+}{\stackrel{k_1^-}{\longrightleftharpoons}} B\\ r_1&= k_1^+ [A] - k_1^- [B]\\ \partial_t[A] &= -r_1\\ \partial_t[B] &= r_1 \end{aligned} $$
Solution via plain ODE problem using OrdinaryDiffEq.jl:
Set the parameters such that the forward reaction is faster then the backward reaction:
p1 = (k_p = 1, k_m = 0.1)(k_p = 1, k_m = 0.1)
Define an ODE function describing the right hand side of the ODE System:
function example1(du, u, p, t) (; k_p, k_m) = p r1 = k_p * u[1] - k_m * u[2] du[1] = -r1 return du[2] = r1endexample1 (generic function with 1 method)
Define some initial value:
u1_ini = [1.0, 0.0]2-element Vector{Float64}:
1.0
0.0
Create an solve the corresponding ODEProblem
prob1 = ODEProblem(example1, u1_ini, (0, 10), p1)�[38;2;86;182;194mODEProblem�[0m with uType �[38;2;86;182;194mVector{Float64}�[0m and tType �[38;2;86;182;194mInt64�[0m. In-place: �[38;2;86;182;194mtrue�[0m
Non-trivial mass matrix: �[38;2;86;182;194mfalse�[0m
timespan: (0, 10)
u0: 2-element Vector{Float64}:
1.0
0.0
sol1 = solve(prob1, DP5())retcode: Success
Interpolation: specialized 4th order "free" interpolation
t: 14-element Vector{Float64}:
0.0
0.0009990005004983772
0.08009897388807512
0.36242896847896533
0.7889034932765261
1.3826824734596057
2.0815676561276426
2.868973617116413
3.734077866058728
4.700873543432611
5.809433391512048
7.1226920180968705
8.732352513231017
10.0
u: 14-element Vector{Vector{Float64}}:
[1.0, 0.0]
[0.9990015481995943, 0.0009984518004057283]
[0.9233283474747286, 0.07667165252527142]
[0.7011010816554537, 0.29889891834454624]
[0.4726178473244975, 0.5273821526755024]
[0.28956223883648974, 0.7104377611635101]
[0.18301805092567067, 0.8169819490743292]
[0.12966409764760656, 0.8703359023523933]
[0.105885904986734, 0.8941140950132659]
[0.09608994134358326, 0.9039100586564166]
[0.09244782850888557, 0.9075521714911143]
[0.09127926304985184, 0.9087207369501481]
[0.0909785408820833, 0.9090214591179167]
[0.09092657481004998, 0.9090734251899499]
doplots && plot(sol1; size = (600, 200))
Mass conservation: adding the two reaction eqauations results in
$$ \partial_t ([A]+[B]) = 0,$$
therefore \([A]+[B]\) must be constant:
all(s -> isapprox(s, sum(u1_ini)), sum(sol1; dims = 1))true
Catalyst.@reaction_network
Catalyst.jl provides a convenient way to define a reaction network, and the resulting reaction system. So we use this to build the same system:
rn1 = @reaction_network rn1 begin @combinatoric_ratelaws false k_p, A --> B k_m, B --> Aend$$\begin{align*} \mathrm{A} &\xrightleftharpoons[k_{m}]{k_{p}} \mathrm{B} \end{align*}$$
The corresponding ODE system is:
ode1n = complete(ode_model(rn1))$$\begin{align} \frac{\mathrm{d} ~ A\left( t \right)}{\mathrm{d}t} &= - A\left( t \right) ~ \mathtt{k\_p} + B\left( t \right) ~ \mathtt{k\_m} \\ \frac{\mathrm{d} ~ B\left( t \right)}{\mathrm{d}t} &= A\left( t \right) ~ \mathtt{k\_p} - B\left( t \right) ~ \mathtt{k\_m} \end{align}$$
Catalyst.jl adds a new method to the ODEProblem constructor which allows to pass a reaction nerwork:
prob1n = ODEProblem( ode1n, Dict((unknowns(ode1n) .=> u1_ini)..., pairs(p1)...), (0, 10.0))�[38;2;86;182;194mODEProblem�[0m with uType �[38;2;86;182;194mVector{Float64}�[0m and tType �[38;2;86;182;194mFloat64�[0m. In-place: �[38;2;86;182;194mtrue�[0m
Initialization status: �[38;2;86;182;194mFULLY_DETERMINED�[0m
Non-trivial mass matrix: �[38;2;86;182;194mfalse�[0m
timespan: (0.0, 10.0)
u0: 2-element Vector{Float64}:
1.0
0.0
sol1n = solve(prob1n, Rosenbrock23())retcode: Success
Interpolation: specialized 2nd order "free" stiffness-aware interpolation
t: 32-element Vector{Float64}:
0.0
0.00011338642121160358
0.0061274109218615895
0.015580783474880556
0.03550070053557128
0.06358937273264251
0.10874252713661804
⋮
5.167165856404695
5.7943634838502325
6.554738769612671
7.516581299222496
8.811607339291228
10.0
u: 32-element Vector{Vector{Float64}}:
[1.0, 0.0]
[0.9998866206494874, 0.00011337935051261405]
[0.993893182023697, 0.006106817976302993]
[0.984551924300634, 0.015448075699365984]
[0.9651831060947237, 0.034816893905276314]
[0.9385822274798095, 0.0614177725201905]
[0.897504072677057, 0.10249592732294295]
⋮
[0.09389888951514043, 0.9061011104848595]
[0.09238689566886153, 0.9076131043311384]
[0.09153220709702888, 0.9084677929029711]
[0.09111309738379972, 0.9088869026162002]
[0.09095072811849562, 0.9090492718815043]
[0.09091907429253492, 0.909080925707465]
doplots && plot(sol1n; idxs = [rn1.A, rn1.B, rn1.A + rn1.B], size = (600, 200))
Unraveling @reaction_network
Let us try to look behind the macro voodoo - this is necessary to build networks programmatically and is behind the translation from python to Julia in CatmapInterface.jl.
It is interesting if there is a "macro - less" way to define variables, parameters and species.
@variables t$$\begin{equation} \left[ \begin{array}{c} t \\ \end{array} \right] \end{equation}$$
@parameters k_p k_m$$\begin{equation} \left[ \begin{array}{c} \mathtt{k\_p} \\ \mathtt{k\_m} \\ \end{array} \right] \end{equation}$$
@species A(t) B(t)$$\begin{equation} \left[ \begin{array}{c} A\left( t \right) \\ B\left( t \right) \\ \end{array} \right] \end{equation}$$
A reaction network can be combined from several reactions:
r1p = Reaction(k_p, [A], [B], [1], [1])k_p, A --> B
r1m = Reaction(k_m, [B], [A], [1], [1])k_m, B --> A
rn1x = complete(ReactionSystem([r1p, r1m], t; name = :example1))$$\begin{align*} \mathrm{A} &\xrightleftharpoons[k_{m}]{k_{p}} \mathrm{B} \end{align*}$$
Once we are here, the rest remains the same.
ode1x = ode_model(rn1x) |> complete$$\begin{align} \frac{\mathrm{d} ~ A\left( t \right)}{\mathrm{d}t} &= - A\left( t \right) ~ \mathtt{k\_p} + B\left( t \right) ~ \mathtt{k\_m} \\ \frac{\mathrm{d} ~ B\left( t \right)}{\mathrm{d}t} &= A\left( t \right) ~ \mathtt{k\_p} - B\left( t \right) ~ \mathtt{k\_m} \end{align}$$
prob1x = ODEProblem(ode1x, merge(Dict(unknowns(ode1x) .=> u1_ini), Dict(pairs(p1))), (0, 10.0))�[38;2;86;182;194mODEProblem�[0m with uType �[38;2;86;182;194mVector{Float64}�[0m and tType �[38;2;86;182;194mFloat64�[0m. In-place: �[38;2;86;182;194mtrue�[0m
Initialization status: �[38;2;86;182;194mFULLY_DETERMINED�[0m
Non-trivial mass matrix: �[38;2;86;182;194mfalse�[0m
timespan: (0.0, 10.0)
u0: 2-element Vector{Float64}:
1.0
0.0
sol1x = solve(prob1x, Rosenbrock23())retcode: Success
Interpolation: specialized 2nd order "free" stiffness-aware interpolation
t: 32-element Vector{Float64}:
0.0
0.00011338642121160358
0.0061274109218615895
0.015580783474880556
0.03550070053557128
0.06358937273264251
0.10874252713661804
⋮
5.167165856404695
5.7943634838502325
6.554738769612671
7.516581299222496
8.811607339291228
10.0
u: 32-element Vector{Vector{Float64}}:
[1.0, 0.0]
[0.9998866206494874, 0.00011337935051261405]
[0.993893182023697, 0.006106817976302993]
[0.984551924300634, 0.015448075699365984]
[0.9651831060947237, 0.034816893905276314]
[0.9385822274798095, 0.0614177725201905]
[0.897504072677057, 0.10249592732294295]
⋮
[0.09389888951514043, 0.9061011104848595]
[0.09238689566886153, 0.9076131043311384]
[0.09153220709702888, 0.9084677929029711]
[0.09111309738379972, 0.9088869026162002]
[0.09095072811849562, 0.9090492718815043]
[0.09091907429253492, 0.909080925707465]
doplots && plot(sol1x; size = (600, 200))
Example 2: A + 2B \(\longrightleftharpoons\) AB_2
rn2 = @reaction_network rn2 begin @combinatoric_ratelaws false k_0A, ∅ --> A k_0B, ∅ --> B (k_1p, k_1m), A + 2B <--> AB_2end$$\begin{align*} \varnothing &\xrightarrow{k_{0A}} \mathrm{A} \\ \varnothing &\xrightarrow{k_{0B}} \mathrm{B} \\ \mathrm{A} + 2 \mathrm{B} &\xrightleftharpoons[k_{1m}]{k_{1p}} \mathrm{AB_{2}} \end{align*}$$
ode1 = ode_model(rn2) |> complete$$\begin{align} \frac{\mathrm{d} ~ A\left( t \right)}{\mathrm{d}t} &= \mathtt{k\_0A} + \mathtt{AB\_2}\left( t \right) ~ \mathtt{k\_1m} - \left( B\left( t \right) \right)^{2} ~ A\left( t \right) ~ \mathtt{k\_1p} \\ \frac{\mathrm{d} ~ B\left( t \right)}{\mathrm{d}t} &= \mathtt{k\_0B} + 2 ~ \mathtt{AB\_2}\left( t \right) ~ \mathtt{k\_1m} - 2 ~ \left( B\left( t \right) \right)^{2} ~ A\left( t \right) ~ \mathtt{k\_1p} \\ \frac{\mathrm{d} ~ \mathtt{AB\_2}\left( t \right)}{\mathrm{d}t} &= - \mathtt{AB\_2}\left( t \right) ~ \mathtt{k\_1m} + \left( B\left( t \right) \right)^{2} ~ A\left( t \right) ~ \mathtt{k\_1p} \end{align}$$
p2 = (k_0A = 0.5, k_0B = 1, k_1p = 0.1, k_1m = 1.0e-5)(k_0A = 0.5, k_0B = 1, k_1p = 0.1, k_1m = 1.0e-5)
u2_ini = (A = 0, B = 0, AB_2 = 0)(A = 0, B = 0, AB_2 = 0)
prob2 = ODEProblem(rn2, Dict(pairs(u2_ini)), (0, 20.0), Dict(pairs(p2)))�[38;2;86;182;194mODEProblem�[0m with uType �[38;2;86;182;194mVector{Float64}�[0m and tType �[38;2;86;182;194mFloat64�[0m. In-place: �[38;2;86;182;194mtrue�[0m
Initialization status: �[38;2;86;182;194mFULLY_DETERMINED�[0m
Non-trivial mass matrix: �[38;2;86;182;194mfalse�[0m
timespan: (0.0, 20.0)
u0: 3-element Vector{Float64}:
0.0
0.0
0.0
sol2 = solve(prob2, Rosenbrock23())retcode: Success
Interpolation: specialized 2nd order "free" stiffness-aware interpolation
t: 39-element Vector{Float64}:
0.0
9.999999999999999e-5
0.0910612026693229
0.11221080247676292
0.19425387784990195
0.22690966636085835
0.3019872287052696
⋮
5.1543272091630765
6.020028262155715
7.249923642074722
9.245870455356135
13.607725183447293
20.0
u: 39-element Vector{Vector{Float64}}:
[0.0, 0.0, 0.0]
[4.9999999999999365e-5, 9.999999999999873e-5, 6.249999998169415e-19]
[0.04553017064322797, 0.09106034128645595, 4.306914334706892e-7]
[0.05610386037425052, 0.11220772074850104, 1.5408641309379785e-6]
[0.09711064296875035, 0.1942212859375007, 1.629595620061784e-5]
[0.11342330016992642, 0.22684660033985285, 3.153301050274919e-5]
[0.1508927569091926, 0.3017855138183852, 0.00010085744344220703]
⋮
[1.0739586115719235, 2.147917223143847, 1.503204993009615]
[1.0763324747622973, 2.1526649495245946, 1.9336816563155603]
[1.0771173447734066, 2.154234689546813, 2.547844476263955]
[1.0772455425930605, 2.154491085186121, 3.5456896850850077]
[1.0772548476080235, 2.154509695216047, 5.726607744115624]
[1.0772790566439387, 2.1545581132878775, 8.922720943356062]
doplots && plot(sol2; legend = :topleft, size = (600, 300))
Example 3: Catalysis for A + 2B \(\rightleftharpoons\) AB_2
The same reaction as in example 2, but now with a catalyst C.
The reaction between A and B takes places when A and B are bound (adsorbed) to the catalyst. So we have adsorption reactions, reactions at the catalyst, and desorption. The overall of free and bound catalyst needs to be constant over time.
rn3 = @reaction_network rn3 begin @combinatoric_ratelaws false k_0A, ∅ --> A k_0B, ∅ --> B (k_1p, k_1m), A + C <--> CA (k_2p, k_2m), B + C <--> CB (k_3p, k_3m), CA + 2CB <--> CAB2 + 2C (k_4p, k_4m), CAB2 <--> AB2 + Cend$$\begin{align*} \varnothing &\xrightarrow{k_{0A}} \mathrm{A} \\ \varnothing &\xrightarrow{k_{0B}} \mathrm{B} \\ \mathrm{A} + \mathrm{C} &\xrightleftharpoons[k_{1m}]{k_{1p}} \mathrm{CA} \\ \mathrm{B} + \mathrm{C} &\xrightleftharpoons[k_{2m}]{k_{2p}} \mathrm{CB} \\ \mathrm{CA} + 2 \mathrm{CB} &\xrightleftharpoons[k_{3m}]{k_{3p}} \mathrm{CAB2} + 2 \mathrm{C} \\ \mathrm{CAB2} &\xrightleftharpoons[k_{4m}]{k_{4p}} \mathrm{AB2} + \mathrm{C} \end{align*}$$
ode3 = ode_model(rn3) |> complete$$\begin{align} \frac{\mathrm{d} ~ A\left( t \right)}{\mathrm{d}t} &= \mathtt{k\_0A} + \mathtt{CA}\left( t \right) ~ \mathtt{k\_1m} - A\left( t \right) ~ C\left( t \right) ~ \mathtt{k\_1p} \\ \frac{\mathrm{d} ~ B\left( t \right)}{\mathrm{d}t} &= \mathtt{k\_0B} + \mathtt{CB}\left( t \right) ~ \mathtt{k\_2m} - B\left( t \right) ~ C\left( t \right) ~ \mathtt{k\_2p} \\ \frac{\mathrm{d} ~ C\left( t \right)}{\mathrm{d}t} &= \mathtt{CA}\left( t \right) ~ \mathtt{k\_1m} + \mathtt{CAB2}\left( t \right) ~ \mathtt{k\_4p} + \mathtt{CB}\left( t \right) ~ \mathtt{k\_2m} - A\left( t \right) ~ C\left( t \right) ~ \mathtt{k\_1p} - \mathtt{AB2}\left( t \right) ~ C\left( t \right) ~ \mathtt{k\_4m} - B\left( t \right) ~ C\left( t \right) ~ \mathtt{k\_2p} - 2 ~ \left( C\left( t \right) \right)^{2} ~ \mathtt{CAB2}\left( t \right) ~ \mathtt{k\_3m} + 2 ~ \left( \mathtt{CB}\left( t \right) \right)^{2} ~ \mathtt{CA}\left( t \right) ~ \mathtt{k\_3p} \\ \frac{\mathrm{d} ~ \mathtt{CA}\left( t \right)}{\mathrm{d}t} &= - \mathtt{CA}\left( t \right) ~ \mathtt{k\_1m} + A\left( t \right) ~ C\left( t \right) ~ \mathtt{k\_1p} + \left( C\left( t \right) \right)^{2} ~ \mathtt{CAB2}\left( t \right) ~ \mathtt{k\_3m} - \left( \mathtt{CB}\left( t \right) \right)^{2} ~ \mathtt{CA}\left( t \right) ~ \mathtt{k\_3p} \\ \frac{\mathrm{d} ~ \mathtt{CB}\left( t \right)}{\mathrm{d}t} &= - \mathtt{CB}\left( t \right) ~ \mathtt{k\_2m} + B\left( t \right) ~ C\left( t \right) ~ \mathtt{k\_2p} + 2 ~ \left( C\left( t \right) \right)^{2} ~ \mathtt{CAB2}\left( t \right) ~ \mathtt{k\_3m} - 2 ~ \left( \mathtt{CB}\left( t \right) \right)^{2} ~ \mathtt{CA}\left( t \right) ~ \mathtt{k\_3p} \\ \frac{\mathrm{d} ~ \mathtt{CAB2}\left( t \right)}{\mathrm{d}t} &= - \mathtt{CAB2}\left( t \right) ~ \mathtt{k\_4p} + \mathtt{AB2}\left( t \right) ~ C\left( t \right) ~ \mathtt{k\_4m} - \left( C\left( t \right) \right)^{2} ~ \mathtt{CAB2}\left( t \right) ~ \mathtt{k\_3m} + \left( \mathtt{CB}\left( t \right) \right)^{2} ~ \mathtt{CA}\left( t \right) ~ \mathtt{k\_3p} \\ \frac{\mathrm{d} ~ \mathtt{AB2}\left( t \right)}{\mathrm{d}t} &= \mathtt{CAB2}\left( t \right) ~ \mathtt{k\_4p} - \mathtt{AB2}\left( t \right) ~ C\left( t \right) ~ \mathtt{k\_4m} \end{align}$$
p3 = ( k_0A = 0.5, k_0B = 1, k_1p = 10, k_1m = 0.1, k_2p = 10, k_2m = 0.1, k_3p = 10, k_3m = 0.1, k_4p = 10, k_4m = 0.1,)(k_0A = 0.5, k_0B = 1, k_1p = 10, k_1m = 0.1, k_2p = 10, k_2m = 0.1, k_3p = 10, k_3m = 0.1, k_4p = 10, k_4m = 0.1)
Cini = 4040
u3_ini = (A = 0, B = 0, CA = 0, CB = 0, CAB2 = 0, AB2 = 0, C = Cini)(A = 0, B = 0, CA = 0, CB = 0, CAB2 = 0, AB2 = 0, C = 40)
t3end = 200200
prob3 = ODEProblem(rn3, Dict(pairs(u3_ini)), (0, t3end), Dict(pairs(p3)))�[38;2;86;182;194mODEProblem�[0m with uType �[38;2;86;182;194mVector{Float64}�[0m and tType �[38;2;86;182;194mInt64�[0m. In-place: �[38;2;86;182;194mtrue�[0m
Initialization status: �[38;2;86;182;194mFULLY_DETERMINED�[0m
Non-trivial mass matrix: �[38;2;86;182;194mfalse�[0m
timespan: (0, 200)
u0: 7-element Vector{Float64}:
0.0
0.0
40.0
0.0
0.0
0.0
0.0
sol3 = solve(prob3, Rosenbrock23())retcode: Success
Interpolation: specialized 2nd order "free" stiffness-aware interpolation
t: 146-element Vector{Float64}:
0.0
6.467843637174543e-6
0.00012098643893717323
0.0002355050342371719
0.00046854144835696093
0.0007020654986982254
0.0010683942019291256
⋮
148.43142335868797
158.44482045107486
169.05108849030577
180.27585135780103
192.14552642316042
200.0
u: 146-element Vector{Vector{Float64}}:
[0.0, 0.0, 40.0, 0.0, 0.0, 0.0, 0.0]
[3.229742998125995e-6, 6.45948599625199e-6, 39.99999998746354, 4.178820461276082e-9, 8.357640922552164e-9, 4.746544394469828e-31, 8.99172643162884e-36]
[5.9057434677841024e-5, 0.00011811486935568205, 39.999995692645626, 1.4357847907455922e-6, 2.871569581491185e-6, 4.324856087335027e-22, 1.4505871228135646e-25]
[0.00011238531367665834, 0.00022477062735331668, 39.99998389838967, 5.367203441927478e-6, 1.0734406883854957e-5, 1.467137581606062e-19, 5.959854713357241e-23]
[0.00021366983615732585, 0.0004273396723146517, 39.999938197335936, 2.0600888021138075e-5, 4.120177604227615e-5, 1.6535953182237662e-17, 1.3764788284957112e-20]
[0.000306121593402705, 0.00061224318680541, 39.99986526653216, 4.491115594608782e-5, 8.982231189217564e-5, 3.1957839669763537e-16, 3.357848527050394e-19]
[0.0004348816015929623, 0.0008697632031859246, 39.9997020535019, 9.931549936625078e-5, 0.00019863099873250155, 5.340888819921991e-15, 8.841053211671128e-18]
⋮
[0.0035692805812069157, 0.007138561162414221, 20.610384484010112, 2.3571545006330097, 4.714309001266831, 12.318152014084903, 59.53683588404366]
[0.003663085165435014, 0.007326170330871041, 20.06379957895806, 2.350299659658327, 4.70059931931873, 12.885301442059736, 63.98314603865246]
[0.0037611366797918626, 0.007522273359585107, 19.51639613498695, 2.3411713119880084, 4.682342623978709, 13.460089929040583, 68.72052186744305]
[0.003863598397665346, 0.007727196795331786, 18.96936861411051, 2.3298187873635094, 4.6596375747290875, 14.041175023791782, 73.76306826934629]
[0.003970639073546245, 0.00794127814709347, 18.423885683654106, 2.316304667607705, 4.632609335217225, 14.62720031351653, 79.12528759138101]
[0.004040742010090244, 0.008081484020181053, 18.081339221610328, 2.306639073240135, 4.613278146481272, 14.998743558663831, 82.69057662608475]
doplots && plot(sol3; legend = :topleft, size = (600, 300))
ctotal = rn3.C + rn3.CA + rn3.CB + rn3.CAB2$$\begin{equation} C\left( t \right) + \mathtt{CA}\left( t \right) + \mathtt{CAB2}\left( t \right) + \mathtt{CB}\left( t \right) \end{equation}$$
sol3[ctotal]146-element Vector{Float64}:
40.0
40.0
40.0
40.0
40.0
40.0
40.00000000000001
⋮
39.999999999994856
39.999999999994856
39.999999999994245
39.99999999999489
39.999999999995566
39.999999999995566
@test sol3[ctotal] ≈ fill(Cini, length(sol3.t))�[32m�[1mTest Passed�[22m�[39m
Example 4: Heterogeneous catalysis
Heterogeneous catalysis assumes that the catalytic reaction takes place at surface. This means that reacting species need to be transported towards or away from the surface, and one has to model coupled transport and surface reaction.
Here we use VoronoiFVM.jl to model transport and Catalyst.jl to create the surface reaction network.
Problem specification
Assume \(\Omega=(0,1)\) where a catalytic reaction takes place at \(x=0\). We assume that the educts A, B, and the product AB2 are bulk species transported to the domain. At \(x=1\) we set Dirichlet boundary conditions providing A,B and removing AB2.
A, B can adsorb at the catalyst at \(x=0\) and react to AB2 while adsorbed. The product desorbs and is released to the bulk. So we have
Mass transport in the interior of \(\Omega\):
$$\begin{aligned} \partial_t c_A + \nabla \cdot D_A \nabla c_A &=0\\ \partial_t c_B + \nabla \cdot D_B \nabla c_B &=0\\ \partial_t c_{AB2} + \nabla \cdot D_{AB2} \nabla c_{AB2} &=0 \end{aligned}$$
Coupled nonlinear robin boundary conditions at \(x=0\):
$$\begin{aligned} D_A\partial_n c_A + r_1 &= 0\\ D_B\partial_n c_A + r_2 &= 0\\ D_{AB2}\partial_n c_{AB2} - r_4 &= 0\\ \end{aligned}$$
\(r_1, r_2\) and \(r_4\) are asorption/desorption reactions:
$$\begin{aligned} r_1&=k_{1p}c_A c_C - k_{1m}c_{CA}\\ r_2&=k_{2p}c_B c_C - k_{2m}c_{CB}\\ r_4&=k_{4p}c_{AB2} - k_{4m}c_{C_C}c_{C_{AB2}}\\ \end{aligned}$$
The free catalyst sites C and the catalyst coverages CA, CB, CAB2 behave according to:
$$\begin{equation} \frac{\mathrm{d} \cdot C\left( t \right)}{\mathrm{d}t} = \mathrm{\mathtt{CA}}\left( t \right) \cdot \mathtt{k_{1m}} + \mathrm{\mathtt{CAB2}}\left( t \right) \cdot \mathtt{k_{4p}} + \mathrm{\mathtt{CB}}\left( t \right) \cdot \mathtt{k_{2m}} - A\left( t \right) \cdot C\left( t \right) \cdot \mathtt{k_{1p}} - \mathrm{\mathtt{AB2}}\left( t \right) \cdot C\left( t \right) \cdot \mathtt{k_{4m}} - B\left( t \right) \cdot C\left( t \right) \cdot \mathtt{k_{2p}} - 2 \cdot \left( C\left( t \right) \right)^{2} \cdot \mathrm{\mathtt{CAB2}}\left( t \right) \cdot \mathtt{k_{3m}} + 2 \cdot \left( \mathrm{\mathtt{CB}}\left( t \right) \right)^{2} \cdot \mathrm{\mathtt{CA}}\left( t \right) \cdot \mathtt{k_{3p}} \end{equation}$$
$$\begin{equation} \frac{\mathrm{d} \cdot \mathrm{\mathtt{CA}}\left( t \right)}{\mathrm{d}t} = - \mathrm{\mathtt{CA}}\left( t \right) \cdot \mathtt{k_{1m}} + A\left( t \right) \cdot C\left( t \right) \cdot \mathtt{k_{1p}} + \left( C\left( t \right) \right)^{2} \cdot \mathrm{\mathtt{CAB2}}\left( t \right) \cdot \mathtt{k_{3m}} - \left( \mathrm{\mathtt{CB}}\left( t \right) \right)^{2} \cdot \mathrm{\mathtt{CA}}\left( t \right) \cdot \mathtt{k_{3p}} \end{equation}$$
$$\begin{equation} \frac{\mathrm{d} \cdot \mathrm{\mathtt{CB}}\left( t \right)}{\mathrm{d}t} = - \mathrm{\mathtt{CB}}\left( t \right) \cdot \mathtt{k_{2m}} + B\left( t \right) \cdot C\left( t \right) \cdot \mathtt{k_{2p}} + 2 \cdot \left( C\left( t \right) \right)^{2} \cdot \mathrm{\mathtt{CAB2}}\left( t \right) \cdot \mathtt{k_{3m}} - 2 \cdot \left( \mathrm{\mathtt{CB}}\left( t \right) \right)^{2} \cdot \mathrm{\mathtt{CA}}\left( t \right) \cdot \mathtt{k_{3p}} \end{equation}$$
$$\begin{equation} \frac{\mathrm{d} \cdot \mathrm{\mathtt{CAB2}}\left( t \right)}{\mathrm{d}t} = - \mathrm{\mathtt{CAB2}}\left( t \right) \cdot \mathtt{k_{4p}} + \mathrm{\mathtt{AB2}}\left( t \right) \cdot C\left( t \right) \cdot \mathtt{k_{4m}} - \left( C\left( t \right) \right)^{2} \cdot \mathrm{\mathtt{CAB2}}\left( t \right) \cdot \mathtt{k_{3m}} + \left( \mathrm{\mathtt{CB}}\left( t \right) \right)^{2} \cdot \mathrm{\mathtt{CA}}\left( t \right) \cdot \mathtt{k_{3p}} \end{equation}$$
Dirichlet boundary conditions at \(x=1\) :
$$\begin{aligned} c_A&=1\\ c_B&=1\\ c_{AB2}&=0 \end{aligned}$$
Finally, we set all initial concentrations to zero besides of the catalyst concenration (number of catalyst sites) \(c_C|_{t=0}=C_0=1\).
Implementation
Surface reaction network
Define a reaction network under the assumption that the supply of A and B comes from the transport and does not need to be specified.
rnv = @reaction_network rnv begin @combinatoric_ratelaws false (k_1p, k_1m), A + C <--> CA (k_2p, k_2m), B + C <--> CB (k_3p, k_3m), CA + 2CB <--> CAB2 + 2C (k_4p, k_4m), CAB2 <--> AB2 + Cend$$\begin{align*} \mathrm{A} + \mathrm{C} &\xrightleftharpoons[k_{1m}]{k_{1p}} \mathrm{CA} \\ \mathrm{B} + \mathrm{C} &\xrightleftharpoons[k_{2m}]{k_{2p}} \mathrm{CB} \\ \mathrm{CA} + 2 \mathrm{CB} &\xrightleftharpoons[k_{3m}]{k_{3p}} \mathrm{CAB2} + 2 \mathrm{C} \\ \mathrm{CAB2} &\xrightleftharpoons[k_{4m}]{k_{4p}} \mathrm{AB2} + \mathrm{C} \end{align*}$$
odesys = ode_model(rnv)$$\begin{align} \frac{\mathrm{d} ~ A\left( t \right)}{\mathrm{d}t} &= \mathtt{CA}\left( t \right) ~ \mathtt{k\_1m} - A\left( t \right) ~ C\left( t \right) ~ \mathtt{k\_1p} \\ \frac{\mathrm{d} ~ C\left( t \right)}{\mathrm{d}t} &= \mathtt{CA}\left( t \right) ~ \mathtt{k\_1m} + \mathtt{CAB2}\left( t \right) ~ \mathtt{k\_4p} + \mathtt{CB}\left( t \right) ~ \mathtt{k\_2m} - A\left( t \right) ~ C\left( t \right) ~ \mathtt{k\_1p} - \mathtt{AB2}\left( t \right) ~ C\left( t \right) ~ \mathtt{k\_4m} - B\left( t \right) ~ C\left( t \right) ~ \mathtt{k\_2p} - 2 ~ \left( C\left( t \right) \right)^{2} ~ \mathtt{CAB2}\left( t \right) ~ \mathtt{k\_3m} + 2 ~ \left( \mathtt{CB}\left( t \right) \right)^{2} ~ \mathtt{CA}\left( t \right) ~ \mathtt{k\_3p} \\ \frac{\mathrm{d} ~ \mathtt{CA}\left( t \right)}{\mathrm{d}t} &= - \mathtt{CA}\left( t \right) ~ \mathtt{k\_1m} + A\left( t \right) ~ C\left( t \right) ~ \mathtt{k\_1p} + \left( C\left( t \right) \right)^{2} ~ \mathtt{CAB2}\left( t \right) ~ \mathtt{k\_3m} - \left( \mathtt{CB}\left( t \right) \right)^{2} ~ \mathtt{CA}\left( t \right) ~ \mathtt{k\_3p} \\ \frac{\mathrm{d} ~ B\left( t \right)}{\mathrm{d}t} &= \mathtt{CB}\left( t \right) ~ \mathtt{k\_2m} - B\left( t \right) ~ C\left( t \right) ~ \mathtt{k\_2p} \\ \frac{\mathrm{d} ~ \mathtt{CB}\left( t \right)}{\mathrm{d}t} &= - \mathtt{CB}\left( t \right) ~ \mathtt{k\_2m} + B\left( t \right) ~ C\left( t \right) ~ \mathtt{k\_2p} + 2 ~ \left( C\left( t \right) \right)^{2} ~ \mathtt{CAB2}\left( t \right) ~ \mathtt{k\_3m} - 2 ~ \left( \mathtt{CB}\left( t \right) \right)^{2} ~ \mathtt{CA}\left( t \right) ~ \mathtt{k\_3p} \\ \frac{\mathrm{d} ~ \mathtt{CAB2}\left( t \right)}{\mathrm{d}t} &= - \mathtt{CAB2}\left( t \right) ~ \mathtt{k\_4p} + \mathtt{AB2}\left( t \right) ~ C\left( t \right) ~ \mathtt{k\_4m} - \left( C\left( t \right) \right)^{2} ~ \mathtt{CAB2}\left( t \right) ~ \mathtt{k\_3m} + \left( \mathtt{CB}\left( t \right) \right)^{2} ~ \mathtt{CA}\left( t \right) ~ \mathtt{k\_3p} \\ \frac{\mathrm{d} ~ \mathtt{AB2}\left( t \right)}{\mathrm{d}t} &= \mathtt{CAB2}\left( t \right) ~ \mathtt{k\_4p} - \mathtt{AB2}\left( t \right) ~ C\left( t \right) ~ \mathtt{k\_4m} \end{align}$$
For coupling with VoronoiFVM we need species numbers which need to correspond to the species in our network:
begin smap = speciesmap(rnv) const iA = smap[rnv.A] const iB = smap[rnv.B] const iC = smap[rnv.C] const iCA = smap[rnv.CA] const iCB = smap[rnv.CB] const iCAB2 = smap[rnv.CAB2] const iAB2 = smap[rnv.AB2]end;Grid:
grid = simplexgrid(0:0.01:1)ExtendableGrids.ExtendableGrid{Float64, Int32}
dim = 1
nnodes = 101
ncells = 100
nbfaces = 2
gridplot(grid, size = (600, 100))
The grid has two boundary regions: region 1 at x=0 and region 2 at x=1.
Reaction parameters:
pcat = ( k_1p = 50, k_1m = 0.1, k_2p = 50, k_2m = 0.1, k_3p = 10, k_3m = 0.1, k_4p = 50, k_4m = 0.1,)(k_1p = 50, k_1m = 0.1, k_2p = 50, k_2m = 0.1, k_3p = 10, k_3m = 0.1, k_4p = 50, k_4m = 0.1)
Parameters for the VoronoiFVM system:
params = ( D_A = 1.0, D_B = 1.0, D_AB2 = 1.0, pcat = pcat,)(D_A = 1.0, D_B = 1.0, D_AB2 = 1.0, pcat = (k_1p = 50, k_1m = 0.1, k_2p = 50, k_2m = 0.1, k_3p = 10, k_3m = 0.1, k_4p = 50, k_4m = 0.1))
Initial values for the reaction network (needed only for the definition of the ODE problem)
C0 = 1.01.0
uv_ini = (A = 0, B = 0, CA = 0, CB = 0, CAB2 = 0, AB2 = 0, C = C0)(A = 0, B = 0, CA = 0, CB = 0, CAB2 = 0, AB2 = 0, C = 1.0)
tvend = 200.0200.0
const probv = ODEProblem(rnv, Dict(pairs(uv_ini)), (0, tvend), Dict(pairs(pcat)))�[38;2;86;182;194mODEProblem�[0m with uType �[38;2;86;182;194mVector{Float64}�[0m and tType �[38;2;86;182;194mFloat64�[0m. In-place: �[38;2;86;182;194mtrue�[0m
Initialization status: �[38;2;86;182;194mFULLY_DETERMINED�[0m
Non-trivial mass matrix: �[38;2;86;182;194mfalse�[0m
timespan: (0.0, 200.0)
u0: 7-element Vector{Float64}:
0.0
1.0
0.0
0.0
0.0
0.0
0.0
Callback functions for VoronoiFVM
First, define flux and storage functions for the bulk process:
function storage(y, u, node, p) y[iA] = u[iA] y[iB] = u[iB] return y[iAB2] = u[iAB2]endstorage (generic function with 1 method)
function flux(y, u, edge, p) (; D_A, D_B, D_AB2) = p y[iA] = D_A * (u[iA, 1] - u[iA, 2]) y[iB] = D_B * (u[iB, 1] - u[iB, 2]) return y[iAB2] = D_A * (u[iAB2, 1] - u[iAB2, 2])endflux (generic function with 1 method)
Storage term for the surface reaction:
function bstorage(y, u, bnode, p) y[iC] = u[iC] y[iCA] = u[iCA] y[iCB] = u[iCB] return y[iCAB2] = u[iCAB2]endbstorage (generic function with 1 method)
Catalytic reaction. Here we use the right hand side function of the ODE problem generated above. In VoronoiFVM, reaction term are a the left hand side, so we need to multiply by -1.
Note that we need to pass the parameter record as generated for the ODE problem instead of pcat.
function catreaction(f, u, bnode, p) probv.f(f, u, probv.p, bnode.time) for i in 1:length(f) f[i] = -f[i] end returnendcatreaction (generic function with 1 method)
Define the Dirichlet boundary condition at x=1 (region 2):
function bulkbc(f, u, bnode, p) v = ramp(bnode.time; du = (0.0, 1.0), dt = (0.0, 0.01)) boundary_dirichlet!(f, u, bnode; species = iA, value = v, region = 2) boundary_dirichlet!(f, u, bnode; species = iB, value = v, region = 2) return boundary_dirichlet!(f, u, bnode; species = iAB2, value = 0, region = 2)endbulkbc (generic function with 1 method)
Dispatch the boundary conditions
function breaction(f, u, bnode, p) return if bnode.region == 1 catreaction(f, u, bnode, p) else bulkbc(f, u, bnode, p) endendbreaction (generic function with 1 method)
Coupled transport-reaction system
Define a VoronoiFVM system from grid, params and the callback functions and enable the bulk and boundary species. unknown_storage = :sparse means that the solution is stored as a nspecies x nnodes sparse matrix in order to take into account that the surface species are non-existent in the bulk. unknown_storage = :dense would store a full matrix and solve dummy equations for the surface species values in the bulk.
begin sys = VoronoiFVM.System( grid; data = params, flux, breaction, bstorage, storage, unknown_storage = :sparse ) enable_species!(sys, iA, [1]) enable_species!(sys, iB, [1]) enable_species!(sys, iAB2, [1]) enable_boundary_species!(sys, iC, [1]) enable_boundary_species!(sys, iCA, [1]) enable_boundary_species!(sys, iCB, [1]) enable_boundary_species!(sys, iCAB2, [1])end;Define an initial value for sys:
begin u0 = VoronoiFVM.unknowns(sys; inival = 0) u0[iC, 1] = C0end;Solution
Solve the time evolution
tsol = solve(sys; inival = u0, times = (1.0e-4, tvend));t:
let t_plot = round(10^log_t_plot; sigdigits = 3) vis = GridVisualizer(; size = (600, 300), flimits = (0, 1), title = "Bulk concentrations: t=$t_plot", legend = :lt) sol = tsol(t_plot) scalarplot!(vis, grid, sol[iA, :]; color = :red, label = "A") scalarplot!(vis, grid, sol[iB, :]; color = :green, label = "B", clear = false) scalarplot!(vis, grid, sol[iAB2, :]; color = :blue, label = "AB2", clear = false) reveal(vis)end
Ctotalv = tsol[iC, 1, :] + tsol[iCA, 1, :] + tsol[iCB, 1, :] + tsol[iCAB2, 1, :]248-element Vector{Float64}:
1.0
1.0
0.9999999999999999
1.0
1.0
0.9999999999999999
0.9999999999999999
⋮
1.0000000000000004
1.0000000000000004
1.0000000000000004
1.0000000000000004
1.0000000000000004
1.0000000000000004
@test Ctotalv ≈ ones(length(tsol.t))�[32m�[1mTest Passed�[22m�[39m
let vis = GridVisualizer(; size = (600, 300), xlabel = "t", flimits = (0, 1), xlimits = (1.0e-3, tvend), legend = :lt, title = "Concentrations at x=0", xscale = :log10 ) t = tsol.t scalarplot!(vis, t, tsol[iA, 1, :]; color = :darkred, label = "A") scalarplot!(vis, t, tsol[iCA, 1, :]; color = :red, label = "CA") scalarplot!(vis, t, tsol[iB, 1, :]; color = :darkgreen, label = "B") scalarplot!(vis, t, tsol[iCB, 1, :]; color = :green, label = "CB") scalarplot!(vis, t, tsol[iAB2, 1, :]; color = :darkblue, label = "AB2") scalarplot!(vis, t, tsol[iCAB2, 1, :]; color = :blue, label = "CAB2") scalarplot!(vis, t, tsol[iC, 1, :] / C0; color = :orange, label = "C/C0") scalarplot!(vis, t, Ctotalv / C0; color = :darkorange, label = "Ctot/C0") reveal(vis)end