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

AIM:

1. Steady state analysis & Transient State Analysis

to solve the 2D heat conduction equation by using the point iterative techniques that were taught in the class. You will implement the following methods

1. Jacobi

2. Gauss-seidel

3. Successive over-relaxation

Your absolute error criteria is 1e-4

2. You will write code for both implicit and explicit schemes (This is for the Transient Analysis only)

Theory:

The hear condution equation is given as

at steady state dT/dt is omitted.

In Jacobi, each iteration is done using a set of previous values, even if new values become available in the middle of an iteration. In Gauss-Seidel, as soon as you have a new iterate for a particular component, you use it for all subsequent computation.In SOR you extrapolate the value by using a omega term over a value obtained by gauss- siedel method to reduce the iteration required for convergence but if omega value is high the system might start diverging and oscillate resulting in more iteration

Explicit method solves the equation to find the temperature value in new step (tn+1) by using known values from present time step(tn)

Implicit method solves the equation to find the temperature value in new step (tn+1) by using unkown values in the new time step(tn+1)

Its preferred to have cfl value less than 0.5 ,cfl=k1+k2

CODE:

2D Heat conduction equation

clear all
close all
clc
nx=10;
ny=nx;
x=linspace(0,1,10);
y=linspace(1,0,10);
dx=x(2)-x(1)
dy=y(2)-y(1)
dt=0.001;
[xx yy]=meshgrid(x,y)
t=400*ones(nx,ny);
t(:,1)=400;
t(:,10)=800;
t(1,:)=600;
t(10,:)=900;
t(1,1)=(t(1,2)+t(2,1))/2;
t(1,10)=(t(1,9)+t(2,10))/2;
t(10,1)=(t(9,1)+t(10,2))/2;
t(10,10)=(t(10,9)+t(9,10))/2;
alpha=1.4;
nt=1400;
k1=alpha*(dt/(dx)^2)
k2=alpha*(dt/(dy)^2)
cfl=k1+k2
told=t;
tprev=t;
ijacobi=1;
igs=1;
isor=1;
ejacobi=1;
egs=1;
esor=1;
sjacobi=1;
sgs=1;
ssor=1;
tolr=1e-4;
error=9e9;
w=1.1;
%Steady state =1 , Transient state=2;
for n=1:2
    state=n;
 %steady state   
if state==1
    for d=1:3
    iteration=d
%steady state jacobi   
if iteration==1
tic
[t_s_j,sjacobi]=steadyjacobi(nx,ny,t,told,tolr,sjacobi,error,k1,k2,dt,dx,dy,nt)
t_s_jinsec=toc
figure(11)
contourf(xx,yy,t_s_j,'ShowText','on')
colormap('jet')
colorbar
xlabel('x-axis')
ylabel('y-axis')
title('steadyjacobi')
end
%steady state gauss-siedel    
if iteration==2
tic
[t_s_gs,sgs]=steadygs(nx,ny,t,told,tolr,sgs,error,k1,k2,dt,dx,dy,nt)
t_s_gsinsec=toc
figure(12)
contourf(xx,yy,t_s_gs,'ShowText','on')
colormap('jet')
colorbar
xlabel('x-axis')
ylabel('y-axis')
title('steadygs')
end
%steady state successive over relaxation method  
if iteration==3
tic
[t_s_sor,ssor]=steadysor(w,nx,ny,t,told,tolr,ssor,error,k1,k2,dt,dx,dy,nt)
t_s_sorinsec=toc
figure(13)
contourf(xx,yy,t_s_sor,'ShowText','on')
colormap('jet')
colorbar
xlabel('x-axis')
ylabel('y-axis')
title('steadysor')
end
 
end
 % comparing time taken for jacobi , gauss-siedel and sor in steady state  
