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

2-D steady state is solved first using the three methods. The three methods are formulated in MATLAB as follows.

close all
clear all
clc

% Solving steady state 2D heat conduction equation

% Domain is a square of length equal to one

n=20; % Number of nodes in x and y direction

% Discretization of domain
x=linspace(0,1,n);   % 20 nodes in x-direction
dx=x(2)-x(1);        % grid spacing in x-direction
y=linspace(0,1,n);   % 20 nodes in y-direction
dy=y(2)-y(1);        % grid spacing in y-direction

T_old=ones(n,n)*300;   % Initializing temperature as 300 K everywhere - Initial guess values

T_old(:,1)=400;   % Left boundary is at 400 K
T_old(:,n)=800;   % Right boundary is at 800 K
T_old(1,:)=600;   % Top boundary is at 600 K
T_old(n,:)=900;   % Bottom boundary is at 900 K

% Temperature at the corners is average of the intersecting side temperatures
T_old(1,1)=500;   
T_old(1,n)=700;
T_old(n,1)=650; 
T_old(n,n)=850; 

T_new=T_old;   % T_new is temperature for next iteration

tolerance=0.0001;
error=1;
omega=1.72;    % Relaxation factor for SOR method

counter1=0;    % Counter to count iterations for Jacobi method
counter2=0;    % Counter to count iterations for gauss siedal method
counter3=0;    % Counter to count iterations for SOR method

i=2;           % To select the method for solving

% Jacobi method
if (i==1)
    while(error>tolerance)
    for i=2:n-1
        for j=2:n-1
            T_new(i,j)=(0.5/(dx^2+dy^2))*(((T_old(i-1,j)+T_old(i+1,j))*dy^2)+((T_old(i,j-1)+T_old(i,j+1))*dx^2));
            counter1=counter1+1;
        end
    end
    error=max(max(abs(T_old-T_new)));
    T_old=T_new;
    end
end

% Gauss-siedal
if (i==2)
    while(error>tolerance)
    for i=2:n-1
        for j=2:n-1
            T_new(i,j)=(0.5/(dx^2+dy^2))*(((T_new(i-1,j)+T_old(i+1,j))*dy^2)+((T_new(i,j-1)+T_old(i,j+1))*dx^2));
            counter1=counter1+1;
        end
    end
    error=max(max(abs(T_old-T_new)));
    T_old=T_new;
    end
end

% SOR Method on Gauss-siedal
if (i==3)
    while(error>tolerance)
    for i=2:n-1
        for j=2:n-1
            T_new(i,j)=T_old(i,j)*(1-omega)+omega*((0.5/(dx^2+dy^2))*(((T_new(i-1,j)+T_old(i+1,j))*dy^2)+((T_new(i,j-1)+T_old(i,j+1))*dx^2)));
            counter3=counter3+1;
        end
    end
    error=max(max(abs(T_old-T_new)));
    T_old=T_new;
    end
end     

[X,Y]=meshgrid(x,y);
contourf(X,(1-Y),T_old,'ShowText','on');
xlabel('x');
ylabel('y')
title('2D Steady state temperature distribution using Jacobi');

The colormap of the plot is changed to jet and are as follows. 

1. Plot obtained using Jacobi method

2. Plot obtained using Gauss-siedal method

3. Plot obtained using SOR on Gauss-siedal method

For the SOR method, the value of omega is changed to find the number of iterations it takes to solve. The number of iterations is less for omega=1.72 for this problem.

The number of iterations it took to solve steady-state heat transfer is as follows.

Jacobi method, counter1         =267300

Gauss-siedal, counter2           =142884

SOR on gauss-siedal, counter3 =18792

The MATLAB script to solve the unsteady heat transfer using the explicit method is as follows.

close all
clear all
clc

% Solving Unsteady state 2D heat conduction equation by explicit method

% Domain is a square of length equal to one

n=20; % Number of nodes in x and y direction

% Discretization of domain
x=linspace(0,1,n);   % 20 nodes in x-direction
dx=x(2)-x(1);        % grid spacing in x-direction
y=linspace(0,1,n);   % 20 nodes in y-direction
dy=y(2)-y(1);        % grid spacing in y-direction

alpha=1.5;           % Thermal diffusivity

% For stability (alpha*dt*2/dx^2)<=0.5 so dt<=(dx^2/4/alpha) 

dt=0.0003;           % Time step
nt=2100;             % Number of time steps

T_time=ones(n,n)*300;   % Initial temperature as 300 K everywhere at t=0

