BlogFor BusinessTalk to a counsellor
CoursesIIT ProgramsPlacementsLog in
BlogFor Business
Talk to a counsellor
Challenge

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

MMohammad Wasi8 Mar 2021

Steady state analysis

clear all
close all
clc

% Steady state Form

nx = 10 ; % Number of grid points along x-axis
ny = nx ; % Number of grid points along y-axis
 
x = linspace(0,1,nx);
y = linspace(0,1,ny);

dx = x(2)-x(1) ;
dy = y(2)-y(1) ;
[xx,yy] = meshgrid(x,y) ;

% Assing boundry conditions
T = ones(nx,ny) ;
T(1,:)=600;
T(end,:)=900;
T(:,1)=400 ;
T(:,end)=800;
T(1,1)=500;
T(10,1)=650;
T(1,10)=700;
T(10,10)=850;
Told = T;
iterative_solver = [1,2,3];
error = 100 ;
tolerence = 1e-5 ;

k = (2*(dx^2+dy^2)/(dx^2*dy^2)) ;
for m = 1:length(iterative_solver)
% Jacobi method
if iterative_solver(m) == 1
    jacobi_iter = 1;
    while(error>tolerence)
        for i = 2:(nx-1)
            for j = 2:(ny-1)
                TR = Told(i+1,j); % Right
                TL = Told(i-1,j); % Left
                TT = Told(i,j+1); % Top
                TB = Told(i,j-1); % Bottom
                H = (TL + TR)/dx^2;
                V = (TT + TB)/dy^2;
              
              T(i,j) = (H + V)/k;
            end
        end
        error = max(max(abs(Told-T)));
        Told = T;
        jacobi_iter = jacobi_iter + 1;
        

         % plot(x,y)
     figure(1)
        [x,y] = contourf(xx,yy,T)
        colormap(jet)
        colorbar
        clabel(x,y)
        title_text = sprintf('iteration number by jacobi = %d',jacobi_iter);
        title(title_text)
        pause(0.03)
    end
end

      
% Gauss_siedel method
 error_GS = 100;
 tolerence = 1e-5 ;
 T = ones(nx,ny) ;
 T(1,:)=600;
T(end,:)=900;
T(:,1)=400 ;
T(:,end)=800;
T(1,1)=500;
T(10,1)=650;
T(1,10)=700;
T(10,10)=850;

       Told = T;
if iterative_solver(m) == 2
    GS_iter = 1;
    while(error>tolerence)
        for i = 2:(nx-1)
            for j = 2:(ny-1)
                TR = Told(i+1,j); % Right
                TL = T(i-1,j);    % Left
                TT = Told(i,j+1); % Top
                TB = T(i,j-1);    % Bottom
                H = (TL + TR)/dx^2;
                V = (TT + TB)/dy^2;
              
              T(i,j) = (H + V)/k;
            end
        end
        error = max(max(abs(Told-T)));
        Told = T;
        GS_iter = GS_iter + 1;
     
    

 % plot(x,y)
     figure(2)
        [x,y] = contourf(xx,yy,T)
        colormap(jet)
        colorbar
        clabel(x,y)
        title_text = sprintf('iteration number by GS = %d',GS_iter);
        title(title_text)
        pause(0.03)
    end
end


 % Successive over relaxation (SOR) method
 error =100;
 tolerence = 1e-5 ;
T = ones(nx,ny) ;
T(1,:)=600;
T(end,:)=900;
T(:,1)=400 ;
T(:,end)=800;
T(1,1)=500;
T(10,1)=650;
T(1,10)=700;
T(10,10)=850;
Told = T;
 omega=1.3
