I have solved a PDE successfully in OpenModelica, but I have problems with solving it in Julia. The Julia code seems to work for up to 8 slices (control volumes), but mysteriously collapses for 9 and more slices…
I’m considering a model of the temperature in a heated water tank, possibly distributed in the z direction. Water may flow through the tank from the bottom (advection), there is convective heat input from the ambient, heater at a specific height (“mathematical point”), and come quasi buoyancy effect.
I’ve already solved the model in OpenModelica, and the results make sense – at least up to 120 slices. Now, I’m trying to solve the model in Julia/DifferentialEquations. In the code (below), I can vary the number of discretization slices from 1 to whatever. The results from OpenModelica and Julia seem to be identical for N=5 slices. For OpenModelica:
and for Julia (sorry about different colors in lines, as well as different number of curves):
However, with N=10 slices, the Julia model collapses. First, OpenModelica results:
followed by the collapsed result from Julia:
OpenModelica gives sensible results up to at least N=120 (i.e., that is how far I have tested it). Julia appears to give ok results up to N=8, but collapses at N=9.
I have probably done a mistake in my code…??? If anyone sees the problem, I appreciate suggestions. Also, my code is probably not super efficient wrt. memory allocation and speed. The relevant code is included below.
# -- Importing packages
using Plots; pyplot()
using LaTeXStrings
using DifferentialEquations
# -- Defining constants
LW1 = 2.5
LW2 = 1.5
LW3 = 1.0
LS1 = :solid
LS2 = :dot
LS3 = :dashdot;
LA1 = 1.0;
LA2 = 0.6;
LA3 = 0.2;
#
# -- model function
function watertank(dT,T,p,t,u)
# Parameters
N = length(T)
# -- mass flow through tank
md = p.rho*u.VdL(t)*u.uv(t)
#
if N == 1
As = 2*p.A + p.per*p.h
V = p.A*p.h
m = p.rho*V
Cp = m*p.chp
Qd_el = p.Po*u.uP(t)
Qd_a = p.Ux*As*(u.Ta(t)-T[1])
dT[1] = (md*p.chp*(u.Ti(t)-T[1]) + Qd_el + Qd_a)/Cp
else
# N > 1
dz = p.h/N
V = p.A*dz
m = p.rho*V
Cp = m*p.chp
#
# -- heater
Qd_el = zeros(N)
#
jP = Int(ceil(p.zP/dz))
Qd_el[jP] = p.Po*u.uP(t)
#
# -- heat from ambient
Qd_a = zeros(N)
#
Qd_a[1] = p.Ux*(p.per*dz+p.A)*(u.Ta(t)-T[1])
Qd_a[2:N-1] = p.Ux*p.per*dz*(u.Ta(t).-T[2:N-1])
Qd_a[N] = p.Ux*(p.per*dz+p.A)*(u.Ta(t)-T[N])
#
# -- heat flow by diffusion
kb = zeros(N-1)
for j in 1:(N-1)
kb[j] = p.rho*p.chp*p.kappa^2*p.d^2*sqrt(p.g*p.alpha_p*max(-(T[j+1]-T[j]),0)/dz)
end
#
Qd_d = zeros(N)
#
Qd_d[1] = p.A*(p.kt+p.cb*kb[1])*(T[2]-T[1])/dz
Qd_d[2:N-1] = p.A*(p.kt.+p.cb*kb[2:N-1]).*(T[3:N]-T[2:N-1])/dz
-p.A*(p.kt.+p.cb*kb[1:N-2]).*(T[2:N-1]-T[1:N-2])/dz
Qd_d[N] = -p.A*(p.kt+p.cb*kb[N-1])*(T[N]-T[N-1])/dz
#
# -- heat total
Qd = Qd_el + Qd_a + Qd_d
#
# -- temperature equations
dT[1] = (md*p.chp*(u.Ti(t)-T[1]) + Qd[1])/Cp
dT[2:N] = (md*p.chp*(T[1:N-1]-T[2:N]) + Qd[2:N])/Cp
end
#
end
#
# -- function for computing output
function Ti_loop(Te,Ti,uv)
return uv*Te + (1-uv)*Ti
end
;
# Model parameters
# -- constants
g = 9.81y
kappa = 0.41
# -- geometry
h = 1.5
d = 0.5
per = pi*d
A = pi*d^2/4
# -- equipment location
zP = 1.15
z1s = 1.3
z2s = 0.23
# -- thermodynamics
rho = 1.0e3
chp = 4.19e3y
alpha_p = 303e-6
# -- heat transport
kt = 0.6
ha = 3
hw = 50
ki_di = 0.5
Ux = 1/(1/ha + 1/hw + 1/ki_di)
# -- heat source
Po = 15e3
#
p = (g=g,kappa=kappa,h=h,d=d,per=per,A=A,zP=zP,rho=rho,chp=chp,alpha_p=alpha_p,kt=kt,cb=cb,Ux=Ux,Po=Po);
#
# Input functions
#
t_stop = 3600*4.;
#
uP_i = t-> t < t_stop/5 ? 1.0 : 0
uv_i = t-> begin
if t < t_stop/6
return 0.75
elseif t > t_stop/2
return 0.25
else
return 0
end
end
VdL_i = t-> t < t_stop/8 ? 8e-3/60 : 4e-3/60
#
Ta_i = t-> 25.
Ti_i = t-> t < t_stop/3 ? 28. : 20.
#
u = (uP=uP_i,uv=uv_i,Ti=Ti_i,Ta=Ta_i,VdL=VdL_i);
# -- creating function into form suitable for ODEProblem
wt_i = (dx,x,p,t) -> watertank(dx,x,p,t,u);
#
# Setting up Simulations
#
# -- time span
tspan = (0.0,t_stop)
#
# --number of slices
nx = 8
#
# -- initial value of states
T0 = fill(25.,nx)
#
# -- defining and solving problem
prob = ODEProblem(wt_i,T0,tspan,p)
sol = solve(prob,Tsit5(),saveat=60)
#
# -- post processing
Sol = reduce(hcat,sol.u)
Tm_ii_L = u.Ti.(sol.t).*(1.0.-u.uv.(sol.t)) .+ Sol[nx,:].*u.uv.(sol.t)
#
dz = h/nx # slice thickness
j1s = Int(ceil(z1s/dz))
j2s = Int(ceil(z2s/dz))
jP = Int(ceil(p.zP/dz))
#
# PLOTTING RESULTS
#
col = range(RGB(0.6,0.6,1),stop=RGB(1,0.7,0.7),length=nx)
#
plot()
for j in 1:nx
plot!(sol.t/3600,Sol[j,:],lw=LW2,lc=col[j],la=0.6,label="")
end
# Model temperature range
plot!(sol.t[1:5:end]/3600,Sol[nx,1:5:end],st=:scatter,m=:x,mc=col[nx],msc=col[nx],ms=4,label=L"\hat{T}_\mathrm{e}")
plot!(sol.t[1:5:end]/3600,Sol[1,1:5:end],st=:scatter,m=:o,mc=col[1],msc=:blue,ms=4,label=L"\hat{T}_1")
plot!(sol.t[1:5:end]/3600,Sol[jP,1:5:end],st=:scatter,m=:v,mc=col[jP],msc=col[jP],ms=4,label=L"\hat{T}(z_\mathrm{P})")
# Sensor temperatures
plot!(sol.t/3600,Sol[j1s,:],lw=LW0,ls=LS2,lc=:maroon,label=L"\hat{T}_1^\mathrm{s}")
plot!(sol.t/3600,Sol[j2s,:],lw=LW0,ls=LS2,lc=:orange,label=L"\hat{T}_2^\mathrm{s}")
plot!(sol.t/3600,u.Ti.(sol.t),lc=:blue,ls=LS1,lw=LW1,la=0.4,label=L"T_\mathrm{i}")
plot!(sol.t/3600,Tm_ii_L,lw=LW0,lc=:coral,ls=LS2,label=L"\hat{T}_\mathrm{i}^\ell")
#
plot!(ylabel=L"temperature $T$ [$^\circ$C]")
plot!(xlim=tspan./3600,legend=:topright,box=true)
