I’m in the process of building a set of modeling blocks for a gas system, and I’m currently on working on orifices. Despite having registered my functions and symbolics, I get MethodError for what I can only assume to be the functions (the error is a little cryptic).
ERROR: MethodError: no method matching AbstractFloat(::Type{SymbolicUtils.BasicSymbolic{Real}})
A MEW can be found below.
using ModelingToolkit, Plots, DifferentialEquations, NonlinearSolve
using ModelingToolkit: t_nounits as t, D_nounits as D
using Symbolics
using CoolProp
PropsSI(out::AbstractString, name1::AbstractString, value1::Real, name2::AbstractString, value2::Real, fluid::AbstractString) = CoolProp.PropsSI(out, name1, value1, name2, value2, fluid)
@register_symbolic PropsSI(out::AbstractString, name1::AbstractString, value1::Real, name2::AbstractString, value2::Real, fluid::AbstractString)
function mass_flow_gas_orifice(d::Real,P₁::Real,h₁::Real,P₂::Real,h₂::Real,medium::AbstractString)
if P₁ > P₂ # Set flow direction
P_us = P₁ # Upstream
h_us = h₁ # Upstream
P_ds = P₂ # Downstream
h_ds = h₂ # Downstream
dir = 1 # For the sign on the mass flow
else # And vice versa
P_ds = P₁
h_ds = h₁
P_us = P₂
h_us = h₂
dir = -1
end
A = (d/2)^2 * π
γ = PropsSI("ISENTROPIC_EXPANSION_COEFFICIENT","P",P_us,"H",h_us,medium)
ρ_us = PropsSI("D","P",P_us,"H",h_us,medium)
if P_ds/P_us < (2/(γ+1))^(γ/(γ-1)) # Is choked
qₘ = A*√(γ*ρ_us*P_us*(2/(γ+1))^((γ+1)/(γ-1)))
else # Is not choked
qₘ = A*√(2*ρ_us*P_us*γ/(γ-1)*((P_ds/P_us)^(2/γ)-(P_ds/P_us)^((γ+1)/γ)))
end
return dir*qₘ
end
@register_symbolic mass_flow_gas_orifice(d::Real,P₁::Real,h₁::Real,P₂::Real,h₂::Real,medium::AbstractString)
@connector FluidPort begin
P(t)
h_outflow(t), [connect = Stream]
dm(t), [connect = Flow]
end
@mtkmodel P_boundary begin
@parameters begin
P = 1
end
@variables begin
dm(t), [guess = 1e-2]
h(t)
end
@components begin
port = FluidPort()
end
@equations begin
port.P ~ P
instream(port.h_outflow) ~ h
0 ~ port.dm + dm
end
end
@mtkmodel Ph_boundary begin
@parameters begin
P = 1
h = 1
end
@variables begin
dm(t), [guess = 1e-2]
end
@components begin
port = FluidPort()
end
@equations begin
port.P ~ P
instream(port.h_outflow) ~ h
0 ~ port.dm + dm
end
end
@mtkmodel FluidTwoPort begin
@variables begin
dP(t)
h(t)
dm(t)
end
@components begin
port_a = FluidPort()
port_b = FluidPort()
end
@equations begin
dP ~ port_a.P - port_b.P
dm ~ port_a.dm
port_b.h_outflow ~ instream(port_b.h_outflow)
port_a.h_outflow ~ instream(port_a.h_outflow)
port_a.h_outflow ~ port_b.h_outflow
0 ~ port_a.dm + port_b.dm
end
end
@mtkmodel GasOrifice begin
@extend FluidTwoPort()
@parameters begin
d
end
@equations begin
dm ~ mass_flow_gas_orifice(d,port_a.P,port_a.h_outflow,port_b.P,port_b.h_outflow,"AIR")
end
end
@mtkmodel FlowResistor begin
@extend FluidTwoPort()
@parameters begin
R = 1
end
@equations begin
dm ~ √(dP) / R
end
end
@mtkmodel System begin
@parameters begin
P1 = 1
P2 = 1
h1 = 1
#h2 = 1
d1 = 1
d2 = 1
end
@components begin
inlet = Ph_boundary(;P=P1,h=h1)
outlet = P_boundary(;P=P2)
orifice1 = GasOrifice(;d=d1)
orifice2 = GasOrifice(;d=d2)
end
@equations begin
connect(inlet.port, orifice1.port_a)
connect(orifice1.port_b,orifice2.port_a)
connect(orifice2.port_b,outlet.port)
end
end
hin = PropsSI("H","P",5e5,"T",293,"AIR")
@mtkbuild sys = System(;P1=5e5, h1=hin, P2=1e5, d1=2e-3, d2=1e-3)
prob = ODEProblem(sys, [], (0,1))
sol = solve(prob)
Any guidance is much appreciated. I’ve been prodding around for a few days on this on/off, and can’t seem to figure it out. My intuition was to change the types of the register_symbolic functions, but this didn’t do anything.
For what it’s worth, I’m only interested in the steady state solution for now.
Cheers