if iterative_solver(m) ==3
    SOR_iter = 1;
    while(error>tolerence)
        for i = 2:(nx-1)
            for j = 2:(ny-1)
                TR = Told(i+1,j); % Right
                TL = T(i-1,j);    % Left
                TT = Told(i,j+1); % Top
                TB = T(i,j-1);    % Bottom
                H = (TL + TR)/dx^2;
                V = (TT + TB)/dy^2;
              
              TGS(i,j) = (H + V)/k;
              T(i,j)= (1-omega)*Told(i,j) + omega*TGS(i,j);
            end
        end
        error = max(max(abs(Told-T)));
        Told = T;
        SOR_iter = SOR_iter + 1;
     
    

 % plot(x,y)
     figure(3)
        [x,y] = contourf(xx,yy,T)
        colormap(jet)
        colorbar
        clabel(x,y)
        title_text = sprintf('iteration number by SOR = %d',SOR_iter);
        title(title_text)
        pause(0.03)
    end
        
end    
end

Transient State Analysis-Implicit 

clear all
close all
clc

% Trancient State

nx = 10 ; % Number of grid points along x-axis
ny = nx ; % Number of grid points along y-axis
n = 10 ;

x = linspace(0,1,nx);
y = linspace(0,1,ny);

a = 1.4 ;
alpha = a ;

dx = x(2)-x(1) ;
dy = y(2)-y(1) ;
[xx,yy] = meshgrid(x,y) ;

% Assing boundry conditions
T = ones(nx,ny) ;
T_prev = T ;
T(1,:)=600;
T(end,:)=900;
T(:,1)=400 ;
T(:,end)=800;
T(1,1)=500;
T(10,1)=650;
T(1,10)=700;
T(10,10)=850;
Told = T;
iterative_solver = 1;
error = 100 ;
tolerence = 1e-5 ;

% Time
t(1) = 0;       % start
t(2) = 15;      % stop
dt = t(2)-t(1);

k1 = a*(dt/dx^2);
k2 = a*(dt/dy^2);

% Jacobi method
if iterative_solver ==1
    jacobi_iter = 1;
    while(error>tolerence)
        for i = 2:(nx-1)
            for j = 2:(ny-1)
                TR = Told(i+1,j); % Right
                TL = Told(i-1,j); % Left
                TT = Told(i,j+1); % Top
                TB = Told(i,j-1); % Bottom
                %TP = dt(i,j)    % Center
                
              term1 = (1+(2*k1)+(2*k2))^-1;
              term2 = k1*term1;
              term3 = k2*term1;
              H = (TR + TL);
              V = (TT + TB);
              
              T(i,j) = (T_prev(i,j)*term1)+(H*term2)+(V*term3);

            end
        end
        error = max(max(abs(Told-T)));
        Told = T;
        jacobi_iter = jacobi_iter + 1;
        
    
    end
end
         % plot(x,y)
     figure(1)
        [x,y] = contourf(xx,yy,T)
        colormap(jet)
        colorbar
        clabel(x,y)
        title_text = sprintf('iteration number by jacobi = %d',jacobi_iter);
        title(title_text)
        pause(0.03)
               
% Gauss_siedel method
T = ones(nx,ny) ;
T_prev = T ;
T(1,:)=600;
T(end,:)=900;
T(:,1)=400 ;
T(:,end)=800;
T(1,1)=500;
T(10,1)=650;
T(1,10)=700;
T(10,10)=850;
Told = T;
error = 100 ;
tolerence = 1e-5 ;

if iterative_solver ==1
    GS_iter = 1;
    while(error>tolerence)
        for i = 2:(nx-1)
            for j = 2:(ny-1)
                TR = Told(i+1,j); % Right
                TL = T(i-1,j); % Left
                TT = Told(i,j+1); % Top
                TB = T(i,j-1); % Bottom
                %TP = dt(i,j)    % Center
                
              term1 = (1+(2*k1)+(2*k2))^-1;
              term2 = k1*term1;
              term3 = k2*term1;
              H = (TR + TL);
              V = (TT + TB);
              
              T(i,j) = (T_prev(i,j)*term1)+(H*term2)+(V*term3);

            end
        end
        error = max(max(abs(Told-T)));
        Told = T;
        GS_iter = GS_iter + 1;
        
    
    end

         % plot(x,y)
     figure(2)
        [x,y] = contourf(xx,yy,T)
        colormap(jet)
        colorbar
        clabel(x,y)
        title_text = sprintf('iteration number by GS = %d',GS_iter);
        title(title_text)
        pause(0.03)
