Thanks for your reply.
tim and T_exact are imported from an Excel file.
Code with dimensions:
using OrdinaryDiffEq, ModelingToolkit, MethodOfLines, DomainSets,Statistics
using XLSX
using Plots
using Plots.PlotMeasures
custom_font = "Times bold"
@parameters x y t
@variables T(..) c(..)
Dt = Differential(t)
Dx = Differential(x)
Dy = Differential(y)
Dxx = Differential(x)^2
Dyy = Differential(y)^2
x_max = 0.013 #m
y_max = 0.065 #m
ρ = 2231.2 #Density
Cₚ = 1100 #[J/(kg*K)]Heat capacity
kₓ = 0.7 #[W/(m*K)] Thermal conductivity in x dir
ky = 140 #[W/(m*K)] Thermal conductivity in y dir
h_conv = 20 #[W/(m^2*K)] heat transfer oefficient
Tᵢ = 20 + 273.15 #[K] Initial temperature
Tₐ = 250 + 273.15 #[K]
αₓ = kₓ/(ρ*Cₚ)
αy = ky/(ρ*Cₚ)
A = 5.14e25 #[1/s] 5.14E25 1/s
E = 2.7e5 #[J/mol] 2.7E5 J/mol
H = 6.2e5 #[J/kg] 6.2E5 J/kg
W = 334.68 #[kg/m^3] 334.68 kg/m³
R = 8.314
t_max = 100*60
Bi = (h_conv*x_max)/kₓ
Re = -Dt(c(x,y,t))
eq = [ρ*Cₚ*Dt(T(x,y,t)) ~ kₓ*Dxx(T(x,y,t)) + ky*Dyy(T(x,y,t)) + H*W*Re ,
Re ~ A*exp(-E/(R*T(x,y,t)))*c(x,y,t)]
domains = [x ∈ Interval(0, x_max),
y ∈ Interval(0,y_max),
t ∈ Interval(0, t_max)]
bcs = [T(x,y,0) ~ Tᵢ,
Dx(T(0,y,t)) ~ 0,
kₓ*Dx(T(x_max,y,t)) ~ h_conv * (Tₐ - T(x_max,y,t)),
Dy(T(x,0,t)) ~ 0,
Dy(T(x,y_max,t)) ~ 0,
c(x,y,0) ~ 1,
Dt(c(x,y,0)) ~ 0,
Dx(c(0,y,t)) ~ 0,
Dx(c(x_max,y,t)) ~ 0,
Dy(c(x,0,t)) ~ 0,
Dy(c(x,y_max,t)) ~ 0
]
@named pdesys = PDESystem(eq, bcs, domains,[x,y,t],[T(x, y,t),c(x,y,t)])
N = 40
M = 20
order = 2
discretization = MOLFiniteDifference([x=>N, y=>M], t, approx_order=order)
@time prob = discretize(pdesys,discretization)
@time sol = solve(prob, AutoTsit5(Rosenbrock23()), saveat = 1)
discrete_x = sol[x]
discrete_y = sol[y]
discrete_t = sol[t]
solT = sol[T(x, y, t)]
solc = sol[c(x, y, t)]
aveT = mean(solT,dims=[1,2])
aveC = mean(solc,dims=[1,2])
using Plots.PlotMeasures
custom_font = "Times bold"
plot(tim,T_exact , linewidth=3, alpha=1, color=:red,
linestyle=:solid, xlabel="Time (min)",
ylabel="Temperature (K)", guidefontsize=14,tickfontsize=14, legendfontsize=14,label="Real",
guidefont=custom_font, tickfont=custom_font, legendfont=custom_font,framestyle=:box,
left_margin = [6mm 0mm],
bottom_margin = [6mm 0mm],
dpi=300)
plot!(discrete_t/60, aveT[1,1,:] , linewidth=3, alpha=1, color=:black,
linestyle=:solid, xlabel="Time (min)",
ylabel="Temperature (K)", guidefontsize=14,tickfontsize=14, legendfontsize=14,label="MOL",
guidefont=custom_font, tickfont=custom_font, legendfont=custom_font,framestyle=:box,
left_margin = [6mm 0mm],
bottom_margin = [6mm 0mm],
dpi=300)
Result:
Code with scaling:
using OrdinaryDiffEq, ModelingToolkit, MethodOfLines, DomainSets,Statistics
using XLSX
using Plots
using Plots.PlotMeasures
custom_font = "Times bold"
@parameters x y t
@variables T(..) c(..)
Dt = Differential(t)
Dx = Differential(x)
Dy = Differential(y)
Dxx = Differential(x)^2
Dyy = Differential(y)^2
x_max = 0.013 # width m
y_max = 0.065 # length m
ρ = 2231.2 #Density
Cₚ = 1100 #[J/(kg*K)]Heat capacity
kₓ = 0.7 #[W/(m*K)] Thermal conductivity in x dir
ky = 140 #[W/(m*K)] Thermal conductivity in y dir
h_conv = 20 #[W/(m^2*K)] heat transfer oefficient
Tᵢ = 20 + 273.15 #[K] Initial temperature
Tₐ = 250 + 273.15 #[K] Ambient temperature
αₓ = kₓ/(ρ*Cₚ)
αy = ky/(ρ*Cₚ)
A = 5.14e25 #[1/s] 5.14E25 1/s
E = 2.7e5 #[J/mol] 2.7E5 J/mol
H = 6.2e5 #[J/kg] 6.2E5 J/kg
W = 334.68 #[kg/m^3] 334.68 kg/m³
R = 8.314
T_max = 370.72 + 273.15 # characteristic Temperature: maximum temperature from real solution
t_max = 27.835*60 # characteristic time: The time that maximum temperature occure
ΔT = T_max - Tᵢ
θₐ = (Tₐ - Tᵢ) /ΔT # scaled ambient temperature
θᵢ = 0 # scaled initial temperature
tf = 100*60/t_max # scaled final time
Re = -Dt(c(x,y,t))
eq = [Dt(T(x,y,t)) ~ (αₓ/(x_max^2))*t_max*Dxx(T(x,y,t)) + (αy/(y_max^2))*t_max*Dyy(T(x,y,t)) + (H*W*Re)/(ρ*Cₚ*ΔT) ,
Re ~ A*t_max*exp(-E/(R*(T(x,y,t)*ΔT+Tᵢ)))*c(x,y,t)]
domains = [x ∈ Interval(0, 1),
y ∈ Interval(0,1),
t ∈ Interval(0, tf)]
bcs = [T(x,y,0) ~ θᵢ,
Dx(T(0,y,t)) ~ 0,
Dx(T(1,y,t)) ~ ((h_conv*x_max)/kₓ) * (θₐ - T(1,y,t)),
Dy(T(x,0,t)) ~ 0,
Dy(T(x,1,t)) ~ 0,
c(x,y,0) ~ 1,
Dt(c(x,y,0)) ~ 0,
Dx(c(0,y,t)) ~ 0,
Dx(c(1,y,t)) ~ 0,
Dy(c(x,0,t)) ~ 0,
Dy(c(x,1,t)) ~ 0
]
@named pdesys = PDESystem(eq, bcs, domains,[x,y,t],[T(x, y,t),c(x,y,t)])
N = 40
M = 20
order = 2
discretization = MOLFiniteDifference([x=>N, y=>M], t, approx_order=order)
@time prob = discretize(pdesys,discretization)
@time sol = solve(prob, AutoTsit5(Rosenbrock23()), saveat = 0.01)
discrete_x = sol[x]
discrete_y = sol[y]
discrete_t = sol[t]
solT = sol[T(x, y, t)]
solc = sol[c(x, y, t)]
aveT = mean(solT,dims=[1,2])
aveC = mean(solc,dims=[1,2])
plot(tim,T_exact , linewidth=3, alpha=1, color=:red,
linestyle=:solid, xlabel="Time (min)",
ylabel="Temperature [K]", guidefontsize=14,tickfontsize=14, legendfontsize=14,label="Real",
guidefont=custom_font, tickfont=custom_font, legendfont=custom_font,framestyle=:box,
left_margin = [6mm 0mm],
bottom_margin = [6mm 0mm],
dpi=300)
temp = aveT[1,1,:] * ΔT .+ Tᵢ
time_ = discrete_t*t_max/60
plot!(time_, temp , linewidth=3, alpha=1, color=:black,
linestyle=:solid, xlabel="Time (min)",
ylabel="Temperature (K)", guidefontsize=14,tickfontsize=14, legendfontsize=14,label="MOL",
guidefont=custom_font, tickfont=custom_font, legendfont=custom_font,framestyle=:box,
left_margin = [6mm 0mm],
bottom_margin = [6mm 0mm],
dpi=300)
Results:
