How to connect MTK thermal Heatport?

Zhijian · February 22, 2024, 3:22pm

Hi guys,

I am new to Julia and am using ModelingToolkit to try to build a simple heat conduction model (3R2C wall model).

So it will look like this…

I am trying to create a constant ambient temperature then switch that to a hourly weather data. However, I cannot even finish my first step, that is connect the left-resistor to the ambient.

In my code, I am using a heat port, and this is how I try to connect it, which didn’t report any fault but just do nothing to the whole system

@named Resistor_conv_o = ThermalResistor(R = R_conv_o)
@named T_weather = HeatPort(T = 273.15 + 20)
connect(Resistor_conv_o.port_a, T_weather)

I would really appreciate someone’s help!

After figuring this out, I also would like to ask for some insight of how can I switch the constant ambient temeprature data linked with a series of weather data (file/array…). I used to have something like this in my own model:

hour_index = Int(floor(t) % 24 + 1) # in hourly resolution
    T_out = T_out_hourly[hour_index]

Many thanks!

See the full code of the 3R2C model as below:

using ModelingToolkit
using ModelingToolkitStandardLibrary.Thermal, ModelingToolkit, OrdinaryDiffEq, Plots

@parameters t

A = 20.0 #wall area
h = 1.0 #Convective coefficient
C_wall_out = 20000
C_wall_in = 20000
C_zone = 70000

R_conv_o = 1/(R*A)
R_conv_i = 1/(R*A)
R_wall = 0.13
#TempConstant = 20

# Capacitors
@named Capacitor_wall_out = HeatCapacitor(C = C_wall_out, T = 273.15 + 20)
@named Capacitor_wall_in = HeatCapacitor(C = C_wall_in, T = 273.15 + 50)
@named Capacitor_zone = HeatCapacitor(C = C_zone, T = 273.15 + 50)

# Resistors
@named Resistor_conv_o = ThermalResistor(R = R_conv_o)
@named Resistor_conv_i = ThermalResistor(R = R_conv_i)
@named Resistor_wall = ThermalResistor(R = R_wall)

# Heat port aka temperature source
@named T_weather = HeatPort(T = 273.15 + -200)


# Temeprature sensors
@named T_wall_out = TemperatureSensor()
@named T_wall_in = TemperatureSensor()
@named T_zone = TemperatureSensor()


connections = [
    connect(Resistor_conv_o.port_a, T_weather),
    
    connect(Resistor_conv_o.port_b, Capacitor_wall_out.port, Resistor_wall.port_a, T_wall_out.port),

    connect(Resistor_wall.port_b, Capacitor_wall_in.port, Resistor_conv_i.port_a, T_wall_in.port),

    connect(Resistor_conv_i.port_b, Capacitor_zone.port, T_zone.port),

]

@named model = ODESystem(connections, t,
    systems = [Capacitor_wall_out, Capacitor_wall_in, Capacitor_zone, Resistor_conv_o,Resistor_conv_i, Resistor_wall, T_weather, T_wall_out, T_wall_in, T_zone])
sys = structural_simplify(model)
prob = ODEProblem(sys, Pair[], (0, 3600))
sol = solve(prob, Tsit5())


Plot1 = plot(title = "Thermal Conduction Case")
Plot1 = plot!(sol, idxs = [Capacitor_wall_out.T, Capacitor_wall_in.T, Capacitor_zone.T].-273.15, labels = ["Wall out" "Wall in" "Zone"])
display(Plot1)

contradict · February 22, 2024, 5:16pm

A HeatPort is a way to connect components, not a component itself. Your simple case with a fixed temperature should make these two changes.

---    @named T_weather = HeatPort(T = 273.15 + -200)
+++    @named T_weather = FixedTemperature(T = 273.15 + -200)


---        connect(Resistor_conv_o.port_a, T_weather),
+++        connect(Resistor_conv_o.port_a, T_weather.port),

It look like there is no existing component to implement a time-dependent temperature, but you can do something like this:

function TabulatedTemperature(hourly_data; name)
    @named port = HeatPort()
    quations = [
         # t in seconds from your tspan in prob, adjust here if that is not what you meant.
        port.T ~ hourly_data[Int(floor(t) % 3600 + 1)]
    ]
    ODESystem(equations, t; systems=[port], name)
end

Then you can do

    @named T_weather = TabulatedTemperature(T_out_hourly)

and connect as for FixedTemperature.

If you want interpolated response rather than hourly step response, you can replace indexing with an interpolator from DataInterpolations.jl

Zhijian · February 23, 2024, 8:20am

