Strange temperature output for 2D heat equation

3 次查看（过去 30 天）
Steffen B.
Steffen B. 2023-7-19
评论： Steffen B. 2024-1-19
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 个评论
Harald
Harald 2023-7-20
Hi,
could you be more specific about how the results differ from what you expect to see?
Best wishes,
Harald

请先登录，再进行评论。

请先登录，再回答此问题。

采纳的回答

Torsten
Torsten 2023-7-20
编辑：Torsten 2023-7-20
Ran in:
% 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 个评论
Steffen B.
Steffen B. 2024-1-19
Thanks for the helpfull answer. I totally forgot this question, my daugther is keeping me busy since she's able to walk.

请先登录，再进行评论。

更多回答（0 个）

类别

在 Help Center 和 File Exchange 中查找有关 Mathematics and Optimization 的更多信息

标签

产品


版本

R2021b

Community Treasure Hunt

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

Start Hunting!

Translated by