# Strange temperature output for 2D heat equation

7 vues (au cours des 30 derniers jours)
Steffen B. le 19 Juil 2023
Commenté : Steffen B. le 19 Jan 2024
Hello,
I'm working on a script to solve the 2D heat equation with the classic BTCS scheme.
A simple quadratic domain with this temperature BCs:
clear all
close all
clc
% Material properties
rho=7850;
Lambda=50;
cp=477;
alpha=Lambda/(rho*cp);
% Geometry and mesh
L = 0.1;
nx = 10;
ny = nx;
dt = 10;
t = 3220;
nt = t/dt;
x = linspace(0,L,nx);
y = linspace(0,L,ny);
dx = x(2)-x(1);
dy = y(2)-y(1);
% CFL requirement
CFL = alpha*dt*((1/dx^2)+(1/dy^2));
% Initial conditions
TL=400;
TR=800;
TT=600;
TB=900;
T = 300*ones(nx,ny);
T(2:nx-1,1) = TB; %Bottom
T(2:nx-1,ny) = TT; %Top
T(1,ny-1) = TL; %Left
T(nx,2:ny-1) = TR; %Rigth
% Cornerpoints
T(1,1) = (T(1,2)+T(2,1))/2;
T(1,ny) = (T(1,ny-1)+T(2,ny))/2;
T(nx,1) = (T(nx-1,1)+T(nx,2))/2;
T(nx,ny) = (T(nx,ny-1)+T(nx-1,ny))/2;
Tprev = T;
Told = T;
% Equation coefficients
k1 = alpha*dt/dx^2;
k2 = alpha*dt/dy^2;
kxy = 1/(1+2*k1+2*k2);
kx = k1*kxy;
ky = k2*kxy;
time=zeros(nt,1);
tol = 1e-7;
iter = 0;
tic
for k = 1:1:nt %Time loop
time(k) = k*dt;
error = 9e9;
while (error>tol)
for j = 2:ny-1 %Y-Space loop
for i = 2:nx-1 %X-Space loop
Hx = Told(i+1,j)+Told(i-1,j);
Hy = Told(i,j+1)+Told(i,j-1);
T(i,j) = kxy*Tprev(i,j)+kx*Hx+ky*Hy;
end
end
error = max(max(abs(Told-T)));
Told = T;
iter = iter+1;
i_BTCS_Jacobi(iter) = iter;
end
Tprev = T;
f = figure(1);
fontSize = 22;
f.Position(3:4) = [1080 720];
contourf(x,y,T')
clabel(contourf(x,y,T'),'FontSize', fontSize)
cb= colorbar;
colormap(jet);
%caxis([400, 900]);
%set(gca,'ydir','reverse')
set(gca,'FontSize',fontSize)
set(cb,'FontSize',fontSize)
xlabel('X -Axis','FontSize', fontSize)
ylabel('Y -Axis','FontSize', fontSize)
title(sprintf('BTCS, Implicit Scheme Jacobi Method t = %.2f',time(k)),'FontSize', fontSize);
title(sprintf('BTCS, Implicit Scheme Jacobi Method i = %.d',iter),'FontSize', fontSize);
pause(0.000000000000000003)
M(iter)=getframe(gcf);
iter = iter+1;
end
time = toc
The results are looking strange:
I don't see whats wrong. Maybe some has a clue
Greetings
Steffen
##### 1 commentaireAfficher -1 commentaires plus anciensMasquer -1 commentaires plus anciens
Harald le 20 Juil 2023
Hi,
could you be more specific about how the results differ from what you expect to see?
Best wishes,
Harald

Connectez-vous pour commenter.

### Réponse acceptée

Torsten le 20 Juil 2023
Modifié(e) : Torsten le 20 Juil 2023
% Material properties
rho=7850;
Lambda=50;
cp=477;
alpha=Lambda/(rho*cp);
% Geometry and mesh
L = 0.1;
nx = 10;
ny = nx;
dt = 1;
t = 10;
nt = t/dt;
x = linspace(0,L,nx);
y = linspace(0,L,ny);
dx = x(2)-x(1);
dy = y(2)-y(1);
% CFL requirement
CFL = alpha*dt*((1/dx^2)+(1/dy^2));
% Initial conditions
TL=400;
TR=800;
TT=600;
TB=900;
T = 300*ones(nx,ny);
T(2:nx-1,1) = TB; %Bottom
T(2:nx-1,ny) = TT; %Top
T(1,2:ny-1) = TL; %Left
T(nx,2:ny-1) = TR; %Rigth
% Cornerpoints
T(1,1) = (T(1,2)+T(2,1))/2;
T(1,ny) = (T(1,ny-1)+T(2,ny))/2;
T(nx,1) = (T(nx-1,1)+T(nx,2))/2;
T(nx,ny) = (T(nx,ny-1)+T(nx-1,ny))/2;
Tprev = T;
Told = T;
% Equation coefficients
k1 = alpha*dt/dx^2;
k2 = alpha*dt/dy^2;
kxy = 1 + 2*k1 + 2*k2;
time=zeros(nt,1);
tol = 1e-7;
for k = 1:1:nt %Time loop
time(k) = k*dt;
error = 9e9;
iter = 0;
while (error>tol)
for j = 2:ny-1 %Y-Space loop
for i = 2:nx-1 %X-Space loop
Hx = Told(i+1,j)+Told(i-1,j);
Hy = Told(i,j+1)+Told(i,j-1);
T(i,j) = (Tprev(i,j) + k1*Hx + k2*Hy)/kxy;
end
end
error = max(max(abs(Told-T)));
Told = T;
iter = iter+1;
end
i_BTCS_Jacobi(k) = iter;
Tprev = T;
end
i_BTCS_Jacobi
i_BTCS_Jacobi = 1×10
16 16 16 16 16 16 16 16 16 16
contourf(x,y,T')
colorbar
colormap(jet)
##### 1 commentaireAfficher -1 commentaires plus anciensMasquer -1 commentaires plus anciens
Steffen B. le 19 Jan 2024
Thanks for the helpfull answer. I totally forgot this question, my daugther is keeping me busy since she's able to walk.

Connectez-vous pour commenter.

### Catégories

En savoir plus sur Heat Transfer dans Help Center et File Exchange

R2021b

### Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by