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

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

AAkshaykumar Sonawane22 Jun 2020
  1. AIM:
  • To write program in MATLAB solve the ‘Steady State and Transient State analysis 2D Heat Conduction Equation’.
  • Solve the 2D heat conduction equation by using the point iterative techniques.
  1. OBJECTIVE:
  • Solve the 2D heat conduction equation implementing the ‘Jacobi’, ’Gauss Seidel’, & ’Successive over relaxation’
  • To write code in both ‘implicit’ and ‘explicit’ schemes.
  • Plot the graphs for both ‘implicit’ and ‘explicit’ schemes for all methods.
  1. BODY OF THE CONTENT:
  • Solving the ‘Steady State and Transient State analysis 2D Heat Conduction Equation’ for programming prospective.

Steady state

Governing Equation:

`((d^2 T)/(dx^2 ))+((d^2 T)/(dy^2 ))=0 ………. Equation 1`

  • use the Central difference scheme for the second order derivative to solve the equation.

`(d^2 T)/(dx^2 )=(T_((i-1,j))-2T_((i,j) )+T_((i+1,j)))/(dx^2 ) ………. Equation 2`

`(d^2 T)/(dy^2 )=(T_((i,j-1))-2T_((i,j) )+T_((i,j+1)))/(dy^2 ) ………. Equation 3`

put the equation 2 & 3 in 1

`(T_((i-1,j))-2T_((i,j) )+T_((i+1,j)))/(dx^2 )+ (T_((i,j-1))-2T_((i,j) )+T_((i,j+1)))/(dy^2 )=0`

  • After rearrangement of above equation, we get the following equation.

`T_((i,j) )= 1/k [ (T_((i-1,j) )+T_((i+1,j) ))/(dx^2 )+ (T_((i,j-1) )+T_((i,j+1) ))/(dy^2 ) ]`

`Therefore, k=(2(dx^2+dy^2))/(dx^2 dy^2 )`

  • Solved the equation of steady state analysis.
  • Following is the program for steady state analysis in explicit form.

Steady state program using explicit method:

Program Syntax:

% Steady state using explicit
clear all 
close all
clc

% step 1

% inputs
nx=10
ny=10

x=linspace(0,10,nx);
y=linspace(0,10,ny);
dx=x(2)-x(1);
dy=y(2)-y(1);

error= 8e5;
tol = 1e-4;
k = 2*((dx^2+dy^2)/(dx^2*dy^2))
omega=1.2;

% step 2

%BC

T_L=400
T_R=800
T_T=600
T_B=900
T=300*ones(nx,ny);

T(:,1)  =T_L
T(:,nx) =T_R
T(1,:)  =T_B
T(ny,:) =T_T


T(1,1)  =((T_L+T_B)/2)
T(nx,1) =((T_L+T_T)/2)
T(nx,ny)=((T_R+T_T)/2)
T(1,ny) =((T_R+T_B)/2)


[xx,yy]=meshgrid(x,y);
Told=T;

% step 3

iterative_solver=3


% jacobi method
if iterative_solver==1
  jacobi_iter=1
while (error>tol)
    % space loop
    for i=2:nx-1
        for j=2:ny-1
            
         T(i,j)=(1/k)*(((Told(i-1,j)+Told(i+1,j))/dx^2)+((Told(i,j-1)+Told(i,j+1))/dy^2));
         
        end   
    end 
        error=max(max(abs(Told-T)))
        Told=T
        jacobi_iter = jacobi_iter+1;
        
end   
        contourf(x,y,T)
        xlabel('x axis')
        ylabel('y axis')
        title_text=sprintf('iteration number = %d',jacobi_iter)
        title(title_text)
       
end   

     
% Gauss seidel
if iterative_solver==2
  guass_seidel_iter=1
while (error>tol)
    % space loop
    for i=2:nx-1
        for j=2:ny-1
            
         T(i,j)=(1/k)*(((T(i-1,j)+Told(i+1,j))/dx^2)+((T(i,j-1)+Told(i,j+1))/dy^2));
         
        end   
    end 
        error=max(max(abs(Told-T)))
        Told=T
        guass_seidel_iter = guass_seidel_iter+1;
         contourf(x,y,T)
        xlabel('x axis')
        ylabel('y axis')
        title_text=sprintf('iteration number = %d',guass_seidel_iter)
        title(title_text)
end   
       
    
end   


% S.O.R
if iterative_solver==3
  S_O_R=1
while (error>tol)
    % space loop
    for i=2:nx-1
        for j=2:ny-1
            % using gauss seidel method to calculate the S.O.R method
         T(i,j)=Told(i,j)*(1-omega)+(omega)*((1/k)*(((T(i-1,j)+Told(i+1,j))/dx^2)+((T(i,j-1)+Told(i,j+1))/dy^2)));
         
        end   
    end 
        error=max(max(abs(Told-T)))
        Told=T
         S_O_R = S_O_R+1;
        
end   
        contourf(x,y,T)
        xlabel('x axis')
        ylabel('y axis')
        title_text=sprintf('iteration number = %d', S_O_R)
        title(title_text)
   
end

