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

2-D Steady State and Transient Heat Conduction

AAshwanth Ramesh22 Jun 2019
MATLAB

OBJECTIVE : Solving the 2-D Heat conduction equation using iterative solvers employing different numerical schemes.

APPROACH : The Discretization schemes employed for the respective states of the equation, namely Steady state and Transient state, are :

  • Explicit time integration (and)
  • Implicit time integration

And the types of iterative solvers employed are :

  • Jacobi
  • Gauss-Seidal
  • Successive Over Relaxation (S.O.R)

The first order 2D Heat conduction equation can be mathematically written as

`(delT)/(delt) = alpha[(del^2T)/(delx^2) + (del^2T)/(dely^2)] ---> (1)`

Equation (1) can be of Steady state or Transient in nature

(i) Steady State :

A 2-D Heat conduction of in steady state can be written down as

`[(del^2T)/(delt^2) + (del^2T)/(delt^2)] = 0 --->(2)`

Discretizing the elliptical PDE , there by, implementing central differencing scheme

`(T_(i-1,j)-2T_(i,j)+T_(i+1,j))/(Deltax^2)+(T_(i,j-1)-2T_(i,j)+T_(i,j+1))/(Deltay^2)`

Assuming,

  • `T_(i-1,j) = T_L ; T_(i,j)= T_P ; T_(i+1,j)= T_R;T_(i,j-1) = T_B;T_(i,j+1) = T_T ` indicating the left sided, current, right sided,top sided and bottom sided nodes in the numerical stencil
  • `[2/(Deltax^2)+2/(Deltay^2)] =k`

`(T_L+ T_R)/(Deltax^2) + (T_T+T_B)/(Deltay^2) = k.T_P`

Therefore giving the discretized form of the Steady State 2-D Heat conduction eqn:

`T_P = (1/k)[(T_L+T_R)/(Deltax^2)+(T_T+T_B)/(Deltay^2)]--->(3)`

Eqn (3) Using Jacobi solvers would fetch

`T_P^n = (1/k)[(T_L + T_R)/(Deltax^2)+(T_T+T_B)/(Deltay^2)]`

Eqn (3) Using Gauss-Seidal would fetch

`T_P^n = (1/k)[(T_L^n + T_R)/(Deltax^2)+(T_T+T_B^n)/(Deltay^2)]`

Eqn (3) Using SOR would fetch

`T_P^n = {(1-tau)T_P}+{tau(1/k[(T_L^n + T_R)/(Deltax^2)+(T_T+T_B^n)/(Deltay^2)])}`

(ii) Transient Explicit time integration :

`(delT)/(delt) = alpha[(del^2T)/(delx^2) + (del^2T)/(dely^2)]`

`[(T_(i,j)^(n+1)-T_(i,j)^(n))/(Deltat)]=alpha[(T_(i-1,j)-2T_(i,j)+T_(i+1,j))/(Deltax^2)+(T_(i,j-1)-2T_(i,j)+T_(i,j+1))/(Deltay^2)]^n`

`T_(i,j)^(n+1)=T_(i,j)^(n)+alphaDeltat[(T_(i-1,j)-2T_(i,j)+T_(i+1,j))/(Deltax^2)+(T_(i,j-1)-2T_(i,j)+T_(i,j+1))/(Deltay^2)]^n`

Assuming,

  • `T_(i-1,j) = T_L ; T_(i,j)= T_P ; T_(i+1,j)= T_R;T_(i,j-1) = T_B;T_(i,j+1) = T_T`indicating the left sided, current, right sided,top sided and bottom sided nodes in the numerical stencil
  • `alpha(Deltat)/(Deltax^2)=k_1;alpha(Deltat)/(Deltay^2)=k_2`

`T_P^(n+1) = T_p^n+k_1[T_L -2T_P+ T_R]^n+k_2[T_T -2T_P+ T_B]^n`

`T_P^(n+1) = [1-2k_1-2k_2]T_p^n+k_1[T_L + T_R]^n+k_2[T_T + T_B]^n--->(4)`

Eqn (4) Using Jacobi solver would fetch

