Solving 2-D 2nd Order Linear Transient &Steady  Heat Conduction Equation Using Both Explicit And Implicit Numerical Schemes In MATLAB.

Solving 2-D 2nd Order Linear Transient &Steady  Heat Conduction Equation Using Both Explicit And Implicit Numerical Schemes In MATLAB.

 

  • Explicit Scheme : We say that a scheme is explit when the information at time level "n+1" depends on previous time levels i.e., "n , n-1" e.g. the upwind difference schemes for the convection equation is an explicit schemes.
  • Implicit Scheme : In this schemes the information at time level "n+1" is not only dependent on previous time steps "n , n-1" but also on current time step "n+1" e.g.      Crank-Nicolson method, Euler Implicit methos etc.

 

The 2-D Transient Heat Conduction Equation is given by:

 

The numerical difference equation in implicit form is given as:-

where,

The iterative methods implemented to solve in the implicit way are:

  • Gauss-Jacobi method.
  • Gauss-Seidel method.
  • Successive Over Relaxation method For gauss-Seidel method.

The Matlab Codes are:

% Implicit solver method for 2-D Heat Equation.


clear all
close all
clc


% Grid points in x-direction
nx = 30;

% Grid Points in y-direction
ny = nx;
x = linspace(0,1,nx);
y = linspace(0,1,ny);
% Time Step
dt = 1e-1;

% length of the domain
L_x = 1; % length in x - direction
L_y = 1; % Lengthin y- direction

% Mesh size 
dx = x(2)-x(1);
dy = y(2)-y(1);

% Grid 
T = ones(nx,ny);

% % Initial Temperature over the grid = 100K
T = 300*ones(nx,ny);
% T(1,:) = linspace(100,800,nx);
% T(:,1) = 100;
% 
% for k = 2:nx;
%     
%     for l = 2:ny;
%     
%         T(k,l) = T(k-1,l-1)+100;
%         
%     end
%     
% end
% 
% T_initial_before_boundary_condition = T;


% Applying Boundary condition
% Top_Edge = 600K
% Left_Edge = 400K
% Right_Edge = 800k
% Bottom_Edge = 900K

        
T(1,:) = 900;
T(end,:) = 600;

T(:,1) = 400;
T(:,end) =800;   

T_initial = T;
T_old = T;

iteration  = input('enter input 1.Jacobi 2.Gauss-Seidel 3. SOR Gauss-Seidel')
tic;

% Defining constants
% Thermal Diffusivity of the material (aluminium)

alpha = 1e-5;
K_1 = ((alpha*dt)/(dx^2));
K_2 = ((alpha*dt)/(dy^2));
Term_1 = ((1+(2*K_1)+(2*K_2))^-1);
p = K_1*Term_1;
q = K_2*Term_1;

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


jacob_iter = 1;
% Jacobi Iteration 

if iteration ==1;
 for a = 1:20000
         
    error_1 = error;

    
      while(error_1 > tolerance)
        
          for  i = 2:nx-1;

                for j = 2:ny-1;
           

                      T(i,j) =(T_old(i,j)*Term_1) + ((p*(T_initial(i-1,j)+T_initial(i+1,j)))+(q*(T_initial(i,j-1)+T_initial(i,j+1))));
            

                end
              
          end
          
          error_1 = max(max(abs(T_initial-T)))
          T_initial = T;
          jacob_iter = jacob_iter + 1    
      end
      
    T_old = T;
    
end
          
       
      
toc;
F = toc;
   
  %Plotting the curve
  figure(1)
  [m n] = contourf(x,y,T);
  clabel(m,n);
  colormap(jet);
  colorbar;
  zlabel('x-axis');
  ylabel('y-axis');
  [X,Y] = meshgrid(x,y);
  title(sprintf('Transient State Heat Conduction(Gauss-Jacobi):Total Iteration= %d,Total Time = %fs',jacob_iter,F));
 
 end





   % gauss-Seidel Iteration  
   if iteration ==2;
   Gauss_Seidel = 2; 
   tic; 
   
   
      for b = 1:20000
      
      error = 9e9;
      error_2 = error;
      
        while(error_2 > tolerance)
            
            for k = 2:nx-1;
                
                for l = 2:ny-1;
                    
                    
                    T(k,l) = (T_old(k,l)*Term_1) + (p*(T(k-1,l)+T_initial(k+1,l))+(q*(T(k,l-1)+T_initial(k,l+1))));
               
                end
                
            end
            
            error_2 = max(max(abs(T_initial - T)))
            T_initial = T;
            Gauss_Seidel = Gauss_Seidel+1
                        
        end
       
     T_old = T;
     