figure(14)
sec=[t_s_jinsec,t_s_gsinsec,t_s_sorinsec]
bar(sec,'r')
title('STEADY TIME')
% comparing iteration taken for jacobi , gauss-siedel and sor in steady state  
figure(15)
iter=[sjacobi,sgs,ssor]
bar(iter,'b')
title(' STEADY ITERATION')
end

if state==2
%Transient state
for m=1:6
iteration=m;
%Transient state jacobi explicit
if iteration==1
tic
[t_ex_j,ejacobi]=explicitjacobi(nx,ny,t,told,tolr,ejacobi,error,k1,k2,dt,dx,dy,nt)
t_ex_jinsec=toc
figure(21)
contourf(xx,yy,t_ex_j,'ShowText','on')
colormap('jet')
colorbar
xlabel('x-axis')
ylabel('y-axis')
title('explicitjacobi')
end
%Transient state gauss-siedel explicit
if iteration==2
tic
[t_ex_gs,egs]=explicitgs(nx,ny,t,told,tolr,egs,error,k1,k2,dt,dx,dy,nt)
t_ex_gsinsec=toc
figure(22)
contourf(xx,yy,t_ex_gs,'ShowText','on')
colormap('jet')
colorbar
xlabel('x-axis')
ylabel('y-axis')
title('explicitgs')
end
%Transient state sor explicit
if iteration==3
tic
[t_ex_sor,esor]=explicitsor(w,nx,ny,t,told,tolr,esor,error,k1,k2,dt,dx,dy,nt)
t_ex_sorinsec=toc
figure(23)
contourf(xx,yy,t_ex_sor,'ShowText','on')
colormap('jet')
colorbar
xlabel('x-axis')
ylabel('y-axis')
title('explicitsor')
end
%Transient state jacobi implicit
if iteration==4
tic
[t_im_j,ijacobi]=implicitjacobi(nx,ny,t,tprev,told,tolr,ijacobi,error,k1,k2,dt,dx,dy,nt)
t_im_jinsec=toc
figure(24)
contourf(xx,yy,t_im_j,'ShowText','on')
colormap('jet')
colorbar
xlabel('x-axis')
ylabel('y-axis')
title('implicitjacobi')
end
%Transient state gauss-siedel implicit
if iteration==5
tic
[t_im_gs,igs]=implicitgs(nx,ny,t,told,tprev,tolr,igs,error,k1,k2,dt,dx,dy,nt)
t_im_gsinsec=toc
figure(25)
contourf(xx,yy,t_im_gs,'ShowText','on')
colormap('jet')
colorbar
xlabel('x-axis')
ylabel('y-axis')
title('implicitgs')
end
%Transient state sor implicit
if iteration==6
tic
[t_im_sor,isor]=implicitsor(nx,ny,t,told,tprev,tolr,isor,error,k1,k2,w,dt,dx,dy,nt)
t_im_sorinsec=toc
figure(26)
contourf(xx,yy,t_im_sor,'ShowText','on')
colormap('jet')
colorbar
xlabel('x-axis')
ylabel('y-axis')
title('implicitsor')
end

end
 % comparing time taken for jacobi , gauss-siedel and sor in transient state 
figure(27)
sec=[t_ex_jinsec,t_ex_gsinsec,t_ex_sorinsec,t_im_jinsec,t_im_gsinsec,t_im_sorinsec]
bar(sec,'r')
title('TIME')
 % comparing iteration taken for jacobi , gauss-siedel and sor in transient state 
figure(28)
iter=[ejacobi,egs,esor,ijacobi,igs,isor]
bar(iter,'b')
title('ITERATION')
end

end

A grid of 10x10 elements is defined with temperature values at top=600, bottom=900,left=400 and right=800.The values at the intersection of the rows are given by taking the average of the two values.Values for the constants alpha ,w,nt,tolr are initialised. A for loop is used to solve for 2 states( steady and transient state) If statement is used to select the state and each state is solved by means of using iterative solvers jacobi,sor and gauss siedel.The solvers are called as functions and the temperature(t) and iteration value are taken as output. The Tic and toc commands enclose the function to record the time it took for solving. The contourf command is used to generate a plot with x-axis representing xx values and y-axis representing yy values with temperature values plotted on them The colormap command is used to represent the temperature values in appropriate colors.The colorbar command is used to show the relation between the colors and values.