end

% Successive over relaxation (SOR) method
T = ones(nx,ny) ;
T_prev = T ;
T(1,:)=600;
T(end,:)=900;
T(:,1)=400 ;
T(:,end)=800;
T(1,1)=500;
T(10,1)=650;
T(1,10)=700;
T(10,10)=850;
Told = T;
error = 100 ;
tolerence = 1e-5 ;

omega=1.3
if iterative_solver ==1
    SOR_iter = 1;
    while(error>tolerence)
        for i = 2:(nx-1)
            for j = 2:(ny-1)
                TR = Told(i+1,j); % Right
                TL = T(i-1,j); % Left
                TT = Told(i,j+1); % Top
                TB = T(i,j-1); % Bottom
                %TP = dt(i,j)    % Center
                
              term1 = (1+(2*k1)+(2*k2))^-1;
              term2 = k1*term1;
              term3 = k2*term1;
              H = (TR + TL);
              V = (TT + TB);
              
              TGS(i,j) = (T_prev(i,j)*term1)+(H*term2)+(V*term3);
              T(i,j) = (1-omega)*Told(i,j) + omega*TGS(i,j);

            end
        end
        error = max(max(abs(Told-T)));
        Told = T;
        SOR_iter = SOR_iter + 1;
        
    
    end

         % plot(x,y)
     figure(3)
        [x,y] = contourf(xx,yy,T)
        colormap(jet)
        colorbar
        clabel(x,y)
        title_text = sprintf('iteration number by SOR = %d',SOR_iter);
        title(title_text)
        pause(0.03)
end

Transient State Analysis-Explicit

clear all
close all
clc

% Trancient State

nx = 10 ; % Number of grid points along x-axis
ny = nx ; % Number of grid points along y-axis
n = 10 ;

x = linspace(0,1,nx);
y = linspace(0,1,ny);

a = 1.4 ;
alpha = a ;

dx = x(2)-x(1) ;
dy = y(2)-y(1) ;
[xx,yy] = meshgrid(x,y) ;

% Assing boundry conditions
T = ones(nx,ny) ;
T_prev = T ;
T(1,:)=600;
T(end,:)=900;
T(:,1)=400 ;
T(:,end)=800;
T(1,1)=500;
T(10,1)=650;
T(1,10)=700;
T(10,10)=850;
Told = T;
iterative_solver = 1;
error = 100 ;
tolerence = 1e-5 ;

% Time
% t(1) = 0;       % start
% t(2) = 0.002;      % stop
dt = 0.002;

k1 = a*(dt/dx^2); % CFL along x
k2 = a*(dt/dy^2); % CFL along y

CFL = k1 + k2  ;

%CFL(x) + CFL(y) < 'or' = 0.5
% k1 +k2 < 'or' = 0.5


% Jacobi method
if iterative_solver ==1
    jacobi_iter = 1;
    while(error>tolerence)
        for i = 2:(nx-1)
            for j = 2:(ny-1)
                TR = Told(i+1,j); % Right
                TL = Told(i-1,j); % Left
                TT = Told(i,j+1); % Top
                TB = Told(i,j-1); % Bottom
                TP = Told(i,j);
                  H = k1*((TL)-2*(TP)+(TR));
                  V = k2*((TT)-2*(TP)+(TB));
                  
                  T(i,j) = TP+H+V ;
                  
            end
        end
        error = max(max(abs(Told-T)));
        Told = T;
        jacobi_iter = jacobi_iter + 1;
        
    end  
    

         % plot(x,y)
     figure(1)
        [x,y] = contourf(xx,yy,T)
        colormap(jet)
        colorbar
        clabel(x,y)
        title_text = sprintf('iteration number by jacobi = %d',jacobi_iter);
        title(title_text)
        pause(0.03)
    