end
toc;
F  = toc;
%Plotting the curve
figure(1)
[m n] = contourf(x,y,T);
 clabel(m,n);
 colormap(jet);
 colorbar;
 zlabel('x-axis');
 ylabel('y-axis');
 title(sprintf('Transient State Heat Conduction(Gauss-seidel):Total Iteration= %d,Total Time = %fs',Gauss_Seidel,F));
 end
   
   
  
  % Successive Over Relaxation Gauss-Sidel Iteration  
   SOR_Gauss_Seidel = 3;
   tic;
   if iteration ==3;
   omega = 0.9;
   
      for b = 1:20000
      error = 9e9;
      error_3 = error;
      
        while(error_3 > tolerance)
            
            for k = 2:nx-1;
                
                for l = 2:ny-1;
                    
                    
                    T(k,l) = (T_old(k,l)*(1-omega)) + omega*((T_old(k,l)*Term_1) + (p*(T(k-1,l)+T_initial(k+1,l))+(q*(T(k,l-1)+T_initial(k,l+1)))));
               
                end
                
            end
            
            error_3 = max(max(abs(T_initial - T)))
            T_initial = T;
            SOR_Gauss_Seidel = SOR_Gauss_Seidel+1
                        
        end
       
     T_old = T;
     
end
toc;
F  = toc;
%Plotting the curve
 figure(2)
 [C,h] = contourf(x,y,T);
 clabel(C,h);
 colormap(jet);
 colorbar;
 zlabel('x-axis');
 ylabel('y-axis');
 title(sprintf('Transient State Heat Conduction(SOR Gauss-sidel):Total Iteration= %d n,Total Time = %f s',SOR_Gauss_Seidel,F));
 end

   
   
   
   
   
   
   

Results(Plots):

Observation:

  1. Gauss-jacobi methods takes total iterations of "41720" to converge the results.
  2. Whereas, Gauss-Seidel methods takes total iterations of "40570" which is lesser than gauss-Jacobi iteration method.
  3. However, using SOR method for Gauss-Seidel method we introduce a relaxation parameter "omega = 0.2" for faster convergence rate "omega ∈ (0,2)". the result converges at 40003 iteration which faster than the above iteration methods.
  4. The method is stable  due to implicit approach of the solution.

 

The Explicit form of the 2-D Transient Heat equation:

% Explicit Solver for 2-D Heat Equation.

clear all
close all
clc

% Making grid of square matrix

dt = 1e-3;

nx = 10;

ny = nx;

L = 1;
dx =L/nx-1;
dy = L/ny-1;
x = linspace(0,1,nx);
y = linspace(0,1,ny);

dy = dx;

T = 300* ones(nx,ny);

% T = linspace(200,1000,n)
% 
% for a = 2:10
%     
%     T(a,:) = T(a-1)+100
%     
% end


% Applying Boundary condition
% Top_Edge = 600K
% Left_Edge = 400K
% Right_Edge = 800k
% Bottom_Edge = 900K

        
T(1,:) = 900;
T(end,:) = 600;
T(:,1) = 400;
T(:,end) =800;   

T_initial = T;

alpha = 1e-5;
K_1 = ((alpha*dt)/(dx^2));
K_2 = ((alpha*dt)/(dy^2));

% The stability criteria of the method.
CFL = alpha*dt*((1/dx^2) + (1/dy^2))

% Time loop

Exp_iteration = 1;

tic;
for b = 1:10000
    
   % Starting Time Loop
   Exp_iteration = Exp_iteration;
   
   for i = 2:nx-1
       
       
       for j = 2:ny-1   % NESTED LOOPING
           
           T(i,j) =(T_initial(i,j)) + ((K_1*(T_initial(i-1,j)-(2*T_initial(i,j))+(+T_initial(i+1,j)))+(K_2*(T_initial(i,j-1)-(2*T_initial(i,j))+T_initial(i,j+1)))))
            
       
       end
       
   end
   
   T_initial = T;
   Exp_iteration = Exp_iteration +1
   
