 # Modeling Thermodynamic System in MTK

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(ṁ, L, 𝑓, dₕ, ρ, Aᵪ) = ṁ * abs(ṁ) * 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(ṁ, L, 𝑓, dₕ, ρ, Aᵪ) = ṁ * sqrt(L * 𝑓 / (2 * dₕ * ρ * Aᵪ^2))

=#
``````

Full code:

``````using ModelingToolkit, Plots, DifferentialEquations, BenchmarkTools
gr()

ṁ_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(ṁ, L, 𝑓, dₕ, ρ, Aᵪ) = ṁ * abs(ṁ) * L * 𝑓 / (2 * dₕ * ρ * Aᵪ^2) # frictional pressure drop
pressure_drop_pipe_squareroot(ṁ, L, 𝑓, dₕ, ρ, Aᵪ) = ṁ * 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 ~ ṁ_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.p ~ ps[i].p)
push!(eqs, ps.ρ ~ 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 ṁ_pipe_friction(Δp, 𝑓, ρ, L, dₕ, Aᵪ)
@register pressure_drop_pipe(ṁ, L, 𝑓, dₕ, ρ, Aᵪ)
@register pressure_drop_pipe_squareroot(ṁ, 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)
``````

The failure is due to the inappropriate initial condition `mdot=0.0`. The system is initialized at the exact point where `abs(ṁ)`'s derivative is not defined. Writing

``````ODESystem(Equation[], t, [p, mdot, ρ], [], name=name, defaults=[p=>1e5, mdot => -0.1])
``````

solves the problem.

2 Likes