Hi, I’m trying to solve the energy/heat equation with temperature dependent material properties (density, enthalpy and thermal conductivity). Therefore I’m running in the issues that are mentioned as (current) known limitations for MethodOfLines.jl (I’m using version 0.2) MethodOfLines.jl: Automated Finite Difference for Physics-Informed Learning · MethodOfLines.jl
The original equation would look like:
eqs = [
Dt(ρ(T(t,x))*h(T(t,x))) ~ Dx(λ(T(t,x))*Dx(T(t,x)))
]
I’m looking for guidance and best practice how to apply the workarounds mentioned in the docs to my problem.
So far I tired using auxiliary variables:
eqs = [Dt(e(t,x)) ~ Dx(q(t,x)),
e(t,x) ~ ρ(T(t,x))*h(T(t,x)),
q(t,x) ~ λ(T(t,x))*Dx(T(t,x))]
Here I get errors that number of equations and variables are not matching.
The other stategy was applying chain and product rules:
eqs = [ DT(h(T(t,x)))*Dt(T(t,x))*ρ(T(t,x)) + DT(ρ(T(t,x)))*Dt(T(t,x))*h(T(t,x)) ~ λ(T(t,x))*Dxx(T(t,x)) ]
Unfortunately errors still arise:
ERROR: MethodError: no method matching collect_ivs_from_nested_operator!(::Set{Any}, ::SymbolicUtils.Add{Real, Float64, Dict{Any, Number}, Nothing}, ::Type{Differential})
Closest candidates are:
collect_ivs_from_nested_operator!(::Any, ::Term, ::Any) at C:\Users\JentschR\.julia\packages\ModelingToolkit\f218S\src\utils.jl:178
Stacktrace:
[1] collect_ivs_from_nested_operator!(ivs::Set{Any}, x::Term{Real, Nothing}, target_op::Type)
@ ModelingToolkit C:\Users\JentschR\.julia\packages\ModelingToolkit\f218S\src\utils.jl:182
[2] collect_ivs(eqs::Vector{Equation}, op::Type)
@ ModelingToolkit C:\Users\JentschR\.julia\packages\ModelingToolkit\f218S\src\utils.jl:155
[3] collect_ivs
@ C:\Users\JentschR\.julia\packages\ModelingToolkit\f218S\src\utils.jl:149 [inlined]
[4] check_equations(eqs::Vector{Equation}, iv::Sym{Real, Base.ImmutableDict{DataType, Any}})
@ ModelingToolkit C:\Users\JentschR\.julia\packages\ModelingToolkit\f218S\src\utils.jl:169
[...]
I assume it’s because
DT(h(T(t,x)))
is still in there. And I have to use @register and predefine derivatives? Like in Symbolic Calculations and Building Fast Parallel Functions · Symbolics.jl But really looking for some help or further docs on this issue. Thanks in advance!
heat equation
using ModelingToolkit, MethodOfLines
using ModelingToolkit: Interval
using Plots
# Variables, parameters, and derivatives
@parameters t x
# Usage of auxilary variables
#@variables T(..) q(..) e(..)
@variables T(..)
Dt = Differential(t)
Dx = Differential(x)
DT = Differential(T)
Dxx = Differential(x)^2
# Loading parametrized material properties and auxilary functions
# ρ = ρ(ξ(T),T)
# h = h(ξ(T),T)
# λ = λ(ξ(T),T)
# ξ(T) = parametrized 7-th order smoothstep function https://en.wikipedia.org/wiki/Smoothstep
# H(x) = Stepfunction
include("matprops.jl");
# Domain edges
t_min, t_max = (0.,200.0)
x_min, x_max = (0.,5.e-3)
# Discretization parameters
dx=0.1e-3;
order=2;
# Original equation
#eqs = [
# Dt(ρ(T(t,x))*h(T(t,x))) ~ Dx(λ(T(t,x)*Dx(T(t,x))))
#]
# Usage of auxiliary variables
#eqs = [Dt(e(t,x)) ~ Dx(q(t,x)),
# e(t,x) ~ ρ(T(t,x))*h(T(t,x)),
# q(t,x) ~ λ(T(t,x))*Dx(T(t,x))]
#
# Usage of chain and product rules
eqs = [ DT(h(T(t,x)))*Dt(T(t,x))*ρ(T(t,x)) + DT(ρ(T(t,x)))*Dt(T(t,x))*h(T(t,x)) ~ λ(T(t,x))*Dxx(T(t,x)) ]
Tᵢ = 273.15+30.0
Tⱼ = 273.15+50.0
# Initial and boundary conditions
bcs = [T(t_min,x) ~ Tᵢ,
T(t,x_min) ~ Tᵢ+(Tⱼ-Tᵢ)*H(t-10.0),
Dx(T(t,x_max)) ~ 0.]
# Space and time domains
domains = [t ∈ Interval(t_min, t_max),
x ∈ Interval(x_min, x_max)]
# PDE system
@named pdesys = PDESystem(eqs, bcs, domains, [t,x], [T(t,x)])
# Method of lines discretization
discretization = MOLFiniteDifference([x=>dx], t; approx_order=order)
prob = discretize(pdesys, discretization)
sol = solve(prob,Tsit5(),abstol=1e-8,saveat=1)
material properties
include("smoothstep.jl")
"""Gumbel-CDF-Funktion"""
function gumbel(T;μ,β)
return exp(-exp(-(T-μ)/β))
end;
"""Ableitung der Gumbel-Funktion"""
function dgumbel(T;μ=μ₀,β=β₀)
z = (T-μ)/β
return (1/β)*exp(-(z+exp(-z)))
end;
"""Heaviside (Step) Funktion"""
function H(t)
0.5 * (sign(t) + 1)
end
"""Rampen Funktion"""
function R(t;offset=0.0,slope=1.0)
return slope*t*H(t)+offset
end
# choose between smoothstep or gumbel function
xi = "smooth"
if xi == "gumbel"
println("Gumbel Funktion")
μ₀=42.7+273.15
β₀=6.360467e-1
ξ(T) = gumbel(T;μ=μ₀,β=β₀)
dξ(T) = dgumbel(T;μ=μ₀,β=β₀)
else
println("Smoothstep Funktion")
T₀ = 40+273.15
T₁ = 42+273.5
ξ(T) = smoothstep(T;edge1=T₀,edge2=T₁)
dξ(T) = dsmoothstep(T;edge1=T₀,edge2=T₁);
iξ(T) = ismoothstep(T;edge1=T₀,edge2=T₁);
end
"""Temperaturabhängige Dichte"""
function ρ(T;ξ=ξ)
ρₛ = 0.86
ρₗ(T) = -0.000636906389320033 * T + 0.975879174796585
return ξ(T) * ρₗ(T) + (1-ξ(T)) * ρₛ
end;
"""Temperaturabhängige spezifische Wärmekapazität"""
function c(T;ξ=ξ,dξ=dξ)
cₗ = 2.0e3
cₛ = 2.0e3
Δhᵥ = 200e3
return ξ(T) * cₗ + (1-ξ(T)) * cₛ + dξ(T)*Δhᵥ
end
"""Temperaturabhängige Wärmeleitfähigkeit"""
function λ(T;ξ=ξ)
λₗ(T) = 0.00000252159374600214 * T^2 - 0.00178215201558874 * T + 0.461994765195229
λₛ(T) = -0.00000469854522539654 * T^2 + 0.00207972452683292 * T + 0.00000172718360302859
return ξ(T) * λₗ(T) + (1-ξ(T)) * λₛ(T)
end;
"""Ableitung der Temperaturabhängigen Wärmeleitfähigkeit"""
function dλ(T;ξ=ξ)
λₗ(T) = 0.00000252159374600214 * T^2 - 0.00178215201558874 * T + 0.461994765195229
dλₗ(T) = 2*0.00000252159374600214 * T - 0.00178215201558874
λₛ(T) = -0.00000469854522539654 * T^2 + 0.00207972452683292 * T + 0.00000172718360302859
dλₛ(T) = -2*0.00000469854522539654 * T + 0.00207972452683292
return dξ(T)*λₗ(T) + dλₗ(T)*ξ(T) - dξ(T)*λₛ(T) + (1-ξ(T))*dλₛ(T)
end;
function β(T;ρ₀)
ρₛ = 0.86
ρₗ(T) = -0.000636906389320033 * T + 0.975879174796585
dρₗ = -0.000636906389320033
return (-1.0/ρ₀)*(dξ(T)*ρₗ(T)+dρₗ*ξ(T)-dξ(T)*ρₛ)
end;
"""Temperaturabhängige spezifische Enthalpie cp = const."""
function h(T;ξ=ξ,href=0.,Tref=35+273.15)
cₗ = 2.0e3
cₛ = 2.0e3
Δh = 200e3
return href + iξ(T)*cₗ+((T-Tref)-iξ(T))*cₛ+ξ(T)*Δh
end
"""Temperaturleitfähigkeit"""
α(T) = λ(T)/(ρ(T)*c(T));
smoothstep
# Using IfElse package for ModelingToolkit compatability
using IfElse
# Scale, and clamp x to 0..1 range
"""clamp function"""
function clamp(x,edge1,edge2)
return (x - edge1) / ( edge2 - edge1)
end
"""derivative of clamp function"""
function dclamp(x,edge1,edge2)
return 1.0 / ( edge2 - edge1)
end
"""7-th order smoothstep function"""
function smoothstep(x;edge1=0.,edge2=1.)
xx = clamp(x,edge1,edge2)
return IfElse.ifelse(xx < 0., 0.,
IfElse.ifelse(xx > 1., 1.0, -20.0*xx^7.0+70.0*xx^6.0-84.0*xx^5.0+35.0*xx^4.0))
end
"""derivative of 7-th order smoothstep function"""
function dsmoothstep(x;edge1=0.,edge2=1.)
xx = clamp(x,edge1,edge2)
dx = dclamp(x,edge1,edge2)
return IfElse.ifelse(xx < 0., 0.,
IfElse.ifelse(xx > 1., 0., (-140.0*xx^6.0+420.0*xx^5.0-420.0*xx^4.0+140.0*xx^3.0)*dx))
end
"""integral of 7-th order smoothstep function"""
function ismoothstep(x;edge1=0.,edge2=1.)
xx = clamp(x,edge1,edge2)
de = edge2 - edge1
iFct(c,j,de) = de*j^(c+1)/(c+1)
return IfElse.ifelse(xx < 0., 0.,
IfElse.ifelse(xx > 1., (x-edge2)-20.0*de/8.0+70.0*de/7.0-84.0*de/6.0+35.0*de/5.0,
-20.0*iFct(7.0,xx,de)+70.0*iFct(6.0,xx,de)-84.0*iFct(5.0,xx,de)+35.0*iFct(4.0,xx,de)))
end
