Merge pull request #73 from ValentinKaisermayer/fix-der_T

Make thermal models ODEs again

Author: ChrisRackauckas

docs/pages.jl

@@ -4,6 +4,7 @@ pages = [
         "Tutorials" => [
             "RC Circuit" => "tutorials/rc_circuit.md",
             "Custom Components" => "tutorials/custom_component.md",
+            "Thermal Conduction Model" => "tutorials/thermal_model.md",
         ],
 
         "API" => [

docs/src/API/thermal.md

@@ -43,5 +43,7 @@ TemperatureSensor
 
 ```@docs
 FixedHeatFlow
-FixedTemperature 
+FixedTemperature
+PrescribedHeatFlow
+PrescribedTemperature  
 ```

docs/src/index.md

@@ -16,6 +16,7 @@ import Pkg; Pkg.add("ModelingToolkitStandardLibrary")
 
 - [RC Circuit](http://mtkstdlib.sciml.ai/dev/tutorials/rc_circuit/)
 - [Custom Component](http://mtkstdlib.sciml.ai/dev/tutorials/custom_component/)
+- [Thermal Model](http://mtkstdlib.sciml.ai/dev/tutorials/thermal_model/)
 
 ## Libraries
 

docs/src/tutorials/thermal_model.md

@@ -0,0 +1,41 @@
+# Heat Conduction Model
+
+This example demonstrates the thermal response of two masses connected by a conducting element. 
+The two masses have the same heat capacity but different initial temperatures (T1=100 [°C], T2=0 [°C]). 
+The mass with the higher temperature will cool off while the mass with the lower temperature heats up. 
+They will each asymptotically approach the calculated temperature T_final_K that results 
+from dividing the total initial energy in the system by the sum of the heat capacities of each element.
+
+```julia
+using ModelingToolkitStandardLibrary.Thermal, ModelingToolkit, OrdinaryDiffEq, Plots
+
+@parameters t
+
+C1 = 15
+C2 = 15
+@named mass1 = HeatCapacitor(C=C1, T_start=373.15)
+@named mass2 = HeatCapacitor(C=C2, T_start=273.15)
+@named conduction = ThermalConductor(G=10)
+@named Tsensor1 = TemperatureSensor() 
+@named Tsensor2 = TemperatureSensor()
+
+connections = [
+    connect(mass1.port, conduction.port_a),
+    connect(conduction.port_b, mass2.port),
+    connect(mass1.port, Tsensor1.port),
+    connect(mass2.port, Tsensor2.port),
+]
+
+@named model = ODESystem(connections, t, systems=[mass1, mass2, conduction, Tsensor1, Tsensor2])
+sys = structural_simplify(model)
+prob = ODEProblem(sys, Pair[], (0, 5.0))
+sol = solve(prob, Tsit5())
+
+T_final_K = sol[(mass1.T * C1 + mass2.T * C2) / (C1 + C2)]
+
+plot(title = "Thermal Conduction Demonstration")
+plot!(sol, vars = [mass1.T, mass2.T], labels = ["Mass 1 Temperature" "Mass 2 Temperature"])
+plot!(sol.t, T_final_K, label = "Steady-State Temperature")
+savefig("thermal_plot.png")
+```
+![Plot of Temperatures](https://user-images.githubusercontent.com/50108075/173061172-4c7305f5-1193-4b17-9ca4-8f8dcbc0cae7.png)

src/Mechanical/Rotational/Rotational.jl

@@ -4,7 +4,6 @@ Library to model 1-dimensional, rotational mechanical systems
 module Rotational
 
 using ModelingToolkit, Symbolics, IfElse, OrdinaryDiffEq
-using OffsetArrays
 using ...Blocks: RealInput, RealOutput
 
 @parameters t

src/Thermal/HeatTransfer/ideal_components.jl

@@ -3,13 +3,16 @@
 
 Lumped thermal element storing heat
 
-# Parameters:
-- `C`: [J/K] Heat capacity of element (= cp*m)
-- `T_start`: Initial temperature of element
-
 # States:
-- `T`: [K] Temperature of element
-- `der_T`: [K/s] Time derivative of temperature
+- `T`: [`K`] Temperature of element
+- `der_T`: [`K/s`] Time derivative of temperature
+
+# Connectors:
+- `port`
+
+# Parameters:
+- `C`: [`J/K`] Heat capacity of element (= cp*m)
+- `T_start`: [`K`] Initial temperature of element
 """
 function HeatCapacitor(; name, C, T_start=273.15 + 20)    
     @named port = HeatPort()
@@ -22,8 +25,8 @@ function HeatCapacitor(; name, C, T_start=273.15 + 20)
     D = Differential(t)
     eqs = [
         T ~ port.T
-        der_T ~ D(T)
-        D(T) ~ port.Q_flow / C
+        der_T ~ port.Q_flow / C
+        D(T) ~ der_T
     ]
     ODESystem(eqs, t, sts, [C]; systems=[port], name=name)
 end
@@ -33,8 +36,15 @@ end
 
 Lumped thermal element transporting heat without storing it.
 
+# States:
+see [`Element1D`](@ref)
+
+# Connectors:
+`port_a`
+`port_b`
+
 # Parameters:
-- `G`: [W/K] Constant thermal conductance of material
+- `G`: [`W/K`] Constant thermal conductance of material
 """
 function ThermalConductor(;name, G)   
     @named element1d = Element1D()
@@ -51,8 +61,16 @@ end
 
 Lumped thermal element transporting heat without storing it.
 
+# States:
+- `dT`:  [`K`] Temperature difference across the component a.T - b.T
+- `Q_flow`: [`W`] Heat flow rate from port a -> port b
+
+# Connectors:
+- `port_a`
+- `port_b`
+
 # Parameters:
-- `R`: [K/W] Constant thermal resistance of material
+- `R`: [`K/W`] Constant thermal resistance of material
 """
 function ThermalResistor(; name, R)   
     @named element1d = Element1D()
@@ -70,12 +88,16 @@ end
 
 Lumped thermal element for heat convection.
 
+# States:
+- `dT`:  [`K`] Temperature difference across the component solid.T - fluid.T
+- `Q_flow`: [`W`] Heat flow rate from `solid` -> `fluid`
+
+# Connectors:
+- `solid`
+- `fluid`
+
 # Parameters:
 - `G`: [W/K] Convective thermal conductance
-
-# States:
-- `dT`:  [K] Temperature difference across the component solid.T - fluid.T
-- `Q_flow`: [W] Heat flow rate from solid -> fluid
 """
 function ConvectiveConductor(; name, G)
     @named solid = HeatPort()
@@ -96,32 +118,44 @@ end
 
 Lumped thermal element for heat convection.
 
-# Parameters:
-- `R`: [K/W] Constant thermal resistance of material
-
 # States:
-- `dT`:  [K] Temperature difference across the component solid.T - fluid.T
-- `Q_flow`: [W] Heat flow rate from solid -> fluid
+- `dT`:  [`K`] Temperature difference across the component solid.T - fluid.T
+- `Q_flow`: [`W`] Heat flow rate from `solid` -> `fluid`
+
+# Connectors:
+- `solid`
+- `fluid`
+
+# Parameters:
+- `R`: [`K/W`] Constant thermal resistance of material
 """
 function ConvectiveResistor(; name, R)
-    @named solidport = HeatPort()
-    @named fluidport = HeatPort()
+    @named solid = HeatPort()
+    @named fluid = HeatPort()
     @parameters R=R
     sts = @variables Q_flow(t)=0.0 dT(t)=0.0
     eqs = [
-        dT ~ solidport.T - fluidport.T
-        solidport.Q_flow ~ Q_flow
-        fluidport.Q_flow ~ -Q_flow
+        dT ~ solid.T - fluid.T
+        solid.Q_flow ~ Q_flow
+        fluid.Q_flow ~ -Q_flow
         dT ~ R*Q_flow
     ]
-    ODESystem(eqs, t, sts, [R]; systems=[solidport, fluidport], name=name)
+    ODESystem(eqs, t, sts, [R]; systems=[solid, fluid], name=name)
 end
 
 """
     BodyRadiation(; name, G)
 
 Lumped thermal element for radiation heat transfer.
 
+# States:
+- `dT`:  [`K`] Temperature difference across the component a.T - b.T
+- `Q_flow`: [`W`] Heat flow rate from port a -> port b
+
+# Connectors:
+- `port_a`
+- `port_b`
+
 # Parameters:
 - `G`: [m^2] Net radiation conductance between two surfaces
 """
@@ -145,6 +179,13 @@ end
 Collects `m` heat flows
 
 This is a model to collect the heat flows from `m` heatports to one single heatport.
+
+# States:
+
+# Connectors:
+- `port_a1` to `port_am`
+- `port_b`
+
 # Parameters:
 - `m`: Number of heat ports (e.g. m=2: `port_a1`, `port_a2`)
 """

src/Thermal/HeatTransfer/sensors.jl

@@ -6,10 +6,16 @@ Absolute temperature sensor in kelvin.
 This is an ideal absolute temperature sensor which returns the temperature of the connected port in kelvin as an output
 signal. The sensor itself has no thermal interaction with whatever it is connected to. Furthermore, no thermocouple-like 
 lags are associated with this sensor model.
+
+# States:
+- `T(t)`: [`K`] Absolute temperature
+
+# Connectors:
+- `port`
 """
 function TemperatureSensor(; name)
     @named port = HeatPort()
-    @variables T(t) # [K] Absolute temperature
+    @variables T(t)
     eqs = [
         T ~ port.T
         port.Q_flow ~ 0
@@ -24,11 +30,18 @@ Relative Temperature sensor.
 
 The relative temperature `port_a.T - port_b.T` is determined between the two ports of this component and is provided as 
 output signal in kelvin.
+
+# States:
+- `T(t)`: [`K`] Relative temperature `a.T - b.T`
+
+# Connectors:
+- `port_a`
+- `port_b`
 """
 function RelativeTemperatureSensor(; name)
     @named port_a = HeatPort()
     @named port_b = HeatPort()
-    @variables T(t) # [K] Relative temperature a.T - b.T
+    @variables T(t)
     eqs = [
         T ~ port_a.T - port_b.T
         port_a.Q_flow ~ 0
@@ -46,11 +59,18 @@ This model is capable of monitoring the heat flow rate flowing through this comp
 is the amount that passes through this sensor while keeping the temperature drop across the sensor zero. This is an ideal 
 model so it does not absorb any energy and it has no direct effect on the thermal response of a system it is included in.
 The output signal is positive, if the heat flows from `port_a` to `port_b`.
+
+# States:
+- `Q_flow(t)`: [`W`] Heat flow from `port_a` to `port_b`
+
+# Connectors:
+- `port_a`
+- `port_b`
 """
 function HeatFlowSensor(; name)
     @named port_a = HeatPort()
     @named port_b = HeatPort()
-    @variables Q_flow(t) # [W] Heat flow from port a to port b
+    @variables Q_flow(t) 
     eqs = [
         port_a.T ~ port_b.T
         port_a.Q_flow + port_b.Q_flow ~ 0

src/Thermal/HeatTransfer/sources.jl

@@ -7,6 +7,9 @@ This model allows a specified amount of heat flow rate to be "injected" into a t
 The constant amount of heat flow rate `Q_flow` is given as a parameter. The heat flows into the component to which
 the component FixedHeatFlow is connected, if parameter `Q_flow` is positive.
 
+# Connectors:
+- `port`
+
 # Parameters:
 - `Q_flow`: [W] Fixed heat flow rate at port
 - `T_ref`: [K] Reference temperature
@@ -31,7 +34,10 @@ end
 
 Fixed temperature boundary condition in kelvin.
 
-This model defines a fixed temperature T at its port in kelvin, i.e., it defines a fixed temperature as a boundary condition.
+This model defines a fixed temperature `T` at its port in kelvin, i.e., it defines a fixed temperature as a boundary condition.
+
+# Connectors:
+- `port`
 
 # Parameters:
 - `T`: [K] Fixed temperature boundary condition
@@ -44,3 +50,59 @@ function FixedTemperature(; name, T)
     ]
     ODESystem(eqs, t, [], pars; systems=[port], name=name)
 end
+
+
+"""
+    PrescribedHeatFlow(; name, Q_flow=1.0, T_ref=293.15, alpha=0.0)
+
+Prescribed heat flow boundary condition.
+
+This model allows a specified amount of heat flow rate to be "injected" into a thermal system at a given port. 
+The amount of heat is given by the input signal `Q_flow` into the model. The heat flows into the component to which 
+the component `PrescribedHeatFlow` is connected, if the input signal is positive.
+If parameter alpha is > 0, the heat flow is multiplied by `1 + alpha*(port.T - T_ref`) in order to simulate temperature 
+dependent losses (which are given an reference temperature T_ref).
+
+# Connectors:
+- `port`
+- `Q_flow` Input for the heat flow
+
+# Parameters:
+- `T_ref`: [K] Reference temperature
+- `alpha`: [1/K] Temperature coefficient of heat flow rate
+"""
+function PrescribedHeatFlow(; name, T_ref=293.15, alpha=0.0)
+    pars = @parameters begin 
+        T_ref=T_ref 
+        alpha=alpha
+    end
+    @named port = HeatPort()
+    @named Q_flow = RealInput()
+    
+    eqs = [
+        port.Q_flow ~ -Q_flow.u * (1 + alpha * (port.T - T_ref))
+    ]
+    ODESystem(eqs, t, [], pars; systems=[port, Q_flow], name=name)
+end
+
+"""
+    PrescribedTemperature(; name, T)
+
+This model represents a variable temperature boundary condition. 
+
+The temperature in kelvin is given as input signal `T` to the model. The effect is that an instance of 
+this model acts as an infinite reservoir able to absorb or generate as much energy as required to keep 
+the temperature at the specified value.
+
+# Connectors:
+- `port`
+- `T` input for the temperature
+"""
+function PrescribedTemperature(; name)
+    @named port = HeatPort()
+    @named T = RealInput()
+    eqs = [
+        port.T ~ T.u
+    ]
+    ODESystem(eqs, t, [], []; systems=[port, T], name=name)
+end