Media / Numerical solving for 1D thermal modeling with Modelingtoolkit

New to Julia
,
1

Hello everyone,

I’m a master’s course student just starting with Julia and ModelingToolkit. I have a few questions for those experienced in modeling and simulating thermal systems (like Vapor Compression Cycles or ORCs) with these tools.

My research goal is to simulate the energy consumption and COP degradation in a reversible heat pump by developing a 1D dynamic model that includes frost and defrost cycles, referencing the work of Ma et al. (2023). As a first step, I’m working on a simpler 1D plate heat exchanger (PHX) model (refrigerant-to-water) to validate the approach. I’ve run into some challenges regarding media property calculations and numerical stability.

FCV tutorial from review article [Rasmussen, 2012]
FCV tutorial from review article [Rasmussen, 2012]1178×476 46.8 KB

  1. Media Properties: Lookup Tables/ANNs
    Currently, I’m calling CoolProp’s low-level interface (BICUBIC&HEOS) within each control volume to calculate fluid properties from pressure (P) and enthalpy (h). While this works, I’ve noticed that compared to my experience with MATLAB/Simulink, where lookup tables and interpolation offered faster performance, calling CoolProp directly introduces considerable computational overhead. I recently read in Ma et al. (2023) that they employed a Physics-Informed Neural Network (PINN) as a substitute for a lookup table. This approach seems highly compatible with the automatic differentiation capabilities of Julia’s solvers. I’m now planning to replace my current property calculation method. From your experience, which approach—a traditional lookup table or an ANN—integrates more effectively and performs better within the ModelingToolkit ecosystem?
# CoolProp method
import CoolProp as CP
# CoolProp with MTK (Tabular / Consider to change into PINN)

const handle = CP.AbstractState_factory("HEOS&BICUBIC", "R410A")

# ρ, T from P, h
function _R410a_ρ(P::Float64, h::Float64, handle = handle)
    CP.AbstractState_update(handle, CP.get_input_pair_index("HmassP_INPUTS"), h, P)   
    return CP.AbstractState_keyed_output(handle, CP.get_param_index("Dmass"))
end
function _R410a_T(P, h, handle=handle)
    CP.AbstractState_update(handle, CP.get_input_pair_index("HmassP_INPUTS"), h, P)   
    return CP.AbstractState_keyed_output(handle, CP.get_param_index("T"))
end

function _R410a_L(P, handle=handle)
    CP.AbstractState_update(handle, CP.get_input_pair_index("PQ_INPUTS"), P, 0)
    return CP.AbstractState_keyed_output(handle, CP.get_param_index("Hmass")), 
    CP.AbstractState_keyed_output(handle, CP.get_param_index("Dmass")), 
    CP.AbstractState_keyed_output(handle, CP.get_param_index("L")),
    CP.AbstractState_keyed_output(handle, CP.get_param_index("V")),
    CP.AbstractState_keyed_output(handle, CP.get_param_index("Prandtl"))
end
function _R410a_V(P, handle=handle)
    CP.AbstractState_update(handle, CP.get_input_pair_index("PQ_INPUTS"), P, 1)
    return CP.AbstractState_keyed_output(handle, CP.get_param_index("Hmass")),
    CP.AbstractState_keyed_output(handle, CP.get_param_index("Dmass"))
end

function _R410a_dρdP_h(P,h, handle=handle) 
    return _R410a_ρ(P+1, h) - _R410a_ρ(P,h)
end
function _R410a_dρdh_P(P,h, handle=handle)
    return _R410a_ρ(P, h+1) - _R410a_ρ(P,h)
end

R410a_ρ(P, h) =_R410a_ρ(P,h)
R410a_T(P,h) = _R410a_T(P,h)
R410a_L(P) = _R410a_L(P)
R410a_V(P) = _R410a_V(P)
R410a_dρdP_h(P,h) = _R410a_dρdP_h(P,h)
R410a_dρdh_P(P,h) = _R410a_dρdh_P(P,h)

# Symbolics
@register_symbolic R410a_ρ(P, h) 
@register_symbolic R410a_T(P, h)
@register_symbolic R410a_dρdP_h(P,h)
@register_symbolic R410a_dρdh_P(P,h)

