SciML / ModelingToolkitStandardLibrary.jl

A standard library of components to model the world and beyond
https://docs.sciml.ai/ModelingToolkitStandardLibrary/stable/
MIT License
112 stars 36 forks source link

Thermal.FixedHeatFlow cannot actually model a fixed heat flow rate source #234

Open hmdnt opened 10 months ago

hmdnt commented 10 months ago

Here is an MWE that shows this:

using DifferentialEquations, ModelingToolkit
using ModelingToolkitStandardLibrary: Thermal
@variables t
@named temp = Thermal.FixedTemperature(; T=300)
@named heatflow = Thermal.FixedHeatFlow(; Q_flow=-1.0)
@named wall = Thermal.ThermalResistor(; R=1)
eqs = [
    connect(temp.port, wall.port_a),
    connect(wall.port_b, heatflow.port),
]
@named test_model = compose(ODESystem(eqs, t; name=:_test), [temp, heatflow, wall])
sys = structural_simplify(test_model)
sol = solve(ODAEProblem(sys, [], (0, 1.0)))
sol[wall.Q_flow]

The last line gives the following error:

ERROR: The nonlinear solver failed with the return code MaxIters.
...

The problem stems from the temperature dependence of the heat flow implemented in Thermal.FixedHeatFlow. observed(sys) gives the following list for the observed variables.

julia> observed(sys)
9-element Vector{Equation}:
 temp₊port₊T(t) ~ temp₊T
 heatflow₊port₊Q_flow(t) ~ wall₊Q_flow(t)
 wall₊port_a₊Q_flow(t) ~ wall₊Q_flow(t)
 wall₊port_b₊Q_flow(t) ~ -wall₊Q_flow(t)
 temp₊port₊Q_flow(t) ~ -wall₊Q_flow(t)
 wall₊port_a₊T(t) ~ temp₊port₊T(t)
 heatflow₊port₊T(t) ~ (-heatflow₊port₊Q_flow(t) + heatflow₊Q_flow*(-1 + heatflow₊T_ref*heatflow₊alpha)) / (heatflow₊Q_flow*heatflow₊alpha)
 wall₊dT(t) ~ temp₊port₊T(t) - heatflow₊port₊T(t)
 wall₊port_b₊T(t) ~ heatflow₊port₊T(t)

The culprit is the 7th equation in which alpha is in denominator and for constant heat source is zero. Also if the source is modeling an insulated point where Q_flow = 0, the same symptom occurs regardless of the value of alpha.

The remedy is to correctly define the Thermal.FixedHeatFlow equation as port.Q_flow ~ -Q_flow and ditch the temperature dependence, which is meaningless for a fixed heat flow source anyway. Note that the same argument applies to Thermal.PrescribedHeatFlow for which the correct equation should be port.Q_flow ~ -Q_flow.u.

PS. I know that Modelica Standard Library defines the above sources the same way that they are currently implemented in ModelingToolkit Standard Library, but I think the correct way to define them is as discussed above. When a temperature dependent heat flow source is desired, it can be achieved by combining Thermal.PrescribedHeatFlow, Thermal.TemperatureSensor and some function block.

ValentinKaisermayer commented 10 months ago

This is from the Modelica stlib https://github.com/modelica/ModelicaStandardLibrary/blob/d3ad60c312774a1bf84a74e830ccd55abcfb09d2/Modelica/Thermal/HeatTransfer/Sources/FixedHeatFlow.mo#L12

ValentinKaisermayer commented 10 months ago

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 with respect to reference temperature T_ref).

ValentinKaisermayer commented 10 months ago

In Modelica it works with alpha=0. Here, for alpha=0, the correct solution is that the 7th equation should no be here at all, and wall.port_b.T has to be determined from the heat flow, the resistance and wall.port_a.T.

I see this more as an issue with the observed generation.

hmdnt commented 10 months ago

Yes I know that MSL defines this the way it is. May be as you said there is a problem with observed generation. But defining a temperature dependent heat flow for a fixed heat source is illogical unless MTSL is going to be an exact replica of MSL.