end
toc;
F= toc;
figure(2)
  [m n] = contourf(x,y,T);
  clabel(m,n);
  colormap(jet);
  colorbar;
  xlabel('x-axis');
  ylabel('y-axis');
  title('Steady State Heat conduction Explicit Method');
  

Results(plot):

EFFECT OF TIME STEP ON RESULTS FOR EXPLICIT SOLVER:

OBSERVATIONS:

  • Explicit approach is usually not as stable as that of Implicit approach as a result in explicit appraoch we need to use smaller time step in order to maintain the stability whereas we can use higher time step and maintain stability in implicit approach. 
  • The above approach is a time marching of temperature over a grid(plate) via explicit method in which new value of temperature is assigned to a node of time"n+1" from a previous time "n , n-1..". 
  • At time step = 1e+1, we see that the truncation error obtained = 0.0181 and in the plots the temperature ranges are not uniformly distributed where as for time step = 1e-1, the truncation error is =0.0013 and temperature ranges are more closer and more uniformly distributed.Hence, smaller the time step generates smaller truncation error and more accurate results are obtained. 
  • Also if the time step is increased beyond the optimum limit given as CFL (Courent Friedrichs-Lewy Number) then the solution gets unstable.

Steady State Heat Equation Analysis By Implicit Method:

MATLAB CODES:

% Steady state heat equation
clear all
close all
clc

% defining the number of grid points in X,Y direction.
nx = 100;
ny = nx;


% Defining the domain length and mesh size.
L = 1;
dx = L/(nx-1);
dy = L/(ny-1);
x = linspace(0,L,nx);
y = linspace(0,L,ny);
[X Y] = meshgrid(x,y);

% Defining the Temperature distribition over the grid.
T = zeros(nx,ny);
T(1,:) = linspace(100,800,nx);
for a = 2:nx;
    
    for b = 1:ny;
        
        T(a,b) = T(a-1,b);
        
    end
end
    
T_Before_boundary_condition = T; 
T_initial = T;


% Applying Boundary Condition
% Top_Edge = 600K
% Left_Edge = 400K
% Right_Edge = 800k
% Bottom_Edge = 900K

T(1,:)   = 600;
T(end,:) = 900;
T(:,1)   = 400;
T(:,end) = 800;
iteration = input('1. Jacobi_Iteration  2. Gauss-sidel_iteration  3.SOR_Gauss-sidel_iteration');
k = (2*(dx^2 + dy^2)/(dx^2)*(dy^2));

% # Jacobi Iteration for the Steady State Heat Equation.
% # Initializing Convergence Loop
% # Tolerance for steaady state convergence.
Tolerance = 1e-10;

tic;
jacob_iter = 1;
if iteration == 1;
    
    error   = 9e9;
    error_1 = error;
 
    while (error_1 > Tolerance);
        
        
        for i = 2:nx-1;
            
            for j = 2:ny-1;
                
                T(i,j) =((T_initial(i-1,j) + T_initial(i+1,j))/k*(dx^2)) + ((T_initial(i,j-1) + T_initial(i,j+1))/k*(dy^2));
                
            end
            
        end
        
        error_1 = max(max(abs(T_initial - T)))
        T_initial = T;
        jacob_iter = jacob_iter + 1
    
    end
    toc;
    F = toc;
    
 %Plotting the curve
  figure(1)
  [m n] = contourf(x,y,T);
  clabel(m,n);
  colormap(jet);
  colorbar;
  zlabel('x-axis');
  ylabel('y-axis');
  title(sprintf('Steady State heat Conduction(Jacobi):total iteration= %d, Time taken = %fs',jacob_iter,F));
end

  
  
% initializing Gauss-sidel Iteration
tic;
Gauss_Sidel_iter = 2;
if iteration == 2;
    
    error   = 9e9;
    error_2 = error;
    while (error_2 > Tolerance);
        
        for i = 2:nx-1;
            
            for j = 2:ny-1;
                
                T(i,j) = ((T(i-1,j) + T_initial(i+1,j))/k*(dx^2)) + ((T_initial(i,j+1) + T(i,j-1))/k*(dy^2));
                
            end
            
            
        end
        
        error_2 = max(max(abs(T_initial-T)))
        T_initial = T;
        Gauss_Sidel_iter = Gauss_Sidel_iter + 1
        
    end
    toc;
    F = toc;
   %Plotting the curve
  figure(2)
  [m n] = contourf(x,y,T);
  clabel(m,n);
  colormap(jet);
  colorbar;
  zlabel('x-axis');
  ylabel('y-axis');
  title(sprintf('Stady State heat Condition: Total iteration = %d, Time Taken = %fs',Gauss_Sidel_iter,F));