`T_P^(n+1) = [1-2k_1-2k_2]T_p^n+k_1[T_L + T_R]^n+k_2[T_T + T_B]^n--->(4)`

Eqn (4) Using Gauss-Seidal solvers would fetch

`T_P^(n+1) = [1-2k_1-2k_2]T_p^n+k_1[T_L^(n+1) + T_R^n]+k_2[T_T^n + T_B^(n+1)]`

Eqn (4) Using SOR solver would fetch

`T_P^(n+1) = {(1-tau)T_P^n}+{tau([1-2k_1-2k_2]T_p^n+k_1[T_L^(n+1) + T_R^n]+k_2[T_T^n + T_B^(n+1)])}`

(iii) Transient Implicit time integration :

`(delT)/(delt) = alpha[(del^2T)/(delx^2) + (del^2T)/(dely^2)]`

`[(T_(i,j)^(n+1)-T_(i,j)^(n))/(Deltat)]=alpha[(T_(i-1,j)-2T_(i,j)+T_(i+1,j))/(Deltax^2)+(T_(i,j-1)-2T_(i,j)+T_(i,j+1))/(Deltay^2)]^(n+1)`

Assuming,

  • `T_(i-1,j) = T_L ; T_(i,j)= T_P ; T_(i+1,j)= T_R;T_(i,j-1) = T_B;T_(i,j+1) = T_T`indicating the left sided, current, right sided,top sided and bottom sided nodes in the numerical stencil
  • `alpha(Deltat)/(Deltax^2)=k_1;alpha(Deltat)/(Deltay^2)=k_2`

`T_P^(n+1) = T_p^n+k_1[T_L -2T_P+ T_R]^(n+1)+k_2[T_T -2T_P+ T_B]^(n+1)`

`T_P^(n+1) = (T_p^n+k_1[T_L + T_R]^(n+1)+k_2[T_T + T_B]^(n+1))/(1-2k_1-2k_2)--->(5)`

Eqn (5) Using Jacobi solver would fetch

`Tprev_P^(n+1) = (T_p^n+k_1[T_L + T_R]^(n+1)+k_2[T_T + T_B]^(n+1))/(1-2k_1-2k_2)`

Eqn (5) Using Gauss-Seidal solver would fetch

`Tprev_P^(n+1) = (T_p^n+k_1[Tprev_L + T_R]^(n+1)+k_2[T_T + Tprev_B]^(n+1))/(1-2k_1-2k_2)`

`Tprev_P^(n+1) ={(1-tau)T_P^n}+ {tau((T_p^n+k_1[Tprev_L + T_R]^(n+1)+k_2[T_T + Tprev_B]^(n+1))/(1-2k_1-2k_2))}`

MATLAB CODE :

(i) Steady State : 

% Steady state
clear all
close all
clc

%%Input Variables
nx = 10;
ny = nx; % 2-D square domain
L = 1;
dx = 1/(nx-1);
dy = dx;
x = linspace(0,1,nx);
y = linspace(0,1,ny);

%Temperature inputs
T0 = 300;        % Initial Guess
T = T0*ones(nx,ny); % Temperature grid
T(:,1)   = 400;     %Left face
T(end,:) = 900;     %Bottom face   %% Boundary conditions1
T(1,:)   = 600;     %Top face
T(:,end) = 800;     %Right face

tou = 2/(1+sin(pi*dx)); % Relaxation factor

%%Stencil numerics
T_old = T;  %Creating a copy of T
w = 2*((dx^2 + dy^2)/(dx^2*dy^2));
tol = 1e-4;
error = 9e9;

%Solver numerics
method = \'Which method would you like to go with? \';
solver = input(method); % Choose 1 for Steady state Jacobi scheme
                        % Choose 2 for Steady state Gauss seidal scheme
                        % Choose 3 for Steady state SOR scheme

