MOLFiniteDifference for 2D transient Heat equation with heat generation

Hi
I am trying to solve the underlying system:

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)]

The Domain is :

domains = [x ∈ Interval(x_min, 0.013),
        y ∈ Interval(0,0.065),
        t ∈ Interval(0, 60*100)]

the system with the parameters:

The initial and Bondary conditions are:

Discretization an solve:

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)

When I solve this problem either scaled or not, The shape of the mean values of my outputs are correct but they have a delay as shown below:

the delay also exists for c(x,y,t).
I have used small N,M and small time steps but the results are the same.
what am I doing wrong in here?
PS: Because there are Neumann and Robin BCs in this problem, The corner points of my problem resulted 0; Does this make sense? Could these 0 corners be the main reason for the delay?

Can you share the full code for recreating that graph with both forms?

You have such widely different parameter values. You can probably rescale t,x and y to reduce it, no?

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:

Can you share the code for how the values from the excel sheet are produced?

So I have 2 Excel files:
one is the mean temperature and the other is the maximum temperature. I used the Maximum file just for the scaling using T_max and t_max. For mean temperature I just plotted the results as below.
The mean Temperature file:

https://docs.google.com/spreadsheets/d/1Tayb-kyOgOtMbRe1cvUWie37lPa3oMNL/edit?usp=sharing&ouid=100996059526622735141&rtpof=true&sd=true

# Load and process data from the Excel file
xc = XLSX.readxlsx("avetemp.xlsx")
sh1 = xc["Sheet1"]
tim = sh1[1,4:end]
tim = vec(tim)
tim = convert(Vector{Float64}, tim)
T_exact = sh1[2,4:end]
T_exact = vec(T_exact)
T_exact = convert(Vector{Float64}, T_exact)
plot(tim,T_exact , linewidth=3, alpha=1, color=:red,label = "Mean  Temperature",
 linestyle=:solid, xlabel="Time (min)",
  ylabel="Temperature (K)", guidefontsize=14,tickfontsize=14, legendfontsize=14,
  guidefont=custom_font, tickfont=custom_font, legendfont=custom_font,framestyle=:box,
  left_margin = [6mm 0mm],
  bottom_margin = [6mm 0mm],
  dpi=300)

Yes but can you share what you’re using as the ground truth for the datagen there?

It’s from COMSOL multiphysics.

The real problem is from the paper below:
https://www.sciencedirect.com/science/article/abs/pii/S037877531830819X

I have validated full model via COMSOL. then For simplicity, I just considered the electrolyte to generate heat.

How do you know it’s not a property of the discretization?

My apologies
I’m not able to understand what you mean.
Is there a problem with MOLFD discretization?
COMSOL uses FEM (shape functions and orders).

Could the cause of the problem be the high exponential term that changes rapidly?
the heat generation caused by Re is highly nonlinear (from 0, it goes up to 4.5e5 in 4 minutes then falls to 0 again). Must I choose Nonelinear solvers?

I think these points are never used in the discretization and cannot be responsible for the delay. But they could be responsible for a small amplitude error when you calculate/diagnose the mean values later on. You should calculate the mean values only on the interior points.

BTW, I notice that nothing in your problem depends on y. Maybe the example provided here is a simplified version of a more complicated setup with an explicit y dependence, otherwise you could work in 1D I suppose?

Thanks for Your reply. Yes it’s simplified. I don’t think the problem is from those zero points because even the maximum temperature has the delay.

Is

match with your T_exact? As the paper’s Fig. 7. A shows the higher heq shifts the curve to left.
BTW, is x dir represents Radial?

The problem is not from parameters. I have checked them out. The problem is the delay. I assume the chosen grids and timesteps make my problem stable and consistent (I used small values and there was no change in the results).

Yes, as it must be.

If it is radial distance, I think the correct diffusion operator should be

kₓ*(Dx(x*Dx(T(x,y,t))))/x

and you get

But if you used exactly the same equation kₓ*Dxx(T(x,y,t)) in your reference solution it still does not explain the delay you observed. If on the other hand your reference solution discretized the correct radial operator, then maybe you have an explanation.

1 Like