Hello everybody,
I am currently developing a ModelingToolkit.jl model to simulate axial flow in a tubular pipe with a chemical reaction, including axial diffusion. This system is currently an initial value problem, and I want to simulate the evolution of several variables along the length of pipe. The general molar balance for species j in the fluid mixture is written as (apologize for the screenshots as I do not know how to write equations in this forum):
Here, the variables that are in bold fond are functions of axial position in the pipe and the normal typeface denotes parameters. As you can see on the LHS, there are variables inside derivatives. I could not get Julia to automatically handle two variables inside a derivative term, so I used a product rule to expand these terms to the equation below:
The RHS term denotes a chemical reaction term A+B → C, which is an algebraic function:
The rest of the equations are presented in the simplified model shown at the end of this post. In this current model, I am running into issues with structural_simplify. Once I run the command and try to create a new ODEProblem with the simplified system, I get the following error message:
ERROR: ArgumentError: SymbolicUtils.BasicSymbolic{Real}[Ctotˍz(z), y_Aˍzz(z), y_Bˍzz(z), y_Cˍzz(z), uzˍz(z)] are missing from the variable map.
As far as I can tell, this is being caused due to the new variables structural_simplify generates during simplification. I looked at this section of the ModelingToolkit.jl FAQ, which suggests running the command
ModelingToolkit.missing_variable_defaults()
for my simplified system, unfortunately, this command makes all initial guesses for my system 0.0, even for the initial conditions I set for the original model. Not sure if I am using this incorrectly for the simplified model? The FAQ shows you can hard code a vector of initial guesses instead of zeroes, but depending on a given model the order of these variables and which variables are created may differ. Is there a way to just add the default IC values for the variables ModelingToolkit generates? I have included initial condition values for Dz(y_A), Dz(y_B), and Dz(y_C), as those are the variables that have a double derivative in the mole balance equation, yet it seems this initial condition is insufficient. Any feedback on how to properly specify the initial conditions for a model that uses structural_simplify would be appreciated!
Finally, as a related note, I am wondering about the need for initial conditions for algebraic variables. Is there a way to make ModelingToolkit calculate the initial values of algebraic variables based on the initial condition for differential variables? In this scenario for example, If I know the initial condition for Ca, T, and P, in principle the values for R_j, Ctot, and R_total are known. Is there a best practice to automatically calculate the initial conditions of these variables? Right now, I initialize all algebraic variables with 0.0 or with hard-coded equations, but perhaps there may be instabilities if these initial conditions are too far-off the values that the actual algebraic equations calculate.
Thank you for any feedback or suggestions on how to proceed! The code is shown below:
using ModelingToolkit, DifferentialEquations, NonlinearSolve, ForwardDiff, Plots
#Parameters
@parameters L ϵ D μ ρ_g dp
@parameters Rg k1ref Ea1
#variables
@variables z Ca(z) Cb(z) Cc(z) Ctot(z)
@variables y_A(z) y_B(z) y_C(z)
@variables uz(z) T(z) P(z)
@variables R_A(z) R_B(z) R_C(z) Rtotal(z) R1(z)
@variables k1(z)
#Differential Operators
Dz = Differential(z)
#Dzz = Differential(z)^2.0
@named pfr = ODESystem([
#Ideal Gas Law
P ~ Ctot * Rg * T
Ca ~ y_A * Ctot
Cb ~ y_B * Ctot
Cc ~ y_C * Ctot
#Component mass balances
# 1.0/L * Dz(Ca.*uz) - D*ϵ/(L^2.0) * Dz(Ctot.*Dz(y_A)) ~ R_A
# 1.0/L * Dz(Cb.*uz) - D*ϵ/(L^2.0) * Dz(Ctot.*Dz(y_B)) ~ R_B
# 1.0/L * Dz(Cc.*uz) - D*ϵ/(L^2.0) * Dz(Ctot.*Dz(y_C)) ~ R_C
1.0/L * (Ca.*Dz(uz) + uz.*Dz(Ca)) - D*ϵ/(L^2.0) * (Ctot.*Dz(Dz(y_A))+Dz(y_A).*Dz(Ctot)) ~ R_A
1.0/L * (Cb.*Dz(uz) + uz.*Dz(Cb)) - D*ϵ/(L^2.0) * (Ctot.*Dz(Dz(y_B))+Dz(y_B).*Dz(Ctot)) ~ R_B
1.0/L * (Cc.*Dz(uz) + uz.*Dz(Cc)) - D*ϵ/(L^2.0) * (Ctot.*Dz(Dz(y_C))+Dz(y_C).*Dz(Ctot)) ~ R_C
#Total Mass balance
#1.0/L * Dz(Ctot*uz) ~ Rtotal
1.0/L * (Ctot*Dz(uz)+ uz.*Dz(Ctot)) ~ Rtotal
#Pressure Drop
1.0/L * Dz(P) ~ - ρ_g * (1-ϵ)/ϵ^3.0 * uz^2.0/(2*dp)
#Energy balance
#Currently assume isothermal
Dz(T) ~ 0
#Kinetics Equations
# Arrhenius
k1 ~ k1ref * exp(-Ea1./Rg./T)
#rates
R1 ~ (k1.*(Ca^0.5))/(1.0 + k1.*Ca^0.5)
#Molar reaction rates
R_A ~ -1.0*R1
R_B ~ -0.5*R1
R_C ~ +1.0*R1
#Total reaction rate
Rtotal ~ R_A + R_B + R_C
])
bcs =[
Ca => 1.74,
Cb => 34.0,
Cc => 0.0,
Ctot => 34.0+1.74,
T => 517.0,
P => (34.0+1.74)*Rg*517.0,
y_A => 1.74/(34.0+1.74),
y_B => 34.0/(34.0+1.74),
y_C => 0.0,
k1 => k1ref.*exp(-Ea1/Rg/517.0),
R1 => 1E-10,
uz => 0.47,
R_A => 0.0,
R_B => 0.0,
R_C => 0.0,
Rtotal => 0.0,
Dz(y_A) => 0.0,
Dz(y_B) => 0.0,
Dz(y_C) => 0.0,
]
#Parameter map
pars = [
Rg => 8.314,
ϵ => 0.4,
Ea1 => 79.496*1000,
k1ref => 125.0E7*0.2,
L => 0.7,
dp => 0.0046,
ρ_g => 1.018,
D => 4.9E-6,
μ => 1E-8
]
Lspan = (0.0,1.0)
simpsys = structural_simplify(pfr)
prob = ODEProblem(simpsys,bcs,Lspan,pars,jac=true)
# prob = ODEProblem(simpsys,ModelingToolkit.missing_variable_defaults(simpsys),Lspan,pars,jac=true) #ODEProblem with default values for ICs
sol = solve(prob,Rodas5P(),reltol=1E-10,abstol=1E-10,isoutofdomain=(u,p,t)->any(x->x<0,u))
