basically, I would like to simulate a gas-solid reaction. u[:,1] is gas reactant, u[:,2] is gas product, and u[:3] is solid product. some times the result is strange.
the code is below
using OrdinaryDiffEq, DiffEqOperators, Plots
pyplot()
nknots = 100 # 100 knots
Lbed = 0.03 # bed length
h = Lbed/(nknots+1) #
knots = collect(range(h, step=h, length=nknots))
ord_approx = 2
usg = 0.1
Deg = 2.2E-4
# Δ1, first derivative, (1st order derivative, 2nd order approxi.)
Δ1 = UpwindDifference(1, ord_approx, h, nknots, -usg)
Δ2 = CenteredDifference(2, ord_approx, h, nknots)
bc1 = RobinBC((1.0, -Deg/usg, 1.0), (0.0, 1.0, 0.0), h, ord_approx)
bc2 = RobinBC((1.0, -Deg/usg, 0.0), (0.0, 1.0, 0.0), h, ord_approx)
function flow!(du, u,p,t)
du[:,1] = Δ1*bc1*u[:,1] + Deg*Δ2*bc1*u[:,1]
du[:,2] = Δ1*bc2*u[:,2] + Deg*Δ2*bc2*u[:,2]
end
function source_term!(du,u,p,t)
du .= 0.
for ii = 1:length(du[:,1])
if u[ii,3] <= 1.0
rA = -0.5*u[ii,1]*(1-u[ii,3])
else
rA = 0.
end
du[ii,1] = rA
du[ii,2] = -rA
du[ii,3] = -rA
end
end
u0 = zeros(length(knots),3)
u0[:,1] .= 0.01
just by running the code below, different results were obtained
probSplit = SplitODEProblem(flow!, source_term!, u0, (t0, 4.))
sol_split = solve(probSplit, KenCarp4(autodiff=:false))
plt1 = surface(sol_split.t, knots, Array(sol_split)[:,1,:], xlabel="t [s]", ylabel="L ", camera=(-45, 30), legend=:false)
plt2 = surface(sol_split.t, knots, Array(sol_split)[:,2,:], xlabel="t [s]", ylabel="L ", camera=(-45, 30), legend=:false)
plt3 = surface(sol_split.t, knots, Array(sol_split)[:,3,:], xlabel="t [s]", ylabel="L ", camera=(-45, 30), legend=:false)
plot(plt1, plt2, plt3, layout=(1,3), size = (900, 300))
a normal result
three abnormal results
error
Can anyone tell me where is the problem?