I’m trying to understand ModelingToolkit (MTK) connectors, components, and now Stream variables in connector.
As an exercise for myself, I am looking at 2 substance, ideal solution liquid tank systems:
I consider 3 cases:
- A single tank case, blue ellipse
- A two tank case, the green ellipse
- A three tank case; all tanks
[m is mass, \dot{m} is mass flow rate — md in code, h is level, p is pressure, x is mass fraction, A is cross sectional area, \rho is density — rho in code.]
I have simulated a “flat” model without using components for (i) the single tank case, and (ii) the three tank case. There is no problem solving them.
I have then created a connector (FluidPort) and components and simulated Case 1 (dashed lines below), which exactly matches the flat model of Case 1 (pale, solid lines below).
I have also used components to simulate Case 2; no problems.
Finally, I have used components to simulate Case 3 (which works with the flat model), but for Case 3 with components, I get an error message:
MethodError: no method matching instream(::Float64)
The function `instream` exists, but no method is defined for this combination of argument types.
Question: Why do I get an error message in the 3 tank case? What am I doing wrong?
Building blocks
Basic code: package import + components
Package import:
using Unitful
using ModelingToolkit
using ModelingToolkit: t_nounits as t, D_nounits as D
using OrdinaryDiffEq
NB: the components are not meant to be production quality!! This is a self-learning project. I have specifically hard-coded the components to be valid for 2 substances.
FluidPort connector:
@connector function FluidPort(; name)
vars = @variables begin
p(t), [guess = 0.0, description="Pressure, Pa"]
md(t), [connect = Flow, guess = 0.0, description="Mass flow, kg/s"]
x1(t), [connect = Stream, guess = 1/2, description="Mass fraction 1, -"]
x2(t), [connect = Stream, guess = 1/2, description="Mass fraction 2, -"]
end
System(Equation[], t, vars, []; name = name)
end
One port sources/sinks
PressureBoundary: specifying pressure and mass fractions – may supply/draw mass, e.g., through a valve.
@component function PressureBoundary(; name, p_set=101325.0, x_set=[0.5, 0.5])
@named a = FluidPort()
x_s = sum(x_set)
x_s > 0 || error("PressureBoundary $name: Sum of x_set must be greater than zero")
x_p = x_set ./ x_s
pars = @parameters begin
p = p_set, [unit = u"Pa"]
x1 = x_p[1], [unit = u"kg/kg"]
x2 = x_p[2], [unit = u"kg/kg"]
end
eqs = [
a.p ~ p
a.x1 ~ x1
a.x2 ~ x2
]
compose(System(eqs, t, [], pars; name = name), a)
end
PressureReference: same as PressureBoundary, but with mass flow forced to zero – used to specify pressure.
@component function PressureReference(; name, p_set=101325.0, x_set=[0.5, 0.5])
@named a = FluidPort()
x_s = sum(x_set)
x_s > 0 || error("PressureBoundary $name: Sum of x_set must be greater than zero")
x_p = x_set ./ x_s
pars = @parameters begin
p = p_set, [unit = u"Pa"]
x1 = x_p[1], [unit = u"kg/kg"]
x2 = x_p[2], [unit = u"kg/kg"]
end
eqs = [
a.p ~ p,
a.md ~ 0, # THE KEY DIFFERENCE from PressureBoundary: Blocks mass flow
a.x1 ~ x1,
a.x2 ~ x2
]
compose(System(eqs, t, [], pars; name = name), a)
end
MassFlowPump: an ideal pump (no pressure dependence) to provide specified total mass flow rate and specified mass fractions.
@component function MassFlowPump(; name, md_set=1.0, x_set=[0.5, 0.5])
@named a = FluidPort()
x_s = sum(x_set)
x_s > 0 || error("MassFlowPump $name: Sum of x_set must be greater than zero")
x_p = x_set ./ x_s
pars = @parameters begin
md = md_set, [unit=u"kg/s", description = "Specified mass flow rate"]
x1 = x_p[1], [unit = u"kg/kg", description = "Specified mass fractions"]
x2 = x_p[2], [unit = u"kg/kg", description = "Specified mass fractions"]
end
eqs = [
# Inflow to system means outflow from component (negative md)
a.md ~ -md,
# Define the stream variables for the port
a.x1 ~ x1,
a.x2 ~ x2
]
compose(System(eqs, t, [], pars; name = name), a)
end
Two port components
Valve: essentially a turbulent flow valve, with flow proportional to the square root of the pressure difference.
@component function Valve(; name, Kv=1.0)
@named a = FluidPort()
@named b = FluidPort()
pars = @parameters begin
# Units must be (mass flow) / (pressure^0.5)
K = Kv, [description = "Flow coefficient"]
p_eps = 1e-4, [description = "Regularization factor"]
end
vars = @variables begin
Dp(t), [unit = u"Pa", description = "Pressure drop a-b"]
end
eqs = [
# 1. Pressure drop
Dp ~ a.p - b.p,
# 2. Flow Equation: Square Root Law (turbulent flow)
# - use 'p_eps' here for stability
a.md ~ K * sign(Dp) * sqrt(abs(Dp) + p_eps),
# 3. Continuity: No storage
a.md + b.md ~ 0,
# 4. Stream Propagation: same composition on both sides
a.x1 ~ instream(b.x1),
a.x2 ~ instream(b.x2),
b.x1 ~ instream(a.x1),
b.x2 ~ instream(a.x2)
]
compose(System(eqs, t, vars, pars; name = name), a, b)
end
Tank: simple liquid tank with specified cross sectional area, and unlimited possible liquid level. The pressure of gas above liquid is specified (PressureReference) where also liquid may be pumped in (port a) and the pressure at the bottom (port b) is this gas level + the hydrostatic pressure of the liquid.
@component function Tank(; name, Ap=1.0, rho_p=1000.0, m0=[1,1])
@named a = FluidPort() # Top Port (Headspace)
@named b = FluidPort() # Bottom Port (Drain)
# Gravity from Unitful
gp = ustrip(u"m/s^2", Unitful.gn)
# Symbolic parameters
pars = @parameters begin
A = Ap, [unit = u"m^2"]
rho = rho_p, [unit = u"kg/m^3"]
g = gp, [unit = u"m/s^2"]
end
# Symbolic variables
vars = @variables begin
m1(t) = m0[1], [unit = u"kg"]
m2(t) = m0[2], [unit = u"kg"]
h(t), [unit = u"m"]
x1(t), [unit = u"kg/kg"]
x2(t), [unit = u"kg/kg"]
x1a(t), [unit = u"kg/kg"]
x2a(t), [unit = u"kg/kg"]
x1b(t), [unit = u"kg/kg"]
x2b(t), [unit = u"kg/kg"]
end
eqs = [
#
# --- Level, assuming incompressible, ideal solution liquid
h ~ (m1 + m2) / (rho * A),
# --- Pressure Logic
# ... pressure at bottom = pressure at the top + liquid head.
b.p ~ a.p + rho * g * h,
# --- Mass fraction
x1 ~ m1 / (m1 + m2),
x2 ~ m2 / (m1 + m2),
# --- Actual mass fractions at ports
x1a ~ ifelse(a.md > 0, instream(a.x1), x1),
x2a ~ ifelse(a.md > 0, instream(a.x2), x2),
x1b ~ ifelse(b.md > 0, instream(b.x1), x1),
x2b ~ ifelse(b.md > 0, instream(b.x2), x2),
# --- Setting stream values ---
a.x1 ~ x1,
a.x2 ~ x2,
b.x1 ~ x1,
b.x2 ~ x2,
# --- Dynamics ---
D(m1) ~ a.md * x1a + b.md * x1b,
D(m2) ~ a.md * x2a + b.md * x2b
]
compose(System(eqs, t, vars, pars; name = name), a, b)
end
Case 1: single tank
Component model:
# Parameters
# - initial masses in tanks
m0 = [2,1] #, [description = "Initial masses tank 1, kg"]
#
# - inlet flow rates and compositions
md_in = 2.0 #, [description = "Inlet mass flow tank 1, kg/s"]
x_in = [0.8,0.2] #, [description = "Inlet mass fractions tank 1, kg/kg"]
# - gas pressure boundary conditions
p_a = 101325.0 #, [description = "Atmospheric pressure
# - tank parameters
A = 0.3 #, [description = "Tank cross section area, m^2"]
rho = 1000.0 #, [description = "Fluid density, kg/m^3"]
# - valve flow coefficients
K = 0.3 #, [description = "Flow coefficient tank 1 to tank 2, kg/(s*Pa^0.5)"]
#
# Building the system
#
@named md_s = MassFlowPump(; md_set=md_in, x_set=x_in)
@named p_si = PressureReference(; p_set=p_a)
@named p_se = PressureBoundary(; p_set=p_a)
@named v_e = Valve(; Kv=K)
@named tnk = Tank(; Ap=A, rho_p=rho, m0=m0)
eqs_t1_2 = [
connect(md_s.a, p_si.a, tnk.a)
connect(tnk.b, v_e.a)
connect(v_e.b, p_se.a)
]
@named __sys_t1_2 = System(eqs_t1_2, t)
@named _sys_t1_2 = compose(__sys_t1_2, [md_s,p_si, p_se,v_e,tnk])
Compare this to Case 1 (blue ellipse). The pumped mass flow (md_s has a single port, a) is connected with the pressure reference (p_si; single port, a) at the top (port a) of the tank tnk.
At the bottom of the tank tnk (port b), the tank is connected to the inlet of the effluent valve (v_e). Finally, the outlet of the valve v_e is connected to the pressure boundary p_se (which may receive mass flow).
The symbolic model is simplified via mtkcompile, a numeric problem is generated via ODEProblem, and the solution is found:
sys_t1_2 = mtkcompile(_sys_t1_2)
tspan = (0.0, 10)
prob_t1_2 = ODEProblem(sys_t1_2, [], tspan)
sol_t1_2 = solve(prob_t1_2)
Note: the mtkcompile(_sys_t1_2) statement produce warnings:
Warning: Internal error: Variable tnk₊a₊md(t) was marked as being in 0 ~ -tnk₊x1a(t) + ifelse(tnk₊a₊md(t) > 0, ModelingToolkitBase.instream_rt(Val{2}(), Val{0}(), md_s₊a₊md(t), p_si₊a₊md(t), md_s₊a₊x1(t), p_si₊a₊x1(t)), tnk₊x1(t)), but was actually zero
└ @ ModelingToolkitTearing ...
I didn’t worry about these, but it would be nice if there was a way to get rid of them.
Case 2: two tanks
Component model (green circle):
# Parameters
# - initial masses in tanks
m0_1 = [2,1] #, [description = "Initial masses tank 1, kg"]
m0_2 = [5,8] #, [description = "Initial masses tank 2, kg"]
#
# - inlet flow rates and compositions
md_in_1 = 2.0 #, [description = "Inlet mass flow tank 1, kg/s"]
x_in_1 = [0.8,0.2] #, [description = "Inlet mass fractions tank 1, kg/kg"]
md_in_2 = 2.5 #, [description = "Inlet mass flow tank 2, kg/s"]
x_in_2 = [0.4,0.6] #, [description = "Inlet mass fractions tank 2, kg/kg"]
# - gas pressure boundary/reference conditions
p_a = 101325.0 #, [description = "Atmospheric pressure
# - tank parameters
A_1 = 0.3 #, [description = "Tank cross section area, m^2"]
A_2 = 0.6 #, [description = "Tank cross section area, m^2"]
rho = 1000.0 #, [description = "Fluid density, kg/m^3"]
# - valve flow coefficients
K_12 = 0.1 #, [description = "Flow coefficient tank 1 to tank 2, kg/(s*Pa^0.5)"]
K_e = 0.4 #, [description = "Flow coefficient tank 2 to tank 3, kg/(s*Pa^0.5)"]
#
@named m_s1 = MassFlowPump(; md_set=md_in_1, x_set=x_in_1)
@named m_s2 = MassFlowPump(; md_set=md_in_2, x_set=x_in_2)
@named p_s1 = PressureReference(; p_set=p_a)
@named p_s2 = PressureReference(; p_set=p_a)
@named p_e = PressureBoundary(; p_set=p_a)
@named v_12 = Valve(; Kv=K_12)
@named v_e = Valve(; Kv=K_e)
@named t_1 = Tank(; Ap=A_1, rho_p=rho, m0=m0_1)
@named t_2 = Tank(; Ap=A_2, rho_p=rho, m0=m0_2)
eqs_t2 = [
connect(m_s1.a, p_s1.a, t_1.a)
connect(m_s2.a, p_s2.a, t_2.a)
connect(t_1.b, v_12.a)
connect(v_12.b, t_2.b,v_e.a)
connect(v_e.b, p_e.a)
]
@named __sys_t2 = System(eqs_t2, t)
@named _sys_t2 = compose(__sys_t2, [m_s1, m_s2, p_s1, p_s2, p_e,
v_12, v_e, t_1, t_2])
sys_t2 = mtkcompile(_sys_t2)
tspan = (0.0, 10)
prob_t2 = ODEProblem(sys_t2, [], tspan)
sol_t2 = solve(prob_t2)
Again, same warning after mtkcompile. But the solution is found, and makes perfect sense:
Because of a change in which tank has the highest level (highest bottom pressure), the flow direction \dot{m}_{1,2} changes direction, with consequences for which mass fractions are transferred:
Case 3: three tanks
Component model:
# Parameters
# - initial masses in tanks
m0_1 = [2,1] #, [description = "Initial masses tank 1, kg"]
m0_2 = [5,8] #, [description = "Initial masses tank 2, kg"]
m0_3 = [3,5] #, [description = "Initial masses tank 3, kg"]
#
# - inlet flow rates and compositions
md_in_1 = 2.0 #, [description = "Inlet mass flow tank 1, kg/s"]
x_in_1 = [0.8,0.2] #, [description = "Inlet mass fractions tank 1, kg/kg"]
md_in_2 = 2.5 #, [description = "Inlet mass flow tank 2, kg/s"]
x_in_2 = [0.4,0.6] #, [description = "Inlet mass fractions tank 2, kg/kg"]
# - gas pressure boundary conditions
p_a = 101325.0 #, [description = "Atmospheric pressure
# - tank parameters
A_1 = 0.3 #, [description = "Tank cross section area, m^2"]
A_2 = 0.6 #, [description = "Tank cross section area, m^2"]
A_3 = 0.4 #, [description = "Tank cross section area, m^2"]
rho = 1000.0 #, [description = "Fluid density, kg/m^3"]
# - valve flow coefficients
K_12 = 0.1 #, [description = "Flow coefficient tank 1 to tank 2, kg/(s*Pa^0.5)"]
K_23 = 0.4 #, [description = "Flow coefficient tank 2 to tank 3, kg/(s*Pa^0.5)"]
K_3e = 0.3 #, [description = "Flow coefficient effluent tank 3, kg/(s*Pa^0.5)"]
#
@named mp_s1 = MassFlowPump(; md_set=md_in_1, x_set=x_in_1)
@named mp_s2 = MassFlowPump(; md_set=md_in_2, x_set=x_in_2)
@named p_r1 = PressureReference(; p_set=p_a)
@named p_r2 = PressureReference(; p_set=p_a)
@named p_r3 = PressureReference(; p_set=p_a)
@named p_be = PressureBoundary(; p_set=p_a)
@named va_12 = Valve(; Kv=K_12)
@named va_23 = Valve(; Kv=K_23)
@named va_3e = Valve(; Kv=K_3e)
@named tnk_1 = Tank(; Ap=A_1, rho_p=rho, m0=m0_1)
@named tnk_2 = Tank(; Ap=A_2, rho_p=rho, m0=m0_2)
@named tnk_3 = Tank(; Ap=A_3, rho_p=rho, m0=m0_3)
eqs_t3 = [
connect(mp_s1.a, p_r1.a, tnk_1.a)
connect(mp_s2.a, p_r2.a, tnk_2.a)
connect(tnk_1.b, va_12.a)
connect(va_12.b, tnk_2.b, va_23.a)
connect(va_23.b, tnk_3.a, p_r3.a)
connect(tnk_2.b, va_3e.a)
connect(va_3e.b, p_be.a)
]
@named __sys_t3 = System(eqs_t3, t)
@named _sys_t3 = compose(__sys_t3, [mp_s1, mp_s2, p_r1, p_r2, p_r3, p_be,
va_12, va_23, va_3e,
tnk_1, tnk_2, tnk_3])
sys_t3 = mtkcompile(_sys_t3)
Same warning as before, but the simplified symbolic model is created. Next, creating numeric problem:
tspan = (0.0, 10)
prob_t3 = ODEProblem(sys_t3, [], tspan)
This works, and creates the expected output:
ODEProblem with uType Vector{Float64} and tType Float64. In-place: true
Initialization status: FULLY_DETERMINED
Non-trivial mass matrix: false
timespan: (0.0, 10.0)
u0: 6-element Vector{Float64}:
5.0
3.0
8.0
5.0
1.0
2.0
Finally, attempting to solve the problem:
sol_t3 = solve(prob_t3)
crashes with the following error message:
MethodError: no method matching instream(::Float64)
The function `instream` exists, but no method is defined for this combination of argument types.
Closest candidates are:
instream(!Matched::SymbolicUtils.BasicSymbolicImpl.var"typeof(BasicSymbolicImpl)"{SymReal})
@ ModelingToolkitBase C:\Users\Be...