# Usage Example
# (Ex). Func for h_flow
# Qiao et al. (2015): Transient modeling of a flash tank vapor injection heat pump system - Part I: Model development
function h_flow_func(P, h)
    h_L, ρ_L, _, _, _ = R410a_L(P)
    h_V, ρ_V = R410a_V(P)
    x = (h - h_L) / (h_V - h_L)
    # Single phase
    if abs(x - 0.5) >= 0.5   
        return h
    end
    ρ = R410a_ρ(P, h)

    # Two-phase
    γ = (ρ - ρ_L) / (ρ_V - ρ_L) 
    # Zivi (1964) : c = 1, q = 1, r = (2/3), s = 0
    # c, q, r = 1.0, 1.0, 2/3
    # Smith (1969) : c = 0.79, q = 0.78, r = 0.58, s = 0
    c, q, r = 0.79, 0.78, 0.58
    x_flow = γ^(1/q) / (γ^(1/q) + ((1/c)*(ρ_L / ρ_V)^(r))^(1/q) * (1 - γ)^(1/q) )

    h_flow = h + (x_flow - x) * (h_V - h_L)
    return h_flow
end
  1. Numerical Stability: Implicit vs. Explicit Formulation
    I’ve learned that minimizing the number of symbolic variables is crucial for numerical performance in ModelingToolkit. In one version of my model, solvers like QNDF and Rodas5p (which I understand are similar to MATLAB’s ode15s/ode23s) perform reasonably well. However, in a second version where I explicitly formulate the equations with derivative terms such as D(P) and D(h[i]), the solvers fail to converge. This issue didn’t seem to arise in Simulink examples, likely because I manually eliminate algebraic constraints by expressing variables like mdot[i] as a function of D(P), D(h[i]), mdot_in and mdot_out, using the mass conservation equations. I’m struggling to understand the underlying numerical difference between these two formulations. Why is the more implicit (less symbolically manipulated) formulation solvable, while the more explicit one is not?
# Solvable with QNDF/Rodas5P [reltol = 1e-3] / Implicit form
HX_Ref_V_elem * (R410a_dρdP_h(P, h[i]) * D(P) + R410a_dρdh_P(P, h[i]) * D(h[i])) ~ mdot[i] - mdot[i+1],
HX_Ref_V_elem * ((h[i] * R410a_dρdP_h(P, h[i]) - 1) * D(P) + (h[i] * R410a_dρdh_P(P, h[i]) + R410a_ρ(P, h[i])) * D(h[i])) ~ mdot[i] * (i == 1 ? h_flow_func(P_in, h_init) : h_flow_func(P, h[i-1])) - mdot[i+1] * h_flow_func(P, h[i]) + Q_r[i],

I. Implicit form with Dummy variables after ODESystem
I. Implicit form with Dummy variables after ODESystem1739×412 69.2 KB

# Failed version (Using inv([A B;C D]) == 1/(AD - BC) * [D -B; -C A] from implicit form)
HX_Ref_V_elem * D(P) ~ ((mdot[i] - mdot[i+1]) * (h[i] * R410a_dρdh_P(P, h[i]) + R410a_ρ(P, h[i])) - (mdot[i] * (i == 1 ? h_flow_func(P_in, h_init) : h_flow_func(P, h[i-1])) - mdot[i+1] * h_flow_func(P, h[i]) + Q_r[i]) * R410a_dρdh_P(P, h[i])) / (R410a_ρ(P, h[i]) * R410a_dρdh_P(P, h[i]) - R410a_dρdP_h(P, h[i])),
HX_Ref_V_elem * D(h[i]) ~ (-(mdot[i] - mdot[i+1]) * (h[i] * R410a_dρdh_P(P, h[i]) - 1) + (mdot[i] * (i == 1 ? h_flow_func(P_in, h_init) : h_flow_func(P, h[i-1])) - mdot[i+1] * h_flow_func(P, h[i]) + Q_r[i]) *  R410a_dρdP_h(P, h[i])) / (R410a_ρ(P, h[i]) * R410a_dρdh_P(P, h[i]) - R410a_dρdP_h(P, h[i])),