Thank you ! That indeed works!

Zhijian · February 23, 2024, 9:16am

Hi I do have one more issue

Now I am implementing the hourly data in that way

function TabulatedTemperature(hourly_data::Vector{Float64}; name)
    @named port = HeatPort()
    equations = [
         # t in seconds
        port.T ~ hourly_data[Int(floor(t/3600) % 24 + 1.0)]
    ]
    ODESystem(equations, t; systems=[port], name)
end

but it reported：

MethodError: no method matching Int64(::Num)

When I change this it can run:

---         port.T ~ hourly_data[Int(floor(t/3600) % 24 + 1.0)]
+++      port.T ~ hourly_data[2]

I feel quite comfused about the data type of this part, because I used to have similar code in ODEs with the t, and it can just run.

Here is the full code:

using ModelingToolkit
using ModelingToolkitStandardLibrary.Thermal, ModelingToolkit, OrdinaryDiffEq, Plots

@parameters t


A = 20.0 #wall area
h = 1.0 #Convective coefficient
C_wall_out = 20000
C_wall_in = 20000
C_zone = 70000

R_conv_o = 1/(R*A)
R_conv_i = 1/(R*A)
R_wall = 0.13
#TempConstant = 20

#To use hourly data in the model

function TabulatedTemperature(hourly_data::Vector{Float64}; name)
    @named port = HeatPort()
    equations = [
         # t in seconds
        port.T ~ hourly_data[Int(floor(t/3600) % 24 + 1.0)]
    ]
    ODESystem(equations, t; systems=[port], name)
end

# Capacitors
@named Capacitor_wall_out = HeatCapacitor(C = C_wall_out, T = 273.15 + 20)
@named Capacitor_wall_in = HeatCapacitor(C = C_wall_in, T = 273.15 + 50)
@named Capacitor_zone = HeatCapacitor(C = C_zone, T = 273.15 + 50)

# Resistors
@named Resistor_conv_o = ThermalResistor(R = R_conv_o)
@named Resistor_conv_i = ThermalResistor(R = R_conv_i)
@named Resistor_wall = ThermalResistor(R = R_wall)

# Heat port aka temperature source
T_out_hourly = [0, 0, 0, 0, 0, 0, 5, 5, 5, 5, 5, 15, 15, 15, 15, 15, 15, 5, 5, 5, 0, 0, 0, 0].+273.15
@named T_weather = TabulatedTemperature(T_out_hourly)


# Temeprature sensors
@named T_amb = TemperatureSensor()
@named T_wall_out = TemperatureSensor()
@named T_wall_in = TemperatureSensor()
@named T_zone = TemperatureSensor()


connections = [
    connect(Resistor_conv_o.port_a, T_weather.port,T_amb.port),
    
    connect(Resistor_conv_o.port_b, Capacitor_wall_out.port, Resistor_wall.port_a, T_wall_out.port),

    connect(Resistor_wall.port_b, Capacitor_wall_in.port, Resistor_conv_i.port_a, T_wall_in.port),

    connect(Resistor_conv_i.port_b, Capacitor_zone.port, T_zone.port),

]

@named model = ODESystem(connections, t,
    systems = [Capacitor_wall_out, Capacitor_wall_in, Capacitor_zone, Resistor_conv_o,Resistor_conv_i, Resistor_wall, T_weather, T_wall_out, T_wall_in, T_zone])
sys = structural_simplify(model)
prob = ODEProblem(sys, Pair[], (0, 36000))
sol = solve(prob, Tsit5())


Plot1 = plot(title = "Thermal Conduction Case")
Plot1 = plot!(sol, idxs = [Capacitor_wall_out.T, Capacitor_wall_in.T, Capacitor_zone.T].-273.15, labels = ["Wall out" "Wall in" "Zone"])
display(Plot1)

contradict · February 23, 2024, 4:35pm

Oops, sorry. Looks like symbolic tracing doesn’t work with indexing. You can tell symbolics not to trace that function and just treat it like a black box.

# This is the part symbolics shouldn't try to trace through
function hourly_index(temps, t)
    temps[Int(floor(t/3600) % 24 + 1)]
end

# register it here,
# Need the type annotation here to request a method for a function which
# accepts a vector, otherwise it only generates scalar methods.
@register_symbolic hourly_index(temps::Vector{Float64}, t)

function TabulatedTemperature(hourly_data; name)
    @parameters t

    @named port = HeatPort()
    equations = [
         # t in seconds
         port.T ~ hourly_index(hourly_data, t)
    ]
    ODESystem(equations, t; systems=[port], name)
end