Week 5.1 - Mid term project - Solving the steady and unsteady 2D heat conduction problem

clear all
close all
clc

%%initial condition
L=1; %length of domain
nx=10;
ny=nx;

x=linspace(0,L,nx);
y=linspace(0,L,ny);
dx=L/(nx-1);%providing the spacing between the grid points along x axis
dy=L/(ny-1);%providing the spacing between the grid points along y axis

error =9e9;
tol=1e-4;

%BC

TL=400;
TR=800;
TT=900;
TB=600;

%Temperature matrix

T=300*ones(nx,ny);
T(1,:)=TB;
T(end,:)=TT;
T(:,1)=TL;
T(:,end)=TR;

%Edge refining

T(1,1)=(600+400)/2;
T(1,end)=(600+800)/2;
T(end,1)=(400+900)/2;
T(end,end)=(800+900)/2;

%Calculation of relaxation factor
%omega=2/(1+sin(pi*dx));

omega=1.2;

%creating a copy of old values
[X,Y]=meshgrid(x,y);
Told=T;
Tinitial = T;

% iterative_solver=1;

% k=2/dx^2 + 2/dy^2;
k = 2*(dx^2 + dy^2)/(dx^2 * dy^2);

% jacobi method
% if iterative_solver ==1
jacobi_iter=1;
while(error>tol)
for i=2:nx-1
for j=2:ny-1
T1=(Told(i-1,j)+Told(i+1,j))/(k*(dx^2));
T2=(Told(i,j-1)+Told(i,j+1))/(k*(dy^2));
T(i,j)=T1+T2;
end
end
error=max(max(abs(Told-T)))
Told=T;
jacobi_iter =jacobi_iter + 1;

end
% end

figure(1)
[A,B]=contourf(X,Y,T);
set(gca,'YDIR','reverse')
colorbar
colormap(jet)
clabel(A,B);
title_text = sprintf('Jacobi iteration= %d',jacobi_iter);
xlabel('X')
ylabel('Y')
title(title_text)
pause(0.03)

%%Gauss seidel method
% if iterative_solver ==2
gs_iter=1;
T = Tinitial;
Told = Tinitial;
error =9e9;
while(error>tol)
for i=2:nx-1
for j=2:ny-1
T1=(T(i-1,j)+Told(i+1,j))/(k*(dx^2));
T2=(T(i,j-1)+Told(i,j+1))/(k*(dy^2));
T(i,j)=T1+T2;
end
end
error=max(max(abs(Told-T)))
Told=T;
gs_iter =gs_iter + 1;

end
% end
figure(2)
[A,B]=contourf(X,Y,T);
set(gca,'YDIR','reverse')
colorbar
colormap(jet)
clabel(A,B);
title_text = sprintf('GS iteration= %d',gs_iter);
xlabel('X')
ylabel('Y')
title(title_text)
pause(0.03)

%%iterative_solver=3;
%SOR method
% if iterative_solver ==3
T = Tinitial;
Told = Tinitial;
error =9e9;
sor_iter=1;
while(error>tol)
for i=2:nx-1
for j=2:ny-1
T1=omega*(T(i-1,j)+Told(i+1,j))/(k*(dx^2));
T2=omega*(T(i,j-1)+Told(i,j+1))/(k*(dy^2));
T(i,j)=(1-omega)*Told(i,j) + T1 + T2;
end
end
error=max(max(abs(Told-T)))
Told=T;
sor_iter =sor_iter + 1;

end
% end

figure(3)
[A,B]=contourf(X,Y,T);
set(gca,'YDIR','reverse')
colorbar
colormap(jet)
clabel(A,B);
title_text = sprintf('SOR iteration= %d',sor_iter);
xlabel('X')
ylabel('Y')
title(title_text)
pause(0.03)