Effacer les filtres
Effacer les filtres

Strange temperature output for 2D heat equation

7 vues (au cours des 30 derniers jours)
Steffen B.
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 commentaire
Harald
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
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 commentaire
Steffen B.
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.

Plus de réponses (0)

Produits


Version

R2021b

Community Treasure Hunt

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

Start Hunting!

Translated by