115: 1D heterogeneous catalysis
Let $\Omega=(0,1)$, $\Gamma_1=\{0\}$, $\Gamma_2=\{1\}$ Regard a system of three species: $A,B,C$ and let $u_A=[A]$, $u_B=[B]$ and $u_C=[C]$ be their corresponding concentrations.
Species $A$ and $B$ exist in the interior of the domain, species $C$ lives a the boundary $\Gamma_1$. We assume a heterogeneous reaction scheme where $A$ reacts to $C$ and $C$ reacts to $B$:
\[\begin{aligned} A &\leftrightarrow C\\ C &\leftrightarrow B \end{aligned}\]
with reaction constants $k_{AC}^\pm$ and k_{BC}^\pm$.
In $\Omega$, both $A$ and $B$ are transported through diffusion:
\[\begin{aligned} \partial_t u_A - \nabla\cdot D_A \nabla u_A & = f_A\\ \partial_t u_B - \nabla\cdot D_B \nabla u_B & = 0\\ \end{aligned}\]
Here, $f(x)$ is a source term creating $A$. On $\Gamma_2$, we set boundary conditions
\[\begin{aligned} D_A \nabla u_A & = 0\\ u_B&=0 \end{aligned}\]
describing no normal flux for $A$ and zero concentration of $B$. On $\Gamma_1$, we use the mass action law to describe the boundary reaction and the evolution of the boundary concentration $C$. We assume that there is a limited amount of surface sites $S$ for species C, so in fact A has to react with a free surface site in order to become $C$ which reflected by the factor $1-u_C$. The same is true for $B$.
\[\begin{aligned} R_{AC}(u_A, u_C)&=k_{AC}^+ u_A(1-u_C) - k_{AC}^-u_C\\ R_{BC}(u_C, u_B)&=k_{BC}^+ u_B(1-u_C) - k_{BC}^-u_C\\ - D_A \nabla u_A + S R_{AC}(u_A, u_C)& =0 \\ - D_B \nabla u_B + S R_{BC}(u_B, u_C)& =0 \\ \partial_t C - R_{AC}(u_A, u_C) - R_{BC}(u_B, u_C) &=0 \end{aligned}\]
module Example115_HeterogeneousCatalysis1Dusing Printfusing VoronoiFVMusing ExtendableGridsusing GridVisualizeusing LinearAlgebrausing OrdinaryDiffEqRosenbrockusing SciMLBase: NoInit# Problem data structure to avoid global variablesmutable struct ProblemData D_A::Float64 # Diffusion coefficient for species A D_B::Float64 # Diffusion coefficient for species B kp_AC::Float64 # Forward reaction constant A->C km_AC::Float64 # Backward reaction constant C->A kp_BC::Float64 # Forward reaction constant B->C km_BC::Float64 # Backward reaction constant C->B S::Float64 # Surface site density iA::Int # Species index for A iB::Int # Species index for B iC::Int # Species index for Cendfunction main(; n = 10, Plotter = nothing, verbose = false, tend = 1, unknown_storage = :sparse, assembly = :edgewise, diffeq = false, switchbc = false ) h = 1.0 / convert(Float64, n) X = collect(0.0:h:1.0) N = length(X) iCat = 1 iBulk = 2 inodeCat = 1 if switchbc iCat, iBulk = iBulk, iCat inodeCat = N end grid = simplexgrid(X) # By default, \Gamma_1 at X[1] and \Gamma_2 is at X[end] # Species numbers iA = 1 iB = 2 iC = 3 # Create problem data structure with all parameters problem_data = ProblemData( 1.0, # D_A 1.0e-2, # D_B 100.0, # kp_AC 1.0, # km_AC 0.1, # kp_BC 1.0, # km_BC 0.01, # S iA, iB, iC ) # Diffusion flux for species A and B function flux!(f, u, edge, data) f[data.iA] = data.D_A * (u[data.iA, 1] - u[data.iA, 2]) f[data.iB] = data.D_B * (u[data.iB, 1] - u[data.iB, 2]) return nothing end # Storage term of species A and B function storage!(f, u, node, data) f[data.iA] = u[data.iA] f[data.iB] = u[data.iB] return nothing end # Source term for species a around 0.5 function source!(f, node, data) x1 = node[1] - 0.5 f[data.iA] = exp(-100 * x1^2) return nothing end # Reaction rate functions using data structure R_AC(u_A, u_C, data) = data.kp_AC * u_A * (1 - u_C) - data.km_AC * u_C R_BC(u_B, u_C, data) = data.kp_BC * u_B * (1 - u_C) - data.km_BC * u_C function breaction!(f, u, node, data) if node.region == iCat f[data.iA] = data.S * R_AC(u[data.iA], u[data.iC], data) f[data.iB] = data.S * R_BC(u[data.iB], u[data.iC], data) f[data.iC] = -R_BC(u[data.iB], u[data.iC], data) - R_AC(u[data.iA], u[data.iC], data) end return nothing end # This is for the term \partial_t u_C at the boundary function bstorage!(f, u, node, data) if node.region == iCat f[data.iC] = u[data.iC] end return nothing end physics = VoronoiFVM.Physics(; breaction = breaction!, bstorage = bstorage!, flux = flux!, storage = storage!, source = source!, data = problem_data ) sys = VoronoiFVM.System(grid, physics; unknown_storage = unknown_storage) # Enable species in bulk resp enable_species!(sys, problem_data.iA, [1]) enable_species!(sys, problem_data.iB, [1]) # Enable surface species enable_boundary_species!(sys, problem_data.iC, [iCat]) # Set Dirichlet bc for species B on \Gamma_2 boundary_dirichlet!(sys, problem_data.iB, iBulk, 0.0) # Initial values inival = unknowns(sys) inival .= 0.0 U = unknowns(sys) tstep = 0.01 time = 0.0 # Data to store surface concentration vs time p = GridVisualizer(; Plotter = Plotter, layout = (3, 1)) if diffeq inival = unknowns(sys, inival = 0) problem = ODEProblem(sys, inival, (0, tend)) # use fixed timesteps just for the purpose of CI odesol = solve(problem, Rosenbrock23(); initializealg = NoInit(), dt = tstep, adaptive = false) tsol = reshape(odesol, sys) else control = fixed_timesteps!(VoronoiFVM.SolverControl(), tstep) tsol = solve(sys; inival, times = [0, tend], control, verbose = verbose) end p = GridVisualizer(; Plotter = Plotter, layout = (3, 1), fast = true) for it in 1:length(tsol.t) time = tsol.t[it] scalarplot!( p[1, 1], grid, tsol[problem_data.iA, :, it]; clear = true, title = @sprintf("[A]: (%.3f,%.3f)", extrema(tsol[problem_data.iA, :, it])...) ) scalarplot!( p[2, 1], grid, tsol[problem_data.iB, :, it]; clear = true, title = @sprintf("[B]: (%.3f,%.3f)", extrema(tsol[problem_data.iB, :, it])...) ) scalarplot!( p[3, 1], tsol.t[1:it], tsol[problem_data.iC, inodeCat, 1:it]; title = @sprintf("[C]"), clear = true, show = true ) end return tsol[problem_data.iC, inodeCat, end]endusing Testfunction runtests() testval = 0.87544440641274 testvaldiffeq = 0.8757307218639448 for unknown_storage in (:sparse, :dense) for assembly in (:edgewise, :cellwise) for switchbc in (false, true) @test isapprox(main(; unknown_storage, assembly, switchbc), testval; rtol = 1.0e-12) @test isapprox(main(; diffeq = true, unknown_storage, assembly, switchbc), testvaldiffeq; rtol = 1.0e-12) end end end return nothingendendThis page was generated using Literate.jl.