Syntax Explanation:

  • We want to solve the 2-D Heat Conduction Equation of steady state analysis.
  • Program of steady state in explicit form divided into three steps.
  • In step 1, To give the input parameter for the program to solve equation. In that we define number nodes for ‘x’ and ‘y’ direction after that divide the nodes equally using linspace command.
  • Define the step size for ‘x’ and ‘y’ & give the error, tolerance and omega value.
  • Omega is relaxation factor using in SOR method.
  • In step 2, defines the boundary conditions for equation.
  • In that we define boundary condition for temperatures for all side and calculating the average temperatures for corners.
  • Use of meshgrid command transform the 1D array into 2D array.
  • Now, in step 3 solving the equation use of different methods.
  • Use of iterative solver we solve the methods to calculate the number of iterations.
  • The screenshots of outputs for the different methods is as follows

Jacobi

Iterations = 208

Gauss Seidel

Iterations = 112

SOR

Iterations = 75

  • When observe above graphs we notice that Jacobi method takes more iterations than Gauss seidel and Gauss seidel method takes more iterations than SOR method to calculate the equation.
  • SOR method has less error other than two methods.

Transient/Unsteady State

Governing Equation:

`dT/dt+∝[(d^2 T)/(dx^2 )+(d^2 T)/(dy^2 )]=0 ………. Equation 1`

  • use the Central difference scheme for the second order derivative to solve the equation.

`dT/dt=((T_(i,j) )^(n+1)-(T_(i,j) )^n)/dt ………. Equation 3`

`(d^2 T)/(dx^2 )=(T_((i-1,j))-2T_((i,j) )+T_((i+1,j)))/(dx^2 ) ………. Equation 3`

`(d^2 T)/(dy^2 )=(T_((i,j-1))-2T_((i,j) )+T_((i,j+1)))/(dy^2 ) ………. Equation 4`

put the equation 2, 3 & 4 in 1

`((T_(i,j) )^(n+1)-(T_(i,j) )^n)/dt +∝[(T_((i-1,j))-2T_((i,j) )+T_((i+1,j)))/(dx^2 )+ (T_((i,j-1))-2T_((i,j) )+T_((i,j+1)))/(dy^2 )]=0`

  • After rearrangement of above equation, we get the following equation.

`(T_(i,j) )^(n+1)=(T_(i,j) )^n+k1(T_((i-1,j) )-2T_((i,j) )+T_((i+1,j) ) )+k2(T_((i,j-1))-2T_((i,j) )+T_((i,j+1)))`

`Therefore, k1=∝dt/(dx^2 ),k2=∝dt/(dy^2 )`

  • Solved the equation of unsteady state analysis.
  • Following is the program for unsteady state analysis

Unsteady state equation program explicit method:

Program Syntax:

% unsteady state equation using explicit method
clear all 
close all
clc

% step 1

% inputs
nx=10
ny=10

x=linspace(0,10,nx);
y=linspace(0,10,ny);
dx=x(2)-x(1);
dy=y(2)-y(1);

error= 8e5;
tol = 1e-4;
omega=1.9;
alpha=1.2;
cfl=0.5;
% CFL = (alpha)*(dt/(dx^2))
dt = ((cfl)*((dx^2))/(alpha))/2
k1 = (alpha)*(dt/(dx^2))
k2 = (alpha)*(dt/(dy^2))

% step 2

%BC
T_L=400
T_R=800
T_T=600
T_B=900
T=300*ones(nx,ny);

T(:,1)  =T_L
T(:,nx) =T_R
T(1,:)  =T_B
T(ny,:) =T_T


T(1,1)  =((T_L+T_B)/2)
T(nx,1) =((T_L+T_T)/2)
T(nx,ny)=((T_R+T_T)/2)
T(1,ny) =((T_R+T_B)/2)


[xx,yy]=meshgrid(x,y);
Told=T;

% step 3

% jacobi method
nt=1000;

  jacobi_iter=1
  
 for k=1:nt;
   
    % space loop
    for i=2:nx-1
        for j=2:ny-1
            
         T(i,j)=((1+(2*(k1))+(2*(k2)))^-1)*((Told(i,j))+((k1)*((Told(i-1,j)+Told(i+1,j)))+((k2)*(Told(i,j+1)+Told(i,j-1)))));
         
        end   
    end 
        error=max(max(abs(Told-T)))
        Told=T
        jacobi_iter = jacobi_iter+1;
 end
        
        contourf(x,y,T)
        xlabel('x axis')
        ylabel('y axis')
        title_text=sprintf('iteration number = %d',jacobi_iter)
        title(title_text)

Syntax Explanation:

  • We want to solve the 2-D Heat Conduction Equation of unsteady state analysis.
  • Program of unsteady state in explicit form divided into three steps.
  • First two steps are same as above as the steady state analysis.
  • In step 3, give the number of time step because in unsteady state analysis we have the time derivative and add the time loop solve the unstable equation.
  • The screenshot of output is as follows

Solved the equation for 1000 iterations.

Unsteady state equation program implicit method:

Program Syntax:

% unsteady state equation using implicit method
clear all 
close all
clc

