I’m working on converting an existing Julia model using DifferentialEquations.jl and solving Stiff ODEs to using MTK. There seems to be some nuance to setting up ODESystems that I don’t understand. I’m modeling pressure drop in a pipe either implicitly, explicitly, or using two separate equations (sqrt of pressure drop then squaring back again). I don’t understand why the implicit and two equation explicit model works, while the explicit model with one equation doesn’t. I’m only a week new to MTK and I’m mostly familiar with the acausal RC circuit example. If there’s other suggestions you have about the model, I’d love to hear them! Thank you!
#=
There are three versions of the PipeResistance function:
PipeResistance_implicit
PipeResistance_explicit
PipeResistance_explict_sqrt
each of these deal with the Δp differently.
PipeResistance_implicit -> does not explicity deal with pressure_drop
PipeResistance_explicit -> added an equation that causes MAXITERS - we do not understand why
Δp ~ pressure_drop_pipe(mdot_ab, L, 𝑓, d_h, ρ, A_flow)
pressure_drop_pipe(ṁ, L, 𝑓, dₕ, ρ, Aᵪ) = ṁ * abs(ṁ) * L * 𝑓 / (2 * dₕ * ρ * Aᵪ^2)
PipeResistance_explict_sqrt -> attempts to fix this by adding yet another equation
sqrt_Δp ~ pressure_drop_pipe_squareroot(mdot_ab, L, 𝑓, d_h, ρ, A_flow)
Δp ~ sqrt_Δp * abs(sqrt_Δp)
pressure_drop_pipe_squareroot(ṁ, L, 𝑓, dₕ, ρ, Aᵪ) = ṁ * sqrt(L * 𝑓 / (2 * dₕ * ρ * Aᵪ^2))
=#
Full code:
using ModelingToolkit, Plots, DifferentialEquations, BenchmarkTools
gr()
ṁ_pipe_friction(Δp, 𝑓, ρ, L, dₕ, Aᵪ) = sign(Δp) * Aᵪ * sqrt(2 * dₕ * ρ / (L * 𝑓)) * sqrt(abs(Δp)) # mass flow based on frictional pressure drop
pressure_drop_pipe(ṁ, L, 𝑓, dₕ, ρ, Aᵪ) = ṁ * abs(ṁ) * L * 𝑓 / (2 * dₕ * ρ * Aᵪ^2) # frictional pressure drop
pressure_drop_pipe_squareroot(ṁ, L, 𝑓, dₕ, ρ, Aᵪ) = ṁ * sqrt(L * 𝑓 / (2 * dₕ * ρ * Aᵪ^2)) # frictional pressure drop
function Port(;name)
@variables p(t) mdot(t) ρ(t)
ODESystem(Equation[], t, [p, mdot, ρ], [], name=name, defaults=[p=>1e5, mdot=>0.0])
end
function PipeResistance_implicit(;name, L = 1.0, A_flow = 0.01, d_h = 1.0)
val_L = L
val_A_flow = A_flow
val_d_h = d_h
@named a = Port()
@named b = Port()
@variables Δp(t) 𝑓(t) ρ(t) mdot_ab(t)
@parameters L A_flow d_h
eqs = [
Δp ~ a.p - b.p
0 ~ a.mdot + b.mdot
ρ ~ (a.ρ + b.ρ) / 2
𝑓 ~ 0.03
mdot_ab ~ a.mdot
mdot_ab ~ ṁ_pipe_friction(Δp, 𝑓, ρ, L, d_h, A_flow)
]
ODESystem(eqs, t, [Δp, 𝑓, ρ, mdot_ab], [L, A_flow, d_h], systems=[a, b], defaults=Dict(L => val_L, A_flow => val_A_flow, d_h => val_d_h), name=name)
end
function PipeResistance_explicit(;name, L = 1.0, A_flow = 0.01, d_h = 1.0)
val_L = L
val_A_flow = A_flow
val_d_h = d_h
@named a = Port()
@named b = Port()
@variables Δp(t) 𝑓(t) ρ(t) mdot_ab(t)
@parameters L A_flow d_h
eqs = [
Δp ~ a.p - b.p
0 ~ a.mdot + b.mdot
ρ ~ (a.ρ + b.ρ) / 2
𝑓 ~ 0.03
mdot_ab ~ a.mdot
Δp ~ pressure_drop_pipe(mdot_ab, L, 𝑓, d_h, ρ, A_flow)
]
ODESystem(eqs, t, [Δp, 𝑓, ρ, mdot_ab], [L, A_flow, d_h], systems=[a, b], defaults=Dict(L => val_L, A_flow => val_A_flow, d_h => val_d_h), name=name)
end
function PipeResistance_explicit_sqrt(;name, L = 1.0, A_flow = 0.01, d_h = 1.0)
val_L = L
val_A_flow = A_flow
val_d_h = d_h
@named a = Port()
@named b = Port()
@variables Δp(t) 𝑓(t) ρ(t) mdot_ab(t) sqrt_Δp(t)
@parameters L A_flow d_h
eqs = [
Δp ~ a.p - b.p
0 ~ a.mdot + b.mdot
ρ ~ (a.ρ + b.ρ) / 2
𝑓 ~ 0.03
mdot_ab ~ a.mdot
sqrt_Δp ~ pressure_drop_pipe_squareroot(mdot_ab, L, 𝑓, d_h, ρ, A_flow)
Δp ~ sqrt_Δp * abs(sqrt_Δp)
]
ODESystem(eqs, t, [Δp, 𝑓, ρ, mdot_ab, sqrt_Δp], [L, A_flow, d_h], systems=[a, b], defaults=Dict(L => val_L, A_flow => val_A_flow, d_h => val_d_h), name=name)
end
function FluidVolume(; name, V = 1.0)
val_V = V
@named a = Port()
@named b = Port()
@variables m(t) p(t) ρ(t)
@parameters V
D = Differential(t)
eqs = [
D(m) ~ a.mdot + b.mdot
ρ ~ m / V
ρ ~ a.ρ # transfer fluid state to port
ρ ~ b.ρ # transfer fluid state to port
p ~ 1e5 # need to add urho_lookup
p ~ a.p # transfer fluid state to port
p ~ b.p # transfer fluid state to port
]
ODESystem(eqs, t, [m, p, ρ], [V], systems=[a, b], defaults=Dict(V=> val_V), name=name)
end
function ConstantReservoir(;name, p, ρ)
val_p = p
val_ρ = ρ
@named a = Port()
@parameters p ρ
eqs = [
a.p ~ p
a.ρ ~ ρ
]
ODESystem(eqs, t, [], [p, ρ], systems=[a], defaults=Dict(p => val_p, p => val_p, ρ => val_ρ), name=name)
end
function connect_ports(ps...)
eqs = [
0 ~ sum(p->p.mdot, ps) # sum(mdot) = 0
]
# Transfer node data to edge data
for i in 2:length(ps)
push!(eqs, ps[1].p ~ ps[i].p)
push!(eqs, ps[1].ρ ~ ps[i].ρ)
end
return eqs
end
V = 1.0
P_low = 1e5
P_high = 1.1e5
ρ = 10
L = 1.0
d_h = 0.0254
A_flow = pi*d_h^2/4
@parameters t
@register ṁ_pipe_friction(Δp, 𝑓, ρ, L, dₕ, Aᵪ)
@register pressure_drop_pipe(ṁ, L, 𝑓, dₕ, ρ, Aᵪ)
@register pressure_drop_pipe_squareroot(ṁ, L, 𝑓, dₕ, ρ, Aᵪ)
function solve_implicit(P_high, ρ, L, A_flow, d_h, P_low, V)
@named source = ConstantReservoir(p=P_high, ρ=ρ)
@named valve = PipeResistance_implicit(L=L, A_flow=A_flow, d_h=d_h)
@named f_volume = FluidVolume(V=V)
@named valve2 = PipeResistance_implicit(L=L, A_flow=A_flow, d_h=d_h)
@named sink = ConstantReservoir(p=P_low, ρ=ρ)
cnnt_eqs = [
connect_ports(source.a, valve.a)
connect_ports(valve.b, f_volume.a)
connect_ports(f_volume.b, valve2.a)
connect_ports(valve2.b, sink.a)
]
@named tes_model = ODESystem(cnnt_eqs, t,
systems=[source, valve, f_volume, valve2, sink])
sys = structural_simplify(tes_model)
u0 = [
f_volume.m => 1.0
f_volume.p => 1e5
]
prob = ODAEProblem(sys, u0, (0, 10.0))
sol = solve(prob, Tsit5())
#plot(sol)
return sol
end
function solve_explicit_sqrt(P_high, ρ, L, A_flow, d_h, P_low, V)
@named source = ConstantReservoir(p=P_high, ρ=ρ)
@named valve = PipeResistance_explicit_sqrt(L=L, A_flow=A_flow, d_h=d_h)
@named f_volume = FluidVolume(V=V)
@named valve2 = PipeResistance_explicit_sqrt(L=L, A_flow=A_flow, d_h=d_h)
@named sink = ConstantReservoir(p=P_low, ρ=ρ)
cnnt_eqs = [
connect_ports(source.a, valve.a)
connect_ports(valve.b, f_volume.a)
connect_ports(f_volume.b, valve2.a)
connect_ports(valve2.b, sink.a)
]
@named tes_model = ODESystem(cnnt_eqs, t,
systems=[source, valve, f_volume, valve2, sink])
sys = structural_simplify(tes_model)
u0 = [
f_volume.m => 1.0
f_volume.p => 1e5
]
prob = ODAEProblem(sys, u0, (0, 10.0))
sol = solve(prob, Tsit5())
#plot(sol)
return sol
end
function solve_explicit(P_high, ρ, L, A_flow, d_h, P_low, V)
@named source = ConstantReservoir(p=P_high, ρ=ρ)
@named valve = PipeResistance_explicit(L=L, A_flow=A_flow, d_h=d_h)
@named f_volume = FluidVolume(V=V)
@named valve2 = PipeResistance_explicit(L=L, A_flow=A_flow, d_h=d_h)
@named sink = ConstantReservoir(p=P_low, ρ=ρ)
cnnt_eqs = [
connect_ports(source.a, valve.a)
connect_ports(valve.b, f_volume.a)
connect_ports(f_volume.b, valve2.a)
connect_ports(valve2.b, sink.a)
]
@named tes_model = ODESystem(cnnt_eqs, t,
systems=[source, valve, f_volume, valve2, sink])
sys = structural_simplify(tes_model)
u0 = [
f_volume.m => 1.0
f_volume.p => 1e5
]
prob = ODAEProblem(sys, u0, (0, 10.0))
sol = solve(prob, Tsit5())
#plot(sol)
return sol
end
@btime solve_implicit(P_high, ρ, L, A_flow, d_h, P_low, V)
@btime solve_explicit_sqrt(P_high, ρ, L, A_flow, d_h, P_low, V)
@btime solve_explicit(P_high, ρ, L, A_flow, d_h, P_low, V)