Hello all, I started using DifferentialEquations package to solve a transit convection-diffusion equation with reaction, but the result looks strange.
Firstly, I define a
du/dt = -usg du/dx + Deg d^2u/dx
the boundary condition is
u(0) = Deg/usg * du/dx + u_in
u_in = 2
du/dx(end) = 0
the code is as follow, modified from the heat equation example in DiffEqOperators.jl
using DiffEqOperators, DifferentialEquations, Plots plotly() nknots = 100 h = 0.1/(nknots+1) knots = collect(range(h, step=h, length=nknots)) ord_deriv = 2 ord_approx = 2 usg = 0.1 Deg = 2.2E-4 Δ1 = UpwindDifference(1, 2, h, nknots, -usg) Δ2 = CenteredDifference(ord_deriv, ord_approx, h, nknots) bc = RobinBC((1.0, -Deg/usg, 2.), (0.0, 1.0, 0.0),knots, 2) #bc = Dirichlet0BC(Float64) t0 = 0.0 t1 = 8 u0 = 0.3*exp.(1 .- knots) #u_analytic.(knots, t0) step(u,p,t) = Δ1*bc*u + Deg*Δ2*bc*u #- 1.2E-1*u prob = ODEProblem(step, u0, (t0, 2.)) sol = solve(prob, alg_hints = [:stiff]) surface(sol.t, knots, Array(sol), xlabel="t [s]", ylabel="L ")
The first interesting thing is that, without rteaction, u at the inlet is less than 2.
and secondly, u at the outlet increased. as can be seen in the figure below
when i replace the boundary condition by:
bc = RobinBC((1.0, 0., 2.), (0.0, 1.0, 0.0),knots, 2)
the first problem is solved, but the second problem is not. Is there anyone can give me some guidance? Thanks a lot!