Steady state Jacobi

function [t,sjacobi]=steadyjacobi(nx,ny,t,told,tolr,sjacobi,error,k1,k2,dt,dx,dy,nt)
k=(2*((dx)^2+(dy)^2)/((dx^2)*(dy^2)));
    while(error>tolr)
      for i=2:nx-1
    for j=2:ny-1
        
        term1=(told(i-1,j)+told(i+1,j))/(dx^2);
        term2=(told(i,j-1)+told(i,j+1))/(dy^2);
        t(i,j)=((term1+term2)/k)
       end
      end
     error=max(max(abs(told-t)));
    sjacobi=sjacobi+1
    told=t;
    end
 end

A spacial loop is created by nx and ny values to find the values of temperature at every element in the grid using the equation.The error values are computed and a while loop is used to check for convergence by comparing derived error value for each set of iteration with tolrerance value.The told values are updated after every spacial loop or else the error value will never end and result in a infinite loop

Steady state Gauss-Siedel

function [t,sgs]=steadygs(nx,ny,t,told,tolr,sgs,error,k1,k2,dt,dx,dy,nt)
k=(2*((dx)^2+(dy)^2)/((dx^2)*(dy^2)));
    while(error>tolr)
      for i=2:nx-1
    for j=2:ny-1
        
        term1=(t(i-1,j)+t(i+1,j))/(dx^2);
        term2=(t(i,j-1)+t(i,j+1))/(dy^2);
        t(i,j)=((term1+term2)/k)
       end
      end
    error=max(max(abs(told-t)));
    sgs=sgs+1
    told=t;
    end
 end

A small change is made told values are replaced by t values so as and when a new value is found the updated value can be used for solving and reduce the iteration required for convergence.

Steady state SOR

function [t,ssor]=steadysor(w,nx,ny,t,told,tolr,ssor,error,k1,k2,dt,dx,dy,nt)
tgs=t;
k=(2*((dx)^2+(dy)^2)/((dx^2)*(dy^2)));
    while(error>tolr)
      for i=2:nx-1
    for j=2:ny-1
        
        term1=(t(i-1,j)+t(i+1,j))/(dx^2);
        term2=(t(i,j-1)+t(i,j+1))/(dy^2);
        tgs(i,j)=((term1+term2)/k)
        t(i,j)=((1-w)*told(i,j))+(w*tgs(i,j));
       end
      end
    error=max(max(abs(told-t)));
    ssor=ssor+1
    told=t;
    end
 end

The Gauss siedel method is used to find the value of temperature at next step and using the slope between the new value and previous value the equation is extrapolated to aid convergence.

Transient state explicit Jacobi

function [t,ejacobi]=explicitjacobi(nx,ny,t,told,tolr,ejacobi,error,k1,k2,dt,dx,dy,nt)
while(error>tolr)
         for i=2:nx-1
    for j=2:ny-1
        term1=told(i,j);
        term2=k1*(told(i-1,j)-2*told(i,j)+told(i+1,j));
        term3=k2*(told(i,j-1)-2*told(i,j)+told(i,j+1));
        t(i,j)=term1+term2+term3
       end
         end
         error=max(max(abs(told-t)));
      ejacobi=ejacobi+1
      told=t;
   end   

end

Transient state explicit Gauss-Siedel

function [t,egs]=explicitgs(nx,ny,t,told,tolr,egs,error,k1,k2,dt,dx,dy,nt)
while(error>tolr)
      for i=2:nx-1
    for j=2:ny-1
        term1=told(i,j);
        term2=k1*(t(i-1,j)-2*t(i,j)+t(i+1,j));
        term3=k2*(t(i,j-1)-2*t(i,j)+t(i,j+1));
        t(i,j)=term1+term2+term3;
       end
      end
      egs=egs+1
      error=max(max(abs(told-t)));
      told=t;
    