%Jacobi scheme
if solver == 1
tic; % Simulation time
iter = 1;  %Calcuation of number of iterations per scheme 
% Spatial marching 
 while(error>tol)            %Convergence loop
     for i=2:nx-1            %Nodal loop 
      for j=2:ny-1        
           T1 = (T_old(i+1,j)+T_old(i-1,j))/(w*dx^2);  %Horizontal component
           T2 = ((T_old(i,j+1)+T_old(i,j-1))/(w*dy^2)); %Vertical component
           T(i,j) = T1 + T2;  %Calculation of T for n+1 timestep
      end
     end
 error = max(max(abs(T - T_old))); %Finding the latest error term in 2D array
 T_old = T;               %Update T for next iteration
 iter = iter+1; 
 end    
 end
 t=toc;
 %Gauss seidal scheme
 if solver == 2
   tic;
   iter = 1;
 while(error>tol)
       for i=2:nx-1
        for j=2:ny-1
        T1 = (T_old(i+1,j)+T(i-1,j))/(w*dx^2);  %Using updated values
        T2 = (T_old(i,j+1)+T(i,j-1))/(w*dy^2);  %for the next iteration of \'i\'
        T(i,j) = T1 + T2;
        end
       end 
 error = max(max(abs(T - T_old)));  
 T_old = T;
 iter = iter+1; 
 end
 end
t=toc;
 %SOR scheme
 if solver == 3
 tic;
 iter = 1;
 while(error>tol)
     for i=2:nx-1
      for j=2:ny-1
       T1 = (T_old(i+1,j)+T(i-1,j))/(w*dx^2);  %Using updated values
       T2 = (T_old(i,j+1)+T(i,j-1))/(w*dy^2);  %for the next iteration of \'i\'
       T(i,j) = tou*(T1 + T2)-((tou-1)*T_old(i,j));
      end
     end 
 error = max(max(abs(T - T_old)));  
 T_old = T;
 iter = iter+1;
 end
t=toc;
 end