% step 1

% inputs
nx=10
ny=10

x=linspace(0,10,nx);
y=linspace(0,10,ny);
dx=x(2)-x(1);
dy=y(2)-y(1);

error= 8e5;
tol = 1e-4;
omega=1.9;
alpha=1.2;
cfl=0.5;
% CFL = (alpha)*(dt/(dx^2))
dt = ((cfl)*((dx^2))/(alpha))/2
k1 = (alpha)*(dt/(dx^2))
k2 = (alpha)*(dt/(dy^2))

% step 2

%BC
T_L=400
T_R=800
T_T=600
T_B=900
T=300*ones(nx,ny);

T(:,1)  =T_L
T(:,nx) =T_R
T(1,:)  =T_B
T(ny,:) =T_T


T(1,1)  =((T_L+T_B)/2)
T(nx,1) =((T_L+T_T)/2)
T(nx,ny)=((T_R+T_T)/2)
T(1,ny) =((T_R+T_B)/2)


[xx,yy]=meshgrid(x,y);
Told=T;
T_prev_dt = Told;

% step 3

iterative_solver=3


% jacobi method
nt=100
if iterative_solver==1
  jacobi_iter=1
 for k = 2:nt
    error= 8e5
   while (error>tol)
   
      % space loop
      for i=2:nx-1
          for j=2:ny-1
            
           T(i,j)=((1+(2*(k1))+(2*(k2)))^-1)*((T_prev_dt(i,j))+((k1)*((Told(i-1,j)+Told(i+1,j)))+((k2)*(Told(i,j+1)+Told(i,j-1)))));
         
          end   
      end 
          error=max(max(abs(Told-T)))
          Told=T
          jacobi_iter = jacobi_iter+1;
    end
          T_prev_dt = Told
 end  
        contourf(x,y,T)
        xlabel('x axis')
        ylabel('y axis')
        title_text=sprintf('iteration number = %d',jacobi_iter)
        title(title_text)
end 
        
          
% Gauss seidel
if iterative_solver==2
  guass_seidel_iter=1
 for k = 2:nt
    error= 8e5
   while (error>tol)
   
      % space loop
      for i=2:nx-1
          for j=2:ny-1
            
           T(i,j)=((1+(2*(k1))+(2*(k2)))^-1)*((T_prev_dt(i,j))+((k1)*((T(i-1,j)+T(i+1,j)))+((k2)*(T(i,j+1)+T(i,j-1)))));
         
          end   
      end 
          error=max(max(abs(Told-T)))
          Told=T
          guass_seidel_iter = guass_seidel_iter+1;
    end
          T_prev_dt = Told
 end  
        contourf(x,y,T)
        xlabel('x axis')
        ylabel('y axis')
        title_text=sprintf('iteration number = %d',guass_seidel_iter)
        title(title_text)
end

% S.O.R
if iterative_solver==3
  S_O_R=1
 for k = 2:nt
    error= 8e5
   while (error>tol)
   
      % space loop
      for i=2:nx-1
          for j=2:ny-1
            
           T(i,j)=T_prev_dt(i,j)*(1-omega)+(omega)*(((1+(2*(k1))+(2*(k2)))^-1)*((T_prev_dt(i,j))+((k1)*((T(i-1,j)+T(i+1,j)))+((k2)*(T(i,j+1)+T(i,j-1))))));
         
          end   
      end 
          error=max(max(abs(Told-T)))
          Told=T
          S_O_R = S_O_R+1;
    end
          T_prev_dt = Told
 end  
        contourf(x,y,T)
        xlabel('x axis')
        ylabel('y axis')
        title_text=sprintf('iteration number = %d',S_O_R)
        title(title_text)
end

Syntax Explanation:

  • We want to solve the 2-D Heat Conduction Equation of unsteady state analysis.
  • Program of unsteady state in implicit form divided into three steps.
  • First two steps are same as above as the steady state analysis.
  • We just include CFL number, CFL number indicates the stability of the equation.
  • In step 3, we are starting to solve the unstable equation using iterative solver for the different methods.
  • for that we take 100-time steps to solve the unstable equation implicitly
  • solve the unstable equation using implicit scheme for Jacobi, gauss seidel, SOR methods.
  • After solving these methods in MATLAB, they give output in form of graph is as follows

JACOBI

Iterations = 1355

Gauss Seidel

Iterations = 956

SOR

Iterations = 674

  • When we observed above graphs SOR method takes less time to solve the equation.
  1. Conclusion:
  • We written program in MATLAB to evaluate the numerical solution for 2-D Heat Conduction Equation.
  • Here we used the different methods (Jacobi, Gauss Seidel, SOR) to solve both steady and unsteady equation using both implicit and explicit schemes.
  • In an explicitnumerical method S would be evaluated in terms of known quantities at the previous time step n. An implicit method, in contrast, would evaluate some or all of the terms in S in terms of unknown quantities at the new time step n+1.
  • In an explicit method we solve the equation for the previous or known time step and in implicit method solve the equation for unknown time step.
  • Hence, we concluded that to obtain the solution in less time and error use the SOR method.
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.