Hello everyone,
I’m currently working on a dynamic of microgas turbine system model in Julia using a system of Differential-Algebraic Equations (DAE).
I want to adapt modelica acausal modelling problem into Julia using modellingtoolkit.jl and i already create a code for the system that look like this :
using ModelingToolkit, DifferentialEquations, CoolProp, Plots
@independent_variables t
D = Differential(t)
const fluid = "Air"
using ModelingToolkit, DifferentialEquations, CoolProp, Plots
@independent_variables t
D = Differential(t)
const fluid = "Air"
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 h(p_in, T_in)
h = PropsSI("H", "P", p_in, "T", T_in, fluid)
return PropsSI("H", "P", p_in, "T", T_in, fluid)
end
function hf(p_in, T_in)
hf = PropsSI("H", "P", p_in, "T", T_in, fluid)
return PropsSI("H", "P", p_in, "T", T_in, fluid)
end
function ρ(p_in, T_in)
ρ = PropsSI("D", "P", p_in, "T", T_in, fluid)
return PropsSI("D", "P", p_in, "T", T_in, fluid)
end
function u(p_in, T_in)
u = PropsSI("D", "P", p_in, "T", T_in, fluid)
return PropsSI("D", "P", p_in, "T", T_in, fluid)
end
function cv(T)
cv = PropsSI("Cvmass", "T", T, "P", 101325, fluid)
return PropsSI("Cvmass", "T", T, "P", 101325, fluid)
end
function cp(T)
cp = PropsSI("Cpmass", "T", T, "P", 101325, fluid)
return PropsSI("Cpmass", "T", T, "P", 101325, fluid)
end
@register_symbolic h(p_in, T_in)
@register_symbolic hf(p_in, T_in)
@register_symbolic ρ(p_in, T_in)
@register_symbolic u(p_in, T_in)
@register_symbolic cv(T)
@register_symbolic cp(T)
@connector function Port(;name)
vars = @variables begin
p(t), [input = true]
T(t), [input = true]
P(t)
mdot(t), [connect = Flow]
end
ODESystem(Equation[], t, vars, [];name=name)
end
function Source(;name,source_temperature, source_pressure,source_mdot, source_power)
@named port = Port()
para = @parameters begin
end
vars = @variables begin
end
eqs = [
port.mdot ~ source_mdot
port.P ~ source_power
port.p ~ source_pressure
port.T ~ source_temperature
]
compose(ODESystem(eqs, t, vars, para;name),port)
end
function Compressor(;name)
@named inport = Port()
@named outport = Port()
para = @parameters begin
p1 = 101020, [description = "Ambient Pressure (Pa)"]
T1 = 300, [description = "Ambient Temperature (K)"]
m_cas_comp = 105, [description = "Mass Casing Compressor (kg)"]
c_casing = 500, [description = "Specific Heat Casing Compressor (J/kg·K)"]
ṁ_comp = 0.4, [description = "Mass Flow Rate Compressor (kg/s)"]
η_c_is = 0.85, [description = "Isentropic Efficiency"]
π_c = 2.0, [description = "Compressor Pressure Ratio"]
κa = cp(inport.T)/cv(inport.T)
γ_a = (κa-1) / κa
end
vars = @variables begin
T2(t), [description = "Temperature outlet compressor (K)"]
end
eqs = [ D(T2) ~ (1 / (m_cas_comp * c_casing)) * (ṁ_comp *(h(outport.p,outport.T) - h(outport.p, T2)))
outport.p ~ π_c * inport.p
outport.T ~ inport.T * (1 + (1 / η_c_is) * (π_c^(γ_a) - 1))
outport.mdot ~ inport.mdot
outport.P ~ inport.P - (ṁ_comp * (h(outport.p, outport.T) - h(inport.p, inport.T)))]
compose(ODESystem(eqs, t, vars, para;name), inport, outport)
end
function Recuperator(; name, As = 3)
@named inport_hot = Port()
@named inport_cold = Port()
@named outport_hot = Port()
@named outport_cold = Port()
para = @parameters begin
A = 10.0 / As, [description = "Heat transfer area per section (m²)"]
cp_h = 1000.0, [description = "Hot fluid heat capacity (J/(kg·K)"]
cp_c = 1000.0, [description = "Cold fluid heat capacity (J/(kg·K)"]
m_dot_h = 0.65, [description = "Hot side mass flow rate (kg/s)"]
m_dot_c = 0.65, [description = "Cold side mass flow rate (kg/s)"]
h_h = 100.0, [description = "Hot side heat transfer coefficient (W/m²K)"]
h_c = 100.0, [description = "Cold side heat transfer coefficient (W/m²K)"]
m_w = 185.0, [description = "Wall mass per section (kg)"]
c_w = 500.0, [description = "Wall heat capacity (J/(kg·K)"]
σ_h = 0.98, [description = "Hot side pressure loss coefficient"]
σ_c = 0.95, [description = "Cold side pressure loss coefficient"]
end
vars = @variables begin
T_h_out(t)[1:As], [description = "Hot fluid outlet temperatures (K)"]
T_c_out(t)[1:As], [description = "Cold fluid outlet temperatures (K)"]
T_w(t)[1:As], [description = "Wall temperatures (K)"]
end
eqs = []
# Temperature equations for each section
for i in 1:As
if i == 1
T_h_in_i = inport_hot.T
else
T_h_in_i = T_h_out[i-1]
end
if i == As
T_c_in_i = inport_cold.T
else
T_c_in_i = T_c_out[i+1]
end
# Hot fluid energy balance
push!(eqs, T_h_out[i] ~ (m_dot_h * cp_h * T_h_in_i - 0.5*h_h*A*T_h_in_i + h_h*A*T_w[i]) / (m_dot_h * cp_h + 0.5*h_h*A))
# Cold fluid energy balance (counter-flow)
push!(eqs, T_c_out[i] ~ (m_dot_c * cp_c * T_c_in_i - 0.5*h_c*A*T_c_in_i + h_c*A*T_w[i]) / (m_dot_c * cp_c + 0.5*h_c*A))
# Wall energy balance (dynamic equation)
Q_h = m_dot_h * cp_h * (T_h_in_i - T_h_out[i])
Q_c = m_dot_c * cp_c * (T_c_out[i] - T_c_in_i)
push!(eqs, D(T_w[i]) ~ (Q_h - Q_c) / (m_w * c_w))
end
# System boundary conditions
push!(eqs,
outport_hot.T ~ T_h_out[end], # Connect last hot section to outlet
outport_cold.T ~ T_c_out[1], # Connect first cold section to outlet
outport_hot.p ~ inport_hot.p * σ_h, # Pressure loss
outport_cold.p ~ inport_cold.p * σ_c,
outport_hot.mdot ~ -m_dot_h, # Mass flow convention
outport_cold.mdot ~ -m_dot_c,
inport_hot.mdot ~ m_dot_h,
inport_cold.mdot ~ m_dot_c)
compose(ODESystem(eqs, t, vars, para;name), [inport_hot, inport_cold], [outport_hot, outport_cold])
end
function CombustionChamber(;name)
@named inport = Port()
@named outport = Port()
pars = @parameters begin
η_comb = 0.97, [description = "Combustor efficiency"]
σ_comb = 0.945, [description = "Artificial Pressure Losses"]
V_comb = 5.5422e-3, [description = "Combustor Volume (m^3)"]
LHV = 40.564e6, [description = "Lower Heating Value (J/kg)"]
m_cas_comb = 10.0, [description = "Combustor casing mass (kg)"]
c_casing = 0, [description = "Specific heat capacity of casing (J/kg·K)"]
Tf = 303.15, [description = "Fuel Temperature (K)"]
pf = 600000, [description = "Fuel Pressure (Pa) (6 bar)"]
ρf = 0.8217, [description = "Fuel Density (kg/m^3)"]
ρ3 = 1.15, [description = "Density in combustor (kg/m^3)"]
R = 287, [description = "Specific gas constant for air (J/kg·K)"]
Rmol = 8.314
n = 3
ṁ_comp = 0.4
ṁ_turb = 0.46
ṁ_fuel = 0.06
end
vars = @variables begin
ρ3(t), [description = "Density in combustor"]
u3(t), [description = "Internal energy in combustor"]
end
# Define equations (steady-state: time derivatives are zero)
eqs = [ D(ρ3) ~ (ṁ_comp + ṁ_fuel - ṁ_turb) / V_comb,
D(u3) ~ (1 / (ρ3 * V_comb)) * (ṁ_comp * h(inport.p, inport.T) + ṁ_fuel * hf(pf, Tf) + (η_comb * ṁ_fuel * LHV) - ṁ_turb * h(outport.p, outport.T) - u3*(ṁ_comp + ṁ_fuel - ṁ_turb)),
u3 ~ u(outport.p, outport.T),
outport.p ~ ρ3 * R * outport.T,
outport.mdot ~ inport.mdot,
outport.P ~ 0]
compose(ODESystem(eqs, t, name=name), inport, outport)
end
function Turbine(;name)
@named inport = Port()
@named outport = Port()
para = @parameters begin
m_cas_turb = 105, [description = "Turbine casing mass (kg)"]
c_casing = 600, [description = "Specific Heat Casing Turbine (J/kg·K)"]
ṁ_t = 0.46, [description = "Mass Flow Turbine (kg/s)"]
η_turb_is = 0.95, [description = "Turbine Isentropic Efficiency"]
π_t = 2.0, [description = "Turbine Pressure Ratio"]
κg = cp(inport.T)/cv(inport.T)
γ_g = (κg - 1) / κg
end
vars = @variables begin
T4(t), [description = "Turbine Outlet Temperature (K)"]
end
eqs = [ D(T4) ~ (1 / (m_cas_turb * c_casing)) * (ṁ_t * (h(outport.p, outport.T) - h(outport.p, T4)))
outport.p ~ π_t * inport.p
outport.T ~ inport.T * (1 - η_turb_is * (1 - π_t^(-γ_g)))
outport.mdot ~ inport.mdot
outport.P ~ inport.P + (ṁ_t * (h(inport.p, inport.T) - h(outport.p, outport.T)))]
compose(ODESystem(eqs, t, vars, para;name), inport, outport)
end
function RotationalSystem(;name)
@named inport = Port()
@named outport = Port()
para = @parameters begin
I = 8.3e-3, [description = "Inertia (kg.m^2)"]
P_lb = 1500, [description = "Power Bearing Losses (W)"]
P_le = 10000, [description = "Power Electrical Losses (W)"]
Pload = 70, [description = "Power Load (W)"]
end
vars = @variables begin
end
eqs = [ outport.p ~ inport.p
outport.T ~ inport.T
outport.mdot ~ inport.mdot
outport.P ~ inport.P - P_lb - P_le]
compose(ODESystem(eqs, t, vars, para;name), inport, outport)
end
function Sink(;name)
@named port = Port()
para = @parameters begin
end
vars = @variables begin
p(t)
T(t)
mdot(t)
P(t)
end
eqs = [
mdot ~ port.mdot
p ~ port.p
T ~ port.T
P ~ port.P
]
compose(ODESystem(eqs, t, vars, para;name),port)
end
@named source = Source(source_temperature = start_T, source_pressure = start_p,source_mdot = mdot, source_power = start_P)
@named compressor = Compressor()
@named recuperator = Recuperator()
@named combustionchamber = CombustionChamber()
@named turbine = Turbine()
@named sink = Sink()
eqs = [
connect(source.port,compressor.inport)
connect(compressor.outport,turbine.inport)
connect(turbine.outport,sink.port)
]
@named model = ODESystem(eqs, t, systems=[source,compressor,turbine,sink])
However, I’m facing an issue:
- In CombustionChamber model, Extracting Output Temperature (output.T):
I need to extract the output temperature T from the combustion chamber using CoolProp, but the problem is that my DAE system requires this T as an input for another equation. This creates a dependency issue, and I’m not sure how to properly structure my equations to resolve it. - In Recuperator model, it causes this error message :
MethodError: no method matching ModelingToolkit.ODESystem(::Vector{Any}, ::Symbolics.Num, ::Vector{Symbolics.Arr{Symbolics.Num, 1}}, ::Vector{Symbolics.Num}; name::Symbol)
The type `ModelingToolkit.ODESystem` exists, but no method is defined for this combination of argument types when trying to construct it.
Closest candidates are:
ModelingToolkit.ODESystem(::Any, ::Any, ::Any, ::Any, ::Any, ::Any, ::Any, ::Any, ::Any, ::Any, ::Any, ::Any, ::Any, ::Any, ::Any, ::Any, ::Any, ::Any, ::Any, ::Any, ::Any, ::Any, ::Any, ::Any, ::Any, ::Any, ::Any, ::Any; ...)
@ ModelingToolkit C:\Users\Andhika\.julia\packages\ModelingToolkit\eiNg3\src\systems\diffeqs\odesystem.jl:184
ModelingToolkit.ODESystem(::Any, ::Any, ::Any, ::Any, ::Any, ::Any, ::Any, ::Any, ::Any, ::Any, ::Any, ::Any, ::Any, ::Any, ::Any, ::Any, ::Any, ::Any, ::Any, ::Any, ::Any, ::Any, ::Any, ::Any, ::Any, ::Any, ::Any, ::Any, ::Any; ...)
@ ModelingToolkit C:\Users\Andhika\.julia\packages\ModelingToolkit\eiNg3\src\systems\diffeqs\odesystem.jl:184
ModelingToolkit.ODESystem(::Any, ::Any, ::Any, ::Any, ::Any, ::Any, ::Any, ::Any, ::Any, ::Any, ::Any, ::Any, ::Any, ::Any, ::Any, ::Any, ::Any, ::Any, ::Any, ::Any, ::Any, ::Any, ::Any, ::Any, ::Any, ::Any, ::Any, ::Any, ::Any, ::Any; ...)"
I would really appreciate any guidance on how to handle these issues. If needed, I can provide more details or share parts of my code. Thank you in advance!