end

  
  
% initializing  Successive Over Relaxation Gauss-sidel Iteration
% Defining Relaxation constant (omega)[0 < omega < 2].
omega = 1.5;
tic;
SOR_Gauss_Sidel_iter = 3;
if iteration ==3;
    
    error   = 9e9;
    error_2 = error;
    while (error_2 > Tolerance);
        
        for i = 2:nx-1;
            
            for j = 2:ny-1;
                
                T(i,j) = (1-omega)*T_initial(i,j) +omega*(((T(i-1,j) + T_initial(i+1,j))/k*(dx^2)) + ((T_initial(i,j+1) + T(i,j-1))/k*(dy^2)));
                
            end
            
            
        end
        
        error_2 = max(max(abs(T_initial-T)))
        T_initial = T;
        SOR_Gauss_Sidel_iter = SOR_Gauss_Sidel_iter + 1
        
    end
    toc;
    F=toc;
   %Plotting the curve
  figure(2)
  [m n] = contourf(x,y,T);
  clabel(m,n);
  colormap(jet);
  colorbar;
  zlabel('x-axis');
  ylabel('y-axis');
  title(sprintf('Steady State Heat conduction:total iteration = %d,Total Time  = %fs',SOR_Gauss_Sidel_iter,F));
 
end
  

 

RESULTS(PLOTS):

OBSERVATIONS:

From the results of steady state heat equation, we can figure out that using Gauss-jacobi method takes more iteration than any other method to determine the temperature distribution. Whereas, using successive over relaxation for Gauss-Seidel method using an optimum relaxation parameter "omega = 0.2" gives a huge margin in the rate of convergence from its predecessor. The final error obtained for following methods are:-

  • Gauss-Jacobi (error)          =  9.9931e-11
  • Gauss-Seidel (error)          =  9.9817e-11
  • SOR Gauss -Seidel (error) =  9.9817e-11

 


Projects by Ankit Chakrabarty

Title:- Geometry cleanup and meshing of an automotive hood for structural analysis   Objective:  Perform geometry cleanup of automotive hood (both inner and outer panel). Perform mid-surface extraction of both the panels and Hem both the surfaces . Perf Read more

Abstract: This project presents the Modelling of the  behaviour of  fluid and its properties when it passes through a convergent-divergent nozzle using a Finite Difference Method. However, to simplify the understanding of the dynamics of the flow in the nozz Read more

 Turbocharger: Literature Review: A turbocharger is mechanical rotating device which forcebly induces air at high pressure into the internal combustion engine, in turn increases the power and effeciency (thermal and volumetric effeciency) of the engine. A turb Read more

  (A) Problem defination:- Solve the following Linear equation using Gauss-Jacobi, Gauss-Seidal and SOR iterative method and understand the following parameters:- Spectral radius and the impact of spectral radius on the convergence of the solution. Effect of d Read more

CASE STUDY:-  A Thermodynamic Analysis Of A Tractor Engine For Use In India With Given Engine Performance Targets. A Brief Literature Review: A tractor is an engineered vehicle specifically designed to deliver Power and traction/tractive effort to haul heavy equi Read more

  CASE STUDY: Effect of Burn Rate on engine performance Of A 6-Cyl-Turbocharged-DI Engine Using A Non-predictive Combustion model (DI WIEBE) Abstract: Analytical functions approximating the burn rate in internal combustion engines are useful and cost-effective to Read more

Solving 1-D Linear Convection Using First-Order Backward Difference And Forward Difference method Using MATLAB LITERATURE REVIEW: WHAT IS CONVECTION? Convection is the sum of bulk transport of the  fluid  and brownian/ osmotic dispersion of fluid constituen Read more

Deriving the 4th order approximation for 2nd order derivative of a function using Taylor table for centered, skewed right-sided stencil and skewed left_sided stencil and compare the error results for the same using MATLAB.   Given function, f(x)    =& Read more

Objective:- (1) Determine the engine parameters running at 1800 rpm and record the same. (2) Further increase the power output obtained at 3600 rpm by 10%.   Model:-  The GT model Of the 1-cylinder engine is shown below:- Observations (@1800rpm):-  Read more


Loading...

The End