% Plot 
solver_used = [\"Steady state Gauss Jacobi\",\"Steady state Gauss Siedal\",\"Steady state SOR\"];
figure (1)
[X,Y] = meshgrid(x,y); %Represent T as a form of X and Y
contourf(X,Y,T,10,\'ShowText\',\'on\'); %Display the height values along the contour lines.
xlabel(\'X-axis\')
ylabel(\'Y-axis\')
title({[(solver_used(solver))],[\'Number of iterations = \',num2str(iter)],[\'Computational time =\',num2str(t),\'secs\']})
colorbar

PLOTS :

%Steady State
Which method would you like to go with? 1
>> error
%Jacobi
error = 9.4868e-05

Which method would you like to go with? 2
>> error
%Gauss-Seidal
error = 9.5865e-05

Which method would you like to go with? 3
>> error
%SOR
error = 8.5742e-05

MATLAB CODE :

(ii) Transient Explicit time integration :

% Explicit transient state
clear all
close all
clc

%%Input Variables
nx = 10;
ny = nx; % 2-D square domain
L = 1;
dx = 1./(nx-1);
dy = dx;
x = linspace(0,1,nx);
y = linspace(0,1,ny);
dt = 0.01;
alpha = 0.3;
tou = 2/(1+sin(pi*dx)); % Relaxation factor

%Temperature inputs
T0 = 300;           %Initial Guess
T = T0*ones(nx,ny); %Temperature grid
T(:,1)   = 400; %Left face
T(end,:) = 600; %Top face   %% Boundary conditions
T(1,:)   = 900; %Bottom face
T(:,end) = 800; %Right face

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

nt = 200; %No of timesteps

%%Stencil numerics
T_old = T;  %Creating a copy of T
tol = 1e-4;
error = 9e9;

%Solver numerics
method = \'Which method would you like to go with? \';
solver = input(method); % Choose 1 for Explicit Jacobi scheme
                        % Choose 2 for Explicit Gauss seidal scheme
                        % Choose 3 for Explicit SOR scheme

%Jacobi scheme
if solver == 1
tic; % Simulation time
iter = 1;  %Calcuation of number of iterations per scheme 
% Spatial marching 
 for t1=1:nt
  while(error>tol)    %Convergence loop
     for i=2:nx-1     %Nodal loop 
      for j=2:ny-1        
           T1 = (1-(2*k1)-(2*k2))*T_old(i,j); 
           T2 = (T_old(i+1,j) + T_old(i-1,j))*k1;  %Horizontal component
           T3 = (T_old(i,j+1) + T_old(i,j-1))*k2;  %Vertical component
          
           T(i,j) = T1 + T2 + T3;  %Calculation of T for n+1 timestep
      end
     end
 error = max(max(abs(T - T_old))); %Finding the latest error term in 2D array
 T_old = T;               %Update T for next iteration
 iter = iter+1; 
 end    
 end
end
 t=toc;
 %Gauss seidal scheme
 if solver == 2
   tic;
   iter = 1;
 for t1=1:nt
  while(error>tol)
       for i=2:nx-1
        for j=2:ny-1
        T1 = (1-(2*k1)-(2*k2))*T_old(i,j);
        T2 = (T_old(i+1,j)+T(i-1,j))*k1;  
        T3 = (T_old(i,j+1)+T(i,j-1))*k2;  
        T(i,j) = T1 + T2 + T3;
        end
       end 
 error = max(max(abs(T - T_old)));  
 T_old = T;
 iter = iter+1; 
 end
 end
 end
t=toc;
 %SOR scheme
 if solver == 3
 tic;
 iter = 1;
 while(error>tol)
     for i=2:nx-1
      for j=2:ny-1
       T1 = (1-tou)*T_old(i,j) + (tou*(1-(2*k1)-(2*k2))*T_old(i,j));
       T2 = (T_old(i+1,j)+T(i-1,j))*k1;  
       T3 = (T_old(i,j+1)+T(i,j-1))*k2;  
      
       T(i,j) = T1+T2+T3;
      end
     end 
 error = max(max(abs(T - T_old)));  
 T_old = T;
 iter = iter+1;
 end
t=toc;
 end
% Plot 
solver_used = [\"Explicit Gauss Jacobi\",\"Explicit Gauss Siedal\",\"Explicit SOR\"];
figure (1)
[X,Y] = meshgrid(x,y); %Represent T as a form of X and Y
contourf(X,Y,T,10,\'ShowText\',\'on\'); %Display the height values along the contour lines.
xlabel(\'X-axis\')
ylabel(\'Y-axis\')
title({[(solver_used(solver))],[\'Number of iterations = \',num2str(iter)],[\'Computational time =\',num2str(t),\'secs\']})
colorbar

PLOTS :

%Explicit
Which method would you like to go with? 1
>> error
%Jacobi
error = 9.8682e-05

Which method would you like to go with? 2
>> error
%Gauss-Seidal
error = 9.0096e-05

Which method would you like to go with? 3
>> error
%SOR
error = 9.3982e-05

MATLAB CODE :

(iii) Transient Implicit time integration :

% Explicit transient state
clear all
close all
clc

%%Input Variables
nx = 10;
ny = nx; % 2-D square domain
L = 1;
dx = 1/(nx-1);
dy = dx;
x = linspace(0,1,nx);
y = linspace(0,1,ny);
dt = 100;
alpha = 0.3;
tou = 2/(1+sin(pi*dx)); % Relaxation factor

%Temperature inputs
T0 = 300;           %Initial Guess
T = T0*ones(nx,ny); %Temperature grid
T(:,1)   = 400;     %Left face
T(end,:) = 900;     %Bottom face   %% Boundary conditions
T(1,:)   = 600;     %Top face
T(:,end) = 800;     %Right face

k1 = (alpha*dt)/(dx^2); 
k2 = (alpha*dt)/(dy^2);
k3 = k1*(1+(2*k1)+(2*k2))^-1;
k4 = k2*(1+(2*k1)+(2*k2))^-1;
nt = 50; %No of timesteps

%%Stencil numerics
T_old = T;  %Creating a copy of T
T_prev_dt = T_old;
tol = 1e-4;
error = 9e9;

%Solver numerics
method = \'Which method would you like to go with? \';
solver = input(method); % Choose 1 for Implicit Jacobi scheme
                        % Choose 2 for Implicit Gauss seidal scheme
                        % Choose 3 for Implicit SOR scheme

%Implicit Jacobi scheme
if solver == 1
tic; %Simulation time
iter = 1;  %Calcuation of number of iterations per scheme 
% Spatial marching 
 for t1=1:nt
  while(error>tol)    %Convergence loop
     for i=2:nx-1     %Nodal loop 
      for j=2:ny-1        
           T1 =  T_prev_dt( i,j)*(1+(2*k1)+(2*k2))^-1;
           T2 = (T_old(i+1,j) + T_old(i-1,j))*k3;     %Horizontal component
           T3 = (T_old(i,j+1) + T_old(i,j-1))*k4;     %Vertical component
           T(i,j) = T1 + T2 + T3;        %Calculation of T for n+1 timestep
      
      end
     end
 iter = iter+1; 
 error = max(max(abs(T - T_old))); %Finding the latest error term in 2D array
  T_old = T;         %Update T for next iteration
 end
 T_prev_dt = T; 
 end
t=toc; 
end
 
 %Implicit Gauss seidal scheme
 if solver == 2
   tic;
   iter = 1;
 while(error>tol)
  for t1=1:nt
       for i=2:nx-1
        for j=2:ny-1
          T1 =  T_prev_dt(i,j)*(1+(2*k1)+(2*k2))^-1;
          T2 = (T_old(i+1,j) + T(i-1,j))*k3;     %Horizontal component
          T3 = (T_old(i,j+1) + T(i,j-1))*k4;     %Vertical component
          T(i,j) = T1 + T2 + T3;
           
        end
      end 
 
 iter = iter+1; 
  error = max(max(abs(T - T_old)));  
  T_old = T;
  end
 T_prev_dt = T; 
 end
 t=toc;
 end

 %Implicit SOR scheme
 if solver == 3
 tic;
 iter = 1; 
 while(error>tol)
  for t1=1:nt  
     for i=2:nx-1  
      for j=2:ny-1
         
        T1 = (T_old(i+1,j)+ (T(i-1,j)))*k3;  
        T2 = (T_old(i,j+1)+ (T(i,j-1)))*k4; 
      
       T(i,j) = T_prev_dt(i,j)*(1+(2*k1)+(2*k2))^-1 + (1-tou)*T_old(i,j) + tou*(T1+T2);
      end
     end 
  error = max(max(abs(T - T_old)));  
  T_old = T;
  iter = iter+1;
  end
  T_prev_dt = T; 
 end
 end
 t=toc;
 % Plot
solver_used = [\"Implicit Gauss Jacobi\",\"Implicit Gauss Siedal\",\"Implicit SOR\"];
figure (1)
[X,Y] = meshgrid(x,y); %Represent T as a form of X and Y
contourf(X,Y,T,10,\'ShowText\',\'on\'); %Display the height values along the contour lines.
xlabel(\'X-axis\')
ylabel(\'Y-axis\')
title({[(solver_used(solver))],[\'Number of iterations = \',num2str(iter)],[\'Computational time =\',num2str(t),\'secs\']})
colorbar

PLOTS :

%Implicit
Which method would you like to go with? 1
>> error
%Jacobi
error = 9.8840e-05

Which method would you like to go with? 2
>> error
%Gauss-Seidal
error = 2.4458e-06

Which method would you like to go with? 3
>> error
%SOR
error = 3.5925e-11

INFERENCE :

  • It can be inferred that the SOR solver for both steady and transient state renders the most time effective iterative solver as compared to the Jacobi and Gauss-Seidal method, due credits to be given to relaxation factor (tau) which employs a linear combination of (quasi) Gauss-Seidal iterative method and the value at the previous timestep.
  • The Explicit SOR solver seems to giving the highest time efficient computation (No. of Iterations =26) , coming close to it is the Steady state SOR solver (No. of Iterations =30). In general, Explicit methods depend greatly on courant number for their stability and in this case a very low CFL number aids in faster convergence.
  • Apart from stability, Accuracy also plays an important role in deciding the convergence and we can see that the Implicit SOR gives the least error upto 11 decimal places inferring that even though implicit schemes takes relatively longer computational time, its much more accurate than rest of the methods.
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.