end   

end

Transient state explicit SOR

function [t,esor]=explicitsor(w,nx,ny,t,told,tolr,esor,error,k1,k2,dt,dx,dy,nt)
tgs=t;
while(error>tolr)
          for i=2:nx-1
    for j=2:ny-1
        term1=told(i,j);
        term2=k1*(t(i-1,j)-2*t(i,j)+t(i+1,j));
        term3=k2*(t(i,j-1)-2*t(i,j)+t(i,j+1));
        tgs(i,j)=term1+term2+term3;
        t(i,j)=((1-w)*told(i,j))+(w*tgs(i,j));
       end
          end
          esor=esor+1
      error=max(max(abs(told-t)));
       told=t;
   end   

end

Transient state implicit Jacobi

In implicit method the Tprevious values are updated at the end ot each time step and does not get changed with every iteration.

function [t,ijacobi]=implicitjacobi(nx,ny,t,tprev,told,tolr,ijacobi,error,k1,k2,dt,dx,dy,nt)
for z=1:nt
    while(error>tolr)
      for i=2:nx-1
    for j=2:ny-1
        term1=((1+(2*k1)+(2*k2))^-1);
        term2=k1*term1;
        term3=k2*term1;
        H=(told(i-1,j)+told(i+1,j));
        V=(told(i,j-1)+told(i,j+1));
        t(i,j)=(tprev(i,j)*term1)+(H*term2)+(V*term3);
       end
      end
    error=max(max(abs(told-t)));
    ijacobi=ijacobi+1
    told=t;
    end
    tprev=t;
    error=9e9;
end   

end

Transient state implicit Gauss Siedel

function [t,igs]=implicitgs(nx,ny,t,told,tprev,tolr,igs,error,k1,k2,dt,dx,dy,nt)
for z=1:nt
    while(error>tolr)
      for i=2:nx-1
    for j=2:ny-1
        term1=((1+(2*k1)+(2*k2))^-1);
        term2=k1*term1;
        term3=k2*term1;
        H=(t(i-1,j)+t(i+1,j));
        V=(t(i,j-1)+t(i,j+1));
        t(i,j)=(tprev(i,j)*term1)+(H*term2)+(V*term3);
       end
      end
    error=max(max(abs(told-t)));
    igs=igs+1
    told=t;
    end
    tprev=t;
    error=9e9;
end   

end

Transient state implicit SOR

function [t,isor]=implicitsor(nx,ny,t,told,tprev,tolr,isor,error,k1,k2,w,dt,dx,dy,nt)
tgs=t;
for z=1:nt
    while(error>tolr)
      for i=2:nx-1
    for j=2:ny-1
        term1=((1+(2*k1)+(2*k2))^-1);
        term2=k1*term1;
        term3=k2*term1;
        H=(t(i-1,j)+t(i+1,j));
        V=(t(i,j-1)+t(i,j+1));
        tgs(i,j)=(tprev(i,j)*term1)+(H*term2)+(V*term3);
        t(i,j)=((1-w)*told(i,j))+(w*tgs(i,j));
       end
      end
    error=max(max(abs(told-t)));
    isor=isor+1
    told=t;
    end
    tprev=t;
    error=9e9;
end   

end

PLOTS:

Conclusion:

The number of iteration required for solving steady state is lesser than transient state.

The steady stote sor method solves the equation at least number of iteration.

In transient state The implicit method require higher iteration to reach convergence compared to explicit method.

The number of iteration required decrease in the following order of method ImplicitJacobi>ImplicitGauss siedel>Implicitsor>ExplicitJacobi>ExplicitGauss siedel>Explicitsor