% Boundary conditions
T_time(:,1)=400;   % Left boundary is at 400 K
T_time(:,n)=800;   % Right boundary is at 800 K
T_time(1,:)=600;   % Top boundary is at 600 K
T_time(n,:)=900;   % Bottom boundary is at 900 K

% Temperature at the corners is average of the intersecting side temperatures
T_time(1,1)=500;   
T_time(1,n)=700;
T_time(n,1)=650; 
T_time(n,n)=850;

T_new=T_time; % Initializing temperature for next time step

for k=1:nt
    for i=2:n-1
        for j=2:n-1
            T_new(i,j)=T_time(i,j)+(alpha*dt/dx^2)*((T_time(i-1,j)-2*T_time(i,j)+T_time(i+1,j))+(T_time(i,j-1)-2*T_time(i,j)+T_time(i,j+1)));
        end
    end
    T_time=T_new;
end    

The MATLAB script to solve unsteady heat transfer using the three different methods (Implicit) is as follows.

close all
clear all
clc

% Solving Unsteady state 2D heat conduction equation by explicit method

% Domain is a square of length equal to one

n=20; % Number of nodes in x and y direction

% Discretization of domain
x=linspace(0,1,n);   % 20 nodes in x-direction
dx=x(2)-x(1);        % grid spacing in x-direction
y=linspace(0,1,n);   % 20 nodes in y-direction
dy=y(2)-y(1);        % grid spacing in y-direction

alpha=1.5;           % Thermal diffusivity

% For stability (alpha*dt*2/dx^2)<=0.5 so dt<=(dx^2/4/alpha) 

dt=0.0003;           % Time step
nt=2100;             % Number of time steps

T_time=ones(n,n)*300;   % Initial temperature as 300 K everywhere at t=0

% Boundary conditions
T_time(:,1)=400;   % Left boundary is at 400 K
T_time(:,n)=800;   % Right boundary is at 800 K
T_time(1,:)=600;   % Top boundary is at 600 K
T_time(n,:)=900;   % Bottom boundary is at 900 K

% Temperature at the corners is average of the intersecting side temperatures
T_time(1,1)=500;   
T_time(1,n)=700;
T_time(n,1)=650; 
T_time(n,n)=850;

T_old=T_time;     % Initializing temperature for next iteration
T_new=T_time;     % T_new is temperature for next iteration

tolerance=0.0001;
error=1;
omega=1.72;    % Relaxation factor for SOR method

counter1=0;    % Counter to count iterations for Jacobi method
counter2=0;    % Counter to count iterations for gauss siedal method
counter3=0;    % Counter to count iterations for SOR method

i=3;           % To select the method for solving

% Constants used in Implicit formulation
c=(alpha*dt)/dx^2;     % constant value
term1=1/(1+4*c);
term2=c*term1;

% Jacobi method
if (i==1)
   for k=1:nt
    while(error>tolerance)
    for i=2:n-1
        for j=2:n-1
            T_new(i,j)=term1*T_time(i,j)+term2*(T_old(i-1,j)+T_old(i+1,j)+T_old(i,j-1)+T_old(i,j+1));
        end
    end
    error=max(max(abs(T_old-T_new)));
    T_old=T_new;
    end
    T_time=T_new;
    error=1;
   end 
end

% Gauss-siedal
if (i==2)
   for k=1:nt
    while(error>tolerance)
    for i=2:n-1
        for j=2:n-1
            T_new(i,j)=term1*T_time(i,j)+term2*(T_new(i-1,j)+T_old(i+1,j)+T_new(i,j-1)+T_old(i,j+1));
        end
    end
    error=max(max(abs(T_old-T_new)));
    T_old=T_new;
    end
    T_time=T_new;
    error=1;
   end  
end

% SOR Method on Gauss-siedal
if (i==3)
    for k=1:nt
    while(error>tolerance)
    for i=2:n-1
        for j=2:n-1
            T_new(i,j)=T_old(i,j)*(1-omega)+omega*(term1*T_time(i,j)+term2*(T_new(i-1,j)+T_old(i+1,j)+T_new(i,j-1)+T_old(i,j+1)));
        end
    end
    error=max(max(abs(T_old-T_new)));
    T_old=T_new;
    end
    T_time=T_new;
    error=1;
   end  
end

For Unsteady state heat transfer, the value of dt is chosen as 0.0003 based on the stability criteria as follows.

`(alpha*dt)/dx^2+(alpha*dt)/dy^2<=0.5`

The value of dx and dy in this problem is equal. i value is used for selecting the method (Jacobi or gauss-siedal or SOR).