end

 % Gauss_siedel method
T = ones(nx,ny) ;
T_prev = T ;
T(1,:)=600;
T(end,:)=900;
T(:,1)=400 ;
T(:,end)=800;
T(1,1)=500;
T(10,1)=650;
T(1,10)=700;
T(10,10)=850;
Told = T;
error = 100 ;
tolerence = 1e-5 ;

if iterative_solver ==1
    GS_iter = 1;
    while(error>tolerence)
        for i = 2:(nx-1)
            for j = 2:(ny-1)
                TR = Told(i+1,j); % Right
                TL = T(i-1,j); % Left
                TT = Told(i,j+1); % Top
                TB = T(i,j-1); % Bottom
                TP = Told(i,j);
                  H = k1*((TL)-2*(TP)+(TR));
                  V = k2*((TT)-2*(TP)+(TB));
                  
                  T(i,j) = TP+H+V ;
                  
            end
        end
        error = max(max(abs(Told-T)));
        Told = T;
        GS_iter = GS_iter + 1;
        
    
    end

         % plot(x,y)
     figure(2)
        [x,y] = contourf(xx,yy,T)
        colormap(jet)
        colorbar
        clabel(x,y)
        title_text = sprintf('iteration number by GS = %d',GS_iter);
        title(title_text)
        pause(0.03)
end

 % Successive over relaxation (SOR) method
T = ones(nx,ny) ;
T_prev = T ;
T(1,:)=600;
T(end,:)=900;
T(:,1)=400 ;
T(:,end)=800;
T(1,1)=500;
T(10,1)=650;
T(1,10)=700;
T(10,10)=850;
Told = T;
error = 100 ;
tolerence = 1e-5 ;

omega=1.3
if iterative_solver ==1
    SOR_iter = 1;
    while(error>tolerence)
        for i = 2:(nx-1)
            for j = 2:(ny-1)
                TR = Told(i+1,j); % Right
                TL = T(i-1,j); % Left
                TT = Told(i,j+1); % Top
                TB = T(i,j-1); % Bottom
                TP = Told(i,j);
                  H = k1*((TL)-2*(TP)+(TR));
                  V = k2*((TT)-2*(TP)+(TB));
                  
                  TGS(i,j) = TP+H+V ;
                  T(i,j) = (1-omega)*TP + omega*TGS(i,j);
                  
                  
            end
        end
        error = max(max(abs(Told-T)));
        Told = T;
        SOR_iter = SOR_iter + 1;
        
    
    end

         % plot(x,y)
     figure(3)
        [x,y] = contourf(xx,yy,T)
        colormap(jet)
        colorbar
        clabel(x,y)
        title_text = sprintf('iteration number by SOR = %d',SOR_iter);
        title(title_text)
        pause(0.03)
end

Explanation :-

For the 2D heat conduction equation by using the point iterative techniques we can conclude that Jacobi method takes more iteration numbers than Gauss-seidel and Successive over-relaxation methods. Here Successive over-relaxation (SOR) method is more efficent as it takes less number of iteration numbers which also leads to less time in solving as it includes Thermal diffusivity. GS and SOR are almost the same just that SOR includes Thermal diffusivity.

Steady state analysis is simpler than Transient State Analysis. Transient State includes Time Deravative which has to be maintained accordingly for CFL value. Like Explicit scheme Time Derivative should be maintained less than or equal to CFL value.

All the iterative values for Steady state analysis and Transient State Analysis for given Jacobi, Gauss-seidel, and Successive over-relaxation are almost the same.

Your portfolio, next

Ship work like this
in your own course.

Talk to a career counsellor

Free 20-minute call. We'll map the right course to where you want to work.

By submitting you agree to our Terms and Privacy Policy.

Talk to a career counsellor

Free 20-minute call. We'll map the right course to where you want to work.

By submitting you agree to our Terms and Privacy Policy.