II. Explicit form with P, h[i]
II. Explicit form with P, h[i]1765×511 173 KB

[Edit, 250901] Since the eqns look a little bit complex, here’s a simplified version for when N_CV is 5.
[Successful case]

# der(y) == dy/dt
A(t)*der(P) + B(t)*der(h[i]) ~ b1(t)
C(t)*der(P) + D(t)*der(h[i]) ~ b2(t)

[Failed case]

der(P) ~ D(t)*b1(t) - B(t)*b2(t) / (A(t)*D(t) - B(t)*C(t))
der(h[i]) ~ b2(t) ~ C(t)*b1(t) - A(t)*b2(t) / (A(t)*D(t) - B(t)*C(t))

2-1. Acausal Modeling with Connectors
I’m considering refactoring my model using the acausal component and connector style (from MTK or a library like Dyad) to separate volume states (P, h, ρ) from flow states (mdot, H_in, H_out). Would the primary benefit of this be improved modularity and ease of model construction, or could this component-based approach also positively impact the numerical stability issues mentioned in my second question? Though I suppose that’s a question best answered by trying it myself. I had a feeling that it would be a great help in implementing a staggered grid scheme.

As a mechanical engineering student, I’m aware that my background in numerical analysis is limited, so I appreciate any insights you can offer. I’m excited about the potential of these tools and if my experience is positive, I’ll definitely encourage my colleagues to consider them as an alternative to Simulink.

Thank you for your time, and I look forward to your response.
Below is my MTK code for solving.

# 1. Params, Variables for Ref
# HX_C::HX_PHX: immutable struct that contains geometry information of PHX
N = HX_C.N_CV
U_Fld_elem = U_func_Fld(HX_C)

# Init Props
P_in_val = 3200e3
P_out_val = 3000e3
h_in_val = 450e3
c1_val = 9e-6
c2_val = c1_val
Δt_val = 3e-3  #[s]

# 1. Params for HX
# 1-0. Init cond
@parameters begin
    c1=c1_val 
    c2=c2_val
    P_in=P_in_val
    P_out=P_out_val
    h_init=h_in_val
    Δt = Δt_val  #[s]
end

# 1-1. Ref side
@parameters begin
    HX_Ref_A_elem = HX_C.Ref.A_elem
    HX_Ref_V_elem = HX_C.Ref.V_elem
end
# 1-2. Wall Side
@parameters begin
    HX_Wall_m_elem = HX_C.Wall.m_elem
    HX_Wall_Cp = HX_C.Wall.Cp
end
# 1-3. Fld Side
@parameters begin
    HX_Fld_mdot_in = HX_C.Fld.mdot
    HX_Fld_T_in = HX_C.Fld.T_in
    HX_Fld_A_elem = HX_C.Fld.A_elem
    HX_Fld_Cp = HX_C.Fld.Cp
    UA_Fld_elem = U_Fld_elem * HX_C.Fld.A_elem # UA const for Fld side
end

# 2. Variables for HX
# 2-1. Ref Side
    # Static
@variables begin
    P(t)
    h(t)[1:N]
    UA_elem(t)[1:N]
    UA_elem_hat(t)[1:N]
end
    # Flow
@variables begin
    mdot(t)[1:N+1]
    Q_r(t)[1:N]
end
# 2-2. Wall Side
@variables T_wall(t)[1:N]
# 2-3. Fld Side
    # Static
@variables begin
    T_Fld(t)[1:N+1]
end
    # Flow
@variables begin
    Q_f(t)[1:N]
end


