Hello All,
I am new to Julia, and I have a lot of difficulty navigating the documentation (especially with functions being imported all at once with “using …”, so I can never tell what function belongs to what).
I’m trying to model the start up (aka not steady state) of a gas phase plug flow reactor. Essentially, this boils down to solving a partial differential algebraic equation. It is a PDE because the concentrations vary in space and time, and it is a DAE because the volumetric flow rate changes along the tube. I am trying to use MethodOfLines.jl to perform finite difference on the spatial variable, and then solve the equation using a DAE solver like Rosenbrock23 or something.
The set of equations I am trying to solve are:
\frac{\partial C_A}{\partial t} = -\frac{\partial (Q C_A)}{\partial V} - 2kP_A^2
\frac{\partial C_B}{\partial t} = -\frac{\partial (Q C_B)}{\partial V} +1kP_A^2
\frac{\partial C_C}{\partial t} = -\frac{\partial (Q C_C)}{\partial V}
0= \frac{\partial (Q)}{\partial V} - (-1(R_gT)*kP_A^2)
Where P_i=R_g T C_i is the partial pressure of species i, and k is the rate constant in units of mol/(L-atm^2-hr)
When I try to run my code, I am getting a LinearAlgebra.SingularException, which I would normally take to mean there is an error in my model, but I can’t seem to find one. I see that the matrix is not singular if I have nonzero initial conditions for Ca and Cb, but the solution is still unstable and that is also not the problem I want to be solving.
One additional pain point is that I can’t really understand the documentation of MethodofLines, (especially the discretize function, it is not even in the documentation!) and the examples are not very rigorous (for example, the heat equation example provides the code, but it has extra variables that go unused like “order”, Solving the Heat Equation · MethodOfLines.jl). Does anyone have any suggestions / a guide on how to set up a PDE and DAE with method of lines? Or are there other methods you recommend to solve a DAE like this.
This equation is first order in space, so it should be solved like an ODE or using forward difference methods. It is not a boundary value problem like a 2nd order pde would be.
Some broader questions I have are:
- Should I be using @parameters vs @variables for my independent variables t and x? I’ve seen both done.
- Do I need to specify the mass matrix for this DAE? If so, where would that be input?
using OrdinaryDiffEq, ModelingToolkit, MethodOfLines, DomainSets
using Plots
#Startup (not steady state) of an isothermal / isobaric plug flow reactor
#2A -> B gas phase reaction
#C is inert gas
@parameters t x
@variables Ca(..) Cb(..) Cc(..) Qtot(..)
# variables are concentration of the species and the total flow rate.
# For method of lines example, it has t and x as parameters, but other examples they are variables.
#t is time
#x is the volume of the tubular reactor passed through, not the length
Dt = Differential(t)
Dx = Differential(x)
Rg = 0.0821 # L-atm / mol K
Temp = 300 # Temperature in Kelvin
Pres = 1.0 # pressure in atm
k1 = 0.05 # rate constant # (mol/(hr-atm^2-L))
#Feed streams
Fa = 0.1 # mol / hr
Fb = 0.0 # mol / hr
Fc = 0.5 # mol / hr
Qin = (Fa+Fb+Fc)*Rg*Temp/Pres # Total Flow Rate L/hr
Pj = Rg*Temp.*[Ca(t,x),Cb(t,x),Cc(t,x)] #Partial pressures
rate = (k1*Pj[1]^2.0) # units of mol/(L-hr)
eq = [
(Dt(Ca(t,x)) ~ -Qtot(t,x)*Dx(Ca(t,x))-Ca(t,x)*Dx(Qtot(t,x)) -2.0*rate),
(Dt(Cb(t,x)) ~ -Qtot(t,x)*Dx(Cb(t,x))-Cb(t,x)*Dx(Qtot(t,x)) +1.0*rate),
(Dt(Cc(t,x)) ~ -Qtot(t,x)*Dx(Cc(t,x))-Cc(t,x)*Dx(Qtot(t,x)) ),
(0 ~ Dx(Qtot(t,x)) - (Rg*Temp/Pres) * (-1.0*rate))] # Total number of moles is changing, need to account for that volumetrically
#Reactor is filled with inert upon startup, then it flows in the species.
bcs =[Ca(0,x) ~ 0.0,
Cb(0,x) ~ 0.0,
Cc(0,x) ~ Pres/(Rg*Temp),
Qtot(0,x) ~ Qin,
Ca(t,0) ~ Fa/Qin,
Cb(t,0) ~ Fb/Qin,
Cc(t,0) ~ Fc/Qin,
Qtot(t,0) ~ Qin]
domains = [t ∈ Interval(0.0,1.0),
x ∈ Interval(0.0, 1.0)]
@named pdesys = PDESystem(eq, bcs, domains, [t, x], [Ca(t, x),Cb(t, x),Cc(t, x), Qtot(t,x)])
# Method of lines discretization
dx = 0.1
discretization = MOLFiniteDifference([x => dx], t)
prob = discretize(pdesys,discretization)
sol = solve(prob, Rosenbrock23(), saveat=0.2)
discrete_x = sol[x]
discrete_t = sol[t]
solu = sol[Ca(t, x)]
plt = plot()
for i in 1:length(discrete_t)
plot!(discrete_x, solu[i, :], label="Numerical, t=$(discrete_t[i])")
end
display(plt)
