Solution of 2D Conduction equation of a wall

Objective: To obtain solution of 2D heat conduction equation using discretization and 3 different iterative methods

Method

Steady state method

clear all
close all
clc
% Defining lenghts of domain
Lx=1;
Ly=1;
% Defining the number of mesh points and w
w=1.56;
nx=ny=10;
% Dividing the lengths along x and y directions
x=linspace(0,Lx,nx);
dx=x(2)-x(1);
y=linspace(0,Ly,ny);
dy=y(2)-y(1);
% Defining tolerance and error
error=1e-3;
tol=1e-4;
T=ones(nx,ny);
% Defining BCs and a copy of T
Told=T;
T(1,:)= 600;
T(end,:)= 900;
T(:,1)= 400;
T(:,end)= 800;
iterative_solver = 3;
%Jacobi method
if iterative_solver == 1
    jacobi_iter = 1;
   while(error>tol)
     for i=2:nx-1
     for j=2:ny-1
     T(i,j)= 0.25*(Told(i-1,j)+Told(i+1,j)+Told(i,j-1)+Told(i,j+1));
     end
     end
     contourf(T);
     title_text=sprintf(\'iteration number=%d\',jacobi_iter);
     title(title_text);
     pause(0.003);
     error=max(max((abs(Told-T))));
     Told=T;
     jacobi_iter = jacobi_iter + 1;
   end
end
%Gauss Sidel method
if iterative_solver ==2
   Gaus_Sidel_iter=1;
   while (error>tol)
    for i=2:nx-1
    for j=2:ny-1
    T(i,j)= 0.25*(T(i-1,j)+Told(i+1,j)+T(i,j-1)+Told(i,j+1));
    end
    end
    contourf(T)
    title_text=sprintf(\'iteration number=%d\',Gaus_Sidel_iter);
    title(title_text);
    pause(0.003);
    error=max(max((abs(Told-T))));
    Told=T;
    Gaus_Sidel_iter=Gaus_Sidel_iter+1;
    end
end 
% SOR method
if iterative_solver ==3
   SOR_iter=1;
   while (error>tol)
    for i=2:nx-1
    for j=2:ny-1
    T(i,j)= (1-w)*Told(i,j)+ w*(0.25*(T(i-1,j)+Told(i+1,j)+T(i,j-1)+Told(i,j+1)));
    end
    end
    contourf(T)
    title_text=sprintf(\'iteration number=%d\',SOR_iter);
    title(title_text);
    pause(0.003);
    error=max(max((abs(Told-T))));
    Told=T;
    SOR_iter=SOR_iter+1;
    end
end 

Transient state implicit

clear all
close all
clc
% Defining lenghts of domain
Lx=1;
Ly=1;
% Defining the number of mesh points and w
w=1;
nx=ny=10;
% Dividing the lengths along x and y directions
x=linspace(0,Lx,nx);
dx=x(2)-x(1);
y=linspace(0,Ly,ny);
dy=y(2)-y(1);
% Defining tolerance and error
error=1e-3;
tol=1e-4;
%Defining constants and properties
alpha=0.004041;                         % thermal coeeficient value for coppet%
dt=1e-1;                                %time step value%
k1= alpha*dt/(dx*dx);
k2= alpha*dt/(dy*dy);
T=ones(nx,ny);
% Defining BCs and a copy of T
Told=T;
T_prev_dt=Told;
T(1,:)=Ttop= 600;
T(end,:)=Tbottom= 900;
T(:,1)= Tleft=400;
T(:,end)=Tright= 800;
iterative_solver = 1;
%Jacobi method
if iterative_solver == 1
    jacobi_iter = 1;
   for k=1:200                          %running the interation for 200 time steps
   while(error>tol)                      
     for i=2:nx-1
     for j=2:ny-1
        term1 = (1+(2*k1)+(2*k2))^-1;
        term2 = k1*term1;
        term3 = k2*term1;
        H = (Told(i-1,j)+Told(i+1,j));
        V = (Told(i,j-1)+Told(i,j+1));
        T(i,j) = (T_prev_dt(i,j)*term1)+(H*term2)+(V*term3);
     end
     end
     error=max(max((abs(Told-T))));
     Told=T;
     jacobi_iter = jacobi_iter + 1;
   end
   T_prev_dt = T;
   end
  contourf(x,y,T,\'showText\',\'on\');
  title_text=sprintf(\'iteration number=%d\',jacobi_iter);
  title(title_text);
  pause(0.003);
   end    
 
%Gauss Sidel method
if iterative_solver ==2
   Gaus_Sidel_iter=1;
   for k=1:200
   while (error>tol)
    for i=2:nx-1
    for j=2:ny-1
        term1 = (1+(2*k1)+(2*k2))^-1;
        term2 = k1*term1;
        term3 = k2*term1;
        H = (T(i-1,j)+Told(i+1,j));
        V = (T(i,j-1)+Told(i,j+1));
        T(i,j) = (T_prev_dt(i,j)*term1)+(H*term2)+(V*term3);
    end
    end
    error=max(max((abs(Told-T))));
    Told=T;
    Gaus_Sidel_iter=Gaus_Sidel_iter+1;
    end
    T_prev_dt = T;
    end
    contourf(x,y,T,\'showText\',\'on\');
    title_text=sprintf(\'iteration number=%d\',Gaus_Sidel_iter);
    title(title_text);
    pause(0.003);
    end 
% SOR method
if iterative_solver ==3
   SOR_iter=1;
   for k=1:200
   while (error>tol)
    for i=2:nx-1
    for j=2:ny-1
        term1 = (1+(2*k1)+(2*k2))^-1;
        term2 = k1*term1;
        term3 = k2*term1;
        H = (T(i-1,j)+Told(i+1,j));
        V = (T(i,j-1)+Told(i,j+1));
        T(i,j) = ((1-w)*(T_prev_dt(i,j)*term1))+w*(H*term2)+(V*term3);
    end
    end    
    error=max(max((abs(Told-T))));
    Told=T;
    SOR_iter=SOR_iter+1;
    end
    T_prev_dt = T;
end 
    contourf(x,y,T,\'showText\',\'on\');
    title_text=sprintf(\'iteration number=%d\',SOR_iter);
    title(title_text);
    pause(0.003);
end

Transient state explicit

clear all
close all
clc
% Defining lenghts of domain
Lx=1;
Ly=1;
% Defining the number of mesh points and w
w=1.56;
nx=ny=10;
% Dividing the lengths along x and y directions
x=linspace(0,Lx,nx);
dx=x(2)-x(1);
y=linspace(0,Ly,ny);
dy=y(2)-y(1);
% Defining tolerance and error
error=1e-3;
tol=1e-4;
%Defining constants and properties
alpha=0.004041;                         % thermal coeeficient value for coppet%
dt=1e-1;                                %time step value%
k1= alpha*dt/(dx*dx);
k2= alpha*dt/(dy*dy);
T=ones(nx,ny);
% Defining BCs and a copy of T
Told=T;
T(1,:)=Ttop= 600;
T(end,:)=Tbottom= 900;
T(:,1)= Tleft=400;
T(:,end)=Tright= 800;
iterative_solver = 3;
%Jacobi method
if iterative_solver == 1
    jacobi_iter = 1
   for k=1:200                          %running the interation for 200 time steps
   while(error>tol)                      
     for i=2:nx-1
     for j=2:ny-1
     H=k1*(Told(i+1,j)-2*Told(i,j)+Told(i-1,j));
     V=k2*(Told(i,j+1)-2*Told(i,j)+Told(i,j-1));
     T(i,j)= Told(i,j)+H+V;
     end
     end
     contourf(T);
     title_text=sprintf(\'iteration number=%d\',jacobi_iter);
     title(title_text);
     pause(0.003);
     error=max(max((abs(Told-T))));
     Told=T;
     jacobi_iter = jacobi_iter + 1;
   end
   end
 end
%Gauss Sidel method
if iterative_solver ==2
   Gaus_Sidel_iter=1;
   for k=1:200
   while (error>tol)
    for i=2:nx-1
    for j=2:ny-1
    T(i,j)= 0.25*(T(i-1,j)+Told(i+1,j)+T(i,j-1)+Told(i,j+1));
    end
    end
    contourf(T)
    title_text=sprintf(\'iteration number=%d\',Gaus_Sidel_iter);
    title(title_text);
    pause(0.003);
    error=max(max((abs(Told-T))));
    Told=T;
    Gaus_Sidel_iter=Gaus_Sidel_iter+1;
    end
end
end 
% SOR method
if iterative_solver ==3
   SOR_iter=1;
   for k=1:200
   while (error>tol)
    for i=2:nx-1
    for j=2:ny-1
    T(i,j)= (1-w)*Told(i,j)+ w*(0.25*(T(i-1,j)+Told(i+1,j)+T(i,j-1)+Told(i,j+1)));
    end
    end
    contourf(T)
    title_text=sprintf(\'iteration number=%d\',SOR_iter);
    title(title_text);
    pause(0.003);
    error=max(max((abs(Told-T))));
    Told=T;
    SOR_iter=SOR_iter+1;
    end
end 
end

Results:

A sample result for transient explicit simulation is provided below which shows temperature profile of the wall run using SOR method

The result below shows steady temperature contours of unit square subjected to different temperatures along its walls

Initially when iteration begins, the four sides of wall are at temoeratures that are imposed by their respective boundary conditions. However once the large number of iterations are crossed whose exact number is decided by type of iterative technique that is used to solve the equation, there occurs steady state and temperature contours no longer change further after this number of iterations

For steady state results, SOR requires lesser number of iterations to converge which is controlled by relaxation factor omega(w) followed by Gauss Siedel and then followed by Jacobi method

However for transient state heat transfer with implicit method, where the temperature evolves with respect to time, there is no much of difference between number of iterations required to converge between different iterative methods.