# 3. Equations
# 3-1. Cond
# 3-1.1 Cond_Ref_side
# Const P method (no friction)
cond_Ref = [
    [
    #= Algebraic eqns for Ref
    ρ[i] == R410a_ρ(P, h[i]),
    T[i] == R410a_T(P, h[i]),
    =#
    #=
    H_in[i] ~ mdot[i] * (i == 1 ? h_flow_func(P_in, h_init) : h_flow_func(P, h[i-1])),
    H_out[i] ~ mdot[i+1] * h_flow_func(P, h[i]),
    =#

    #= Flow vars cal(Q_r[i] : By Static vars (P, h[i], ρ[N]), T)
    (Qiao et al., 2018) On Closure Relations for Dynamic Vapor Compression Cycle models
    Low Pass Filtering of UA_elem
    =#
    UA_elem_hat[i] ~ U_func(P, h[i], mdot[i]) * HX_Ref_A_elem,
    D(UA_elem[i]) ~ (UA_elem_hat[i] - UA_elem[i]) / Δt,
    Q_r[i] ~ UA_elem[i] * (T_wall[i] - R410a_T(P, h[i])),
    
    #=Differential eqns for Ref [Static vars diff by flow vars(H_in, H_out, Q_r)]
    (1) Mass conservation
    HX_Ref_V_elem * D(ρ[i]) == mdot_in[i] - mdot_out[i],
    (2) Energy conservation
    HX_Ref_V_elem * (D(ρ[i] * h[i]) - D(P[i])) == H_in[i] - H_out[i] + Q_r[i]
    (3) Partial eqns
    D(ρ[i]) == dρdP_h(P, h[i]) * D(P) + dρdh_P(P, h[i]) * D(h[i])
    =#

    #= (Failed) why? explicit vs implicit?
    A == R410a_dρdP_h(P, h[i])
    B == R410a_dρdh_P(P, h[i])
    C == (h[i] * R410a_dρdh_P(P, h[i]) - 1)
    D == (h[i] * R410a_dρdh_P(P, h[i]) + R410a_ρ(P, h[i]))
    b1 == (mdot[i] - mdot[i+1])
    b2 == (H_in[i] - H_out[i] + Q_r[i])
    AD - BC == (R410a_ρ(P, h[i]) * R410a_dρdh_P(P, h[i]) - R410a_dρdP_h(P, h[i]))
    HX_Ref_V_elem * D(P) == (b1*D - b2*B) / (A*D - B*C)
    HX_Ref_V_elem * D(h[i]) == (-b1*C + b2*A) / (A*D - B*C)
    
    HX_Ref_V_elem * D(P) ~ ((mdot[i] - mdot[i+1]) * (h[i] * R410a_dρdh_P(P, h[i]) + R410a_ρ(P, h[i])) - (mdot[i] * (i == 1 ? h_flow_func(P_in, h_init) : h_flow_func(P, h[i-1])) - mdot[i+1] * h_flow_func(P, h[i]) + Q_r[i]) * R410a_dρdh_P(P, h[i])) / (R410a_ρ(P, h[i]) * R410a_dρdh_P(P, h[i]) - R410a_dρdP_h(P, h[i])),
    HX_Ref_V_elem * D(h[i]) ~ (-(mdot[i] - mdot[i+1]) * (h[i] * R410a_dρdh_P(P, h[i]) - 1) + (mdot[i] * (i == 1 ? h_flow_func(P_in, h_init) : h_flow_func(P, h[i-1])) - mdot[i+1] * h_flow_func(P, h[i]) + Q_r[i]) *  R410a_dρdP_h(P, h[i])) / (R410a_ρ(P, h[i]) * R410a_dρdh_P(P, h[i]) - R410a_dρdP_h(P, h[i])),
    =#
    
    HX_Ref_V_elem * (R410a_dρdP_h(P, h[i]) * D(P) + R410a_dρdh_P(P, h[i]) * D(h[i])) ~ (mdot[i] - mdot[i+1]),
    HX_Ref_V_elem * ((h[i] * R410a_dρdP_h(P, h[i]) - 1) * D(P) + (h[i] * R410a_dρdh_P(P, h[i]) + R410a_ρ(P, h[i])) * D(h[i])) ~ mdot[i] * (i == 1 ? h_flow_func(P_in, h_init) : h_flow_func(P, h[i-1])) - mdot[i+1] * h_flow_func(P, h[i]) + Q_r[i],
    ]
    for i in 1:N
]
# sqrt, regSqrt options for mdot_initialization / Orifice valve model inputs
push!(cond_Ref, [mdot[1] ~ ifelse(P_in - P > 0.0, c1 * (regSqrt(R410a_ρ(P, h_init) * (P_in - P))), 0.0)])
push!(cond_Ref, [mdot[N+1] ~ ifelse(P - P_out > 0.0, c2 * (regSqrt(R410a_ρ(P, h[N]) * (P - P_out))), 0.0)])

