Strange temperature output for 2D heat equation

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

[Harald]

Hi,

could you be more specific about how the results differ from what you expect to see?

Best wishes,

Harald

Answer 1

[Torsten]

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

[Steffen B.]

Thanks for the helpfull answer. I totally forgot this question, my daugther is keeping me busy since she's able to walk.