Solving PDE + DAE using method of lines?

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])")

I think there is a problem with your boundary conditions. You specify the boundary condition for Ca(0,0) and Cc(0,0) twice with conflicting values at C(0,x) and C(t,0). Maybe you need to replace one of these with Neumann BCs?

That is true, although I’ve gotten a system like that to work with other software, maybe the implementation was different under the hood and it made some correction for my sloppy boundary conditions. Although even fixing the initial conditions to match the boundary conditions gives me an instability warning.

I guess the “right” way would be to use something like a moles or molar flow rate of substance instead of the concentration, but that would make the equations a lot trickier based on the nature of the rate laws and conservation equations.

I rearranged the final line of eq to be
(Dx(Qtot(t,x)) ~ (Rg*Temp/Pres) * (-1.0*rate))
and it solves. MTK runs a simplification algorithm that can reduce a DAE system to an ODE system, and I guess it had trouble performing that reduction on the original equation.

I set dx = 0.01 and got the following result. Ca looks bizarre because of the BC issue at the inlet, but it does solve.