OK… I’m finally trying to learn MTK by re-implementing old models from Modelica (i.e., the models work).
Now, I’m doing a manual discretization of an advection PDE/simple gas pipe. It is not super realistic, but that is not the point.
I want to find the pressure in each volume
I assume that inlet pressure at x=0 is known, and I “stupidly” assume that the fluid velocity is the same throughout the pipe.
Here is my Julia/MTK code attempt:
# Number of volumes
N = 10
#
@variables t m(t)[1:N] n(t)[1:N] md_i(t)[1:N] md_e(t)[1:N] nd_i(t)[1:N] nd_e(t)[1:N] p(t)[1:N] p_i(t)[1:N] Vd(t) p0(t) v(t)
@parameters R M A L T dV dx
Dt = Differential(t)
@register_symbolic u_p(t)
@register_symbolic u_v(t)
eqs = [ Dt(m) ~ md_i - md_e, # mass balance
m ~ n*M, # mass vs. mole
p*dV ~ n*R*T, # ideal gas law, constant temperature
p_i[1] ~ p0, # boundary condition for pressure
p_i[2:N] ~ p[1:N-1], # inlet pressure for volume 2 equals outlet pressure from vol 1, etc.
Vd ~ v*A, # volumetric flow rate vs. velocity
nd_i ~ Vd*p_i/(R*T), # flow version of ideal gas law
nd_e ~ Vd*p/(R*T), # -- " --
md_i ~ nd_i*M, # mass flow vs. molar flow
md_e ~ nd_e*M, # -- " --
p0 ~ u_p(t), # setting inlet pressure equal to a user defined time function
v ~ u_v(t)] # setting velocity equal to a user defined time function
@named pipe = ODESystem(eqs)
This results in an error message:
xes of (broadcast(~, Differential(t)((m(t))[1:10]), broadcast(-, md_i(t), md_e(t)))) not known
Stacktrace:
[1] error(s::String)
@ Base .\error.jl:35
[2] axes(A::Symbolics.Arr{Any, 1})
...
*Question: What is my error (or list of errors :-o )?
I believe the problem is that all variables in this equation are arrays, and Dt and ~ are scalar operators. If you broadcast these operators, it might work
Dt.(m) .~ md_i .- md_e
Normally, I collect(m) all array variables immediately after creating them since there are a lot of issues with them. The call to collect will transform the array variable to an array of variables.
This executes so far at least
N = 10
#
@variables t m(t)[1:N] n(t)[1:N] md_i(t)[1:N] md_e(t)[1:N] nd_i(t)[1:N] nd_e(t)[1:N] p(t)[1:N] p_i(t)[1:N] Vd(t) p0(t) v(t)
m, n, md_i, md_e, nd_i, nd_e, p, p_i = collect.((m, n, md_i, md_e, nd_i, nd_e, p, p_i))
@parameters R M A L T dV dx
Dt = Differential(t)
@register_symbolic u_p(t)
@register_symbolic u_v(t)
eqs = [ Dt.(m) .~ md_i - md_e; # mass balance
m .~ n*M; # mass vs. mole
p*dV .~ n*R*T; # ideal gas law; constant temperature
p_i[1] ~ p0; # boundary condition for pressure
p_i[2:N] .~ p[1:N-1]; # inlet pressure for volume 2 equals outlet pressure from vol 1; etc.
Vd ~ v*A; # volumetric flow rate vs. velocity
nd_i .~ Vd*p_i/(R*T); # flow version of ideal gas law
nd_e .~ Vd*p/(R*T); # -- " --
md_i .~ nd_i*M; # mass flow vs. molar flow
md_e .~ nd_e*M; # -- " --
p0 ~ u_p(t); # setting inlet pressure equal to a user defined time function
v ~ u_v(t)]
Dt.(m) .~ md_i - md_e creates an array, and ; is the array-concatenation operator. The comma is used to create a vector, so using this would create a vector of vectors instead of one long vector.
# Parameters
par = [R => 8.31, M => 28, A => 5, L => 10, dV => A*L/N, T => 300]
# Operating conditions
p0 = 1.2*p_a*ones(N)
x0 = [m .=> p0*(A*L/N)/(R*T)*M] # Specifying initial state(s)
u_p(t) = p0_const(t) # Specifying input function
u_v(t) = v_const(t) # Specifying input function
# Time span
tspan = (0.0,150) # Specifying tuple of initial and final time
… no error messages.
But the following:
prob = ODEProblem(pipe_simp, x0, tspan, par)
fails miserably.
Question: is it still a problem with array variable vs. array of variables?