# 3-1.2 Cond_Wall_Fld_side
if HX_C.isCounterFlow == true
    #  Cond_Wall_side
    cond_Wall = [
        [
            #  Algebraic eqns for Wall
            #  Differential eqns for Wall
            HX_Wall_m_elem * HX_Wall_Cp * D(T_wall[i]) ~ Q_f[N+1-i] - Q_r[i]
        ]
        for i in 1:N
    ]
    #  Cond_Fld_side
    cond_Fld = [
        [
            #  Algebraic eqns for Fld
            Q_f[N+1-i] ~ (HX_Fld_mdot_in * HX_Fld_Cp) * (1 - exp(-UA_Fld_elem / (HX_Fld_mdot_in * HX_Fld_Cp))) * (T_Fld[N+1-i] - T_wall[i]),
            Q_f[N+1-i] ~ (HX_Fld_mdot_in * HX_Fld_Cp) * (T_Fld[N+1-i] - T_Fld[N+2-i])
        ]
        for i in 1:N
    ]
else
    #  Cond_Wall_side
    cond_Wall = [
        [
            #  Algebraic eqns for Wall
            #  Differential eqns for Wall
            HX_Wall_m_elem * HX_Wall_Cp * D(T_wall[i]) ~ Q_f[i] - Q_r[i]
        ]
        for i in 1:N
    ]
    # Cond_Fld_side
    cond_Fld = [
        [
            #  Algebraic eqns for Fld
            Q_f[i] ~ (HX_Fld_mdot_in * HX_Fld_Cp) * (1 - exp(-UA_Fld_elem / (HX_Fld_mdot_in * HX_Fld_Cp))) * (T_Fld[i] - T_wall[i]),
            Q_f[i] ~ (HX_Fld_mdot_in * HX_Fld_Cp) * (T_Fld[i+1] - T_Fld[i])
        ]
        for i in 1:Nc
    ]
end
push!(cond_Fld, [T_Fld[1] ~ HX_Fld_T_in])


# Flatten the cond models
cond_eqns = vcat(cond_Ref..., cond_Wall..., cond_Fld...)
@mtkcompile sys = ODESystem(cond_eqns, t; continuous_events = [P - P_out ~ 0])

# Init
mdot_in_val = c1_val * regSqrt(R410a_ρ(P_in_val, h_in_val) * (P_in_val - P_out_val))
U_init= U_func(P_in_val, h_in_val, mdot_in_val)

u0_P = [P => P_in_val]
u0_h = [h[i] => h_in_val for i in 1:N]
u0_T_wall = [T_wall[i] => HX_C.Fld.T_in for i in 1:N]
u0_UA_elem = [UA_elem[i] => U_init * HX_C.Ref.A_elem for i in 1:N]
u0_map = vcat(u0_P, u0_h, u0_T_wall, u0_UA_elem)

# Guesses
guess_mdot = [mdot[i] => mdot_in_val for i in 1:N+1]
guess_map = vcat(guess_mdot)

# ODEProblem
tspan = (0.0, 10.0)
prob = ODEProblem(sys, u0_map, tspan, guesses = guess_map)

# For ODEProbelm Rosenbrock23, Rodas5P -> FBDF, QNDF  -> AdaptiveRadau [Stiff ODE solvers]
sol = solve(prob, QNDF(autodiff=false))

Reference:
Dynamic modeling for vapor compression systems—Part II: Simulation tutorial, [Rasmussen et al, 2012]
Transient modeling of a flash tank vapor injection heat pump system – Part I: Model development, [Qiao et al, 2015]
Development and validation of a dynamic modeling framework for air-source heat pumps under cycling of frosting and reverse-cycle defrosting, [Ma et al., 2023]