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

AIM-Steady state analysis & Transient State Analysis

OBJECTIVE-To solve the 2D heat conduction equation by using point iterative techniques.

LET AX=B

`[[a,b],[c,d]]*[[x1],[x2]]=[[b1],[b2]]`

A=D+L+U

where

D=diagonal matrix

L=lower triangular matrix

U=upper triangular matrix

JACOBI TECHNIQUE

The solution obtained iteratively

`X^(k+1)=D^-1*(b-(L+U)^k)`

here `x^k`is the kth iteration `x^(k+1)`is the k+1 iteration

`xunderset(i)()^(k+1)=1/aunderset(ij)()(bunderset(i)()-sumaunderset(ij)()xunderset(j)()^k)`

GAUSS SEIDEL

linear equation

`L*x=b-Ux`

`x=L^-1(b-Ux^(k))`

Using forward substitution

`xunderset(i)()^(k+1)=1/aunderset(ij)()(bunderset(i)()-sum_(j=0)^(i-1)aunderset(ij)()xunderset(j)()^(k+1)-sum_(j=i+1)^naunderset(ij)()xunderset(j)()^(k))`

SUCESSIVE OVER RELAXATION

linear solution

`(D+omegaL)x=omegab-(omegau+(omega-1)D)x`

where `omega`is the relaxation factor

The solution obtained iteratively

`(xunderset(i)()^(k+1))=(D+omegaL)^-1*(omegab-[omegau+(omega-1)D]x^k)`

STEADY STATE HEAT CONDUCTION EQUATION

It is the form of conduction that happens when the temperature difference driving the conduction are constant so that the temperature field in te conducting object doesnot change any further.Thus, all the partial derivatives of temperature concernng space may be either zero or have non zero values, but all the derivatives of temperature at any point concerning time are uniformly zero.

`(del^2T)/(delx^2)+(del^2T)/(dely^2)=0`

SOLVING TECHNIQUE

`(del^2T)/(delx^2)+(del^2T)/(dely^2)=0`

`(Tunderset(L)-2*Tunderset(p)+Tunderset(R)())/(Deltax^2)+(Tunderset(T)()-2*Tunderset(p)()+Tunderset(B)())/(Deltay^2)=0`

`2*(Tunderset(p)())/(Deltax)^2+2*(Tunderset(p)())/(Deltay)^2=(Tunderset(L)()+Tunderset(R)())/(Deltax^2)+(Tunderset(T)()+Tunderset(B)())/(Deltay^2)`  

`Tunderset(P)()*(2*(Deltax^2+Deltay^2)/((Deltax^2*Deltay^2)))`=`(Tunderset(L)()+Tunderset(R)())/(Deltax^2)+(Tunderset(T)()+Tunderset(B)())/(Deltay^2)`

`Tunderset(P)()=1/k((Tunderset(L)()+Tunderset(R)())/(Deltax^2)+(Tunderset(T)()+Tunderset(B)())/(Deltay^2))`

where,k=`(2*(Deltax^2+Deltay^2)/((Deltax^2*Deltay^2)))`

here 

`Tunderset(L)()=Tunderset(i-1,j)()`

`Tunderset(R)()=Tunderset(i+1,j)()`

`Tunderset(T)()=Tunderset(i,j+1)()`

`Tunderset(B)()=Tunderset(i,j-1)()`

`Tunderset(P)()=Tunderset(i,,j)()`

CODE

close all
clear 
clc



%Solving steady state heat conduction using iterative methods

%inputs
%no of grid points in x and y axis 
nx=10;
ny=nx;
w=1.5;%relaxation factor

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

dx=x(2)-x(1);
dy=y(2)-y(1);


%Boundary conditions
Tp=300*ones(nx,ny);
Tp(1,:)=600;
Tp(:,end)=800;
Tp(end,:)=900;
Tp(:,1)=400;
Tp(1,1)=500;
Tp(1,end)=700;
Tp(end,end)=850;
Tp(end,1)=650;
k=2*((dx^2+dy^2)/(dx^2*dy^2));


T_old=Tp;


%iterative_solver=1;
%error


tol=1e-4;
error=9e9;
%selection statement
iterative_solver = input('Enter your choice (1 for Jacobi, 2 for GS, 3 for SOR) : ');


if iterative_solver==1
    j_iter=1;
   
    %convergence loop
    while error>tol
        %space loop
         for i=2:nx-1
             for j=2:ny-1
                   xterm=(T_old(i+1,j)+T_old(i-1,j))/(dx^2);
                   yterm=(T_old(i,j+1)+T_old(i,j-1))/(dy^2);
                   Tp(i,j)=(xterm+yterm)/k;
          
             end
        end
      error=max(max(abs(T_old-Tp)));
      T_old=Tp;
      j_iter=j_iter+1;
       
      
    end
    %plotting in graphical form
   
     contourf(x,y,Tp,'showtext','on');
      xlabel('x');
      ylabel('y');
      colorbar;
      axis ij;
    
      a2=sprintf('niteration=%d',j_iter);  
      title(['Gauss Jacobi',a2]);
  
     
end
 

if iterative_solver==2
    gs_iter=1;
    tic
    %convergence loop
    while error>tol
        %space loop
         for i=2:nx-1
             for j=2:ny-1
                   xterm=(Tp(i+1,j)+Tp(i-1,j))/(dx^2);
                   yterm=(Tp(i,j+1)+Tp(i,j-1))/(dy^2);
                   Tp(i,j)=(xterm+yterm)/k;
          
             end
        end
      error=max(max(abs(T_old-Tp)));
      T_old=Tp;
      gs_iter=gs_iter+1;
       
      
    end
    %plotting in graphical form
   
     contourf(x,y,Tp,'showtext','on');
      colorbar;
      axis ij;
      xlabel('x');
      ylabel('y');
      a2=sprintf('niteration=%d',gs_iter);  
      title(['Gauss Seidal',a2]);
  
     
end

if iterative_solver==3
   
    sor_iter=1;
    tic
    %convergence loop loop
    while error>tol
        %space loop
         for i=2:nx-1
             for j=2:ny-1
                   xterm=(Tp(i+1,j)+Tp(i-1,j))/(dx^2);
                   yterm=(Tp(i,j+1)+Tp(i,j-1))/(dy^2);
                   Tp(i,j)=(1-w)*T_old(i,j)+w*(xterm+yterm)/k;
          
             end
        end
      error=max(max(abs(T_old-Tp)));
      T_old=Tp;
      sor_iter=sor_iter+1;
       
      
    end
    %plotting in graphical form
      contourf(x,y,Tp,'showtext','on');
      xlabel('x');
      ylabel('y');
      colorbar;
      axis ij;
      a2=sprintf('niteration=%d',sor_iter);  
      title(['SOR METHOD',a2]);
  
      
end

OUTPUT

CLICKING 1 JACOBI ITERATION IS SOLVED

CLICKING 2 GAUSS SIEDEL ITERATION IS SOLVED

CLICKING 3 SOR ITERATION IS SOLVED

TRANSIENT STATE HEAT CONDUCTION EQUATION

During period in which temperatures changes in time at any place within an object, the mode of thermal energy flow is termed transient heat conduction.

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

solving procedure

Explicit way

`(Tunderset(p)()^(n+1)-Tunderset(p)()^n)/(Deltat)=alpha*((Tunderset(L)()-2*Tunderset(p)()+Tunderset(R)())/(Deltax^2)+(Tunderset(B)()-2*Tunderset(p)()+Tunderset(T)())/(Deltay^2))^(n)`

`Tunderset(p)()^(n+1)=Tunderset(p)()^n+(Deltatalpha)((Tunderset(L)()-2*Tunderset(p)()+Tunderset(R)())/(Deltax^2))+(Deltatalpha)((Tunderset(B)()-2*Tunderset(p)()+Tunderset(T)())/(Deltay^2))^(n)`

CFL number,k1=`CFLunderset(x)()=(alphaDeltat)/(Deltax^2)`

k2=`CFLunderset(y)()=(alphaDeltat)/(Deltay^2)`

if CFL<0.5 stable

if CFL>0.5 unstable

CODE

close all
clear all
clc
%solving unsteady 2D heat conduction equation using explicit scheme

%input
nx=10;
ny=nx;

%defining x and y axis
x=linspace(0,1,nx);
y=linspace(0,1,ny);
dx=x(2)-x(1);
dy=y(2)-y(1);

%time step
nt=1400;
dt=1e-4;

%boundary conditions
T=300*ones(nx,ny);
T(1,:)=600;
T(end,:)=900;
T(:,1)=400;
T(:,end)=800;
T(1,1)=500;
T(1,end)=700;
T(end,end)=850;
T(end,1)=650;
Told=T;
Tpre=T;
iter=1;


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

tol=1e-4;

%iterative statement
if iter==1
     explicit_iter=1;
     
    for k=1:nt
         error=9e9;
        %converging loop
        while error>tol
        %space loop   
        for i=2:nx-1
            
            for j=2:ny-1
                
             xterm=k1*((Told(i-1,j)-2*Told(i,j)+Told(i+1,j)));
             yterm=k2*((Told(i,j-1)-2*Told(i,j)+Told(i,j+1)));
             T(i,j)=Told(i,j)+xterm+yterm;
            end
        end
         
        error=max(max(abs(Told-T)));
        Told=T;
        explicit_iter=explicit_iter+1;
       
        end
    end
end
%plotting
figure(1)
contourf(x,y,T,'showtext','on');
colorbar;
colormap(jet)
set(gca,'YDIR','reverse')
xlabel('x')
ylabel('y')
a2=sprintf('niteration=%d',explicit_iter);  
title(a2);

OUTPUT

IMPLICIT WAY

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

`(Tunderset(p)()^(n+1)-Tunderset(p)()^n)/(Deltat)=alpha*((Tunderset(L)()-2*Tunderset(p)()+Tunderset(R)())/(Deltax^2)+(Tunderset(B)()-2*Tunderset(p)()+Tunderset(T)())/(Deltay^2))^(n+1)`

`Tunderset(p)()^(n+1)=1/(1+2k1+k2)*[Tunderset(p)()^n+kunderset(1)()*(Tunderset(L)()+Tunderset(R)())^(n+1)+k2*(Tunderset(L)()+Tunderset(R)())^(n+1)]`

%solving unsteady heat conduction equation implicit method
close all
clear all
clc

%inputs
nx=10;
ny=nx;



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

dx=x(2)-x(1);
dy=y(2)-y(1);

%relaxation factor
w=1.6;

%time step
nt=1400;
dt=0.001;
tol=1e-4;

%boundary conditions
T=300*ones(nx,ny);
T(1,:)=600;
T(end,:)=900;
T(:,1)=400;
T(:,end)=800;
T(1,1)=500;
T(1,end)=700;
T(end,end)=850;
T(end,1)=650;
Told=T;
Tpre=T;

%tol=1e-4;
iteration=input('enter 1 for Gauss Jacobi 2 for Gauss Seidal 3 for succesive relaxation:');

%thermal diffusion
alpha=1.5;

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

%jacobi iteration

if iteration==1
              tic 
               %iteration loop
              j_iter=1;
    for k=1:nt
        tic 
        error=9e9;
     %converging loop
       while (error>tol)
           %space loop
          for i=2:nx-1
              
                  for j=2:ny-1
                      
                       const=((1+(2*k1)+(2*k2)))^-1;
                       term1=k1*((Told(i+1,j)+Told(i-1,j)));
                       term2=k2*((Told(i,j+1)+Told(i,j-1)));
                      
                       T(i,j)=const*(Told(i,j)+term1+term2);
                       
                       
                 end
          end
            error=max(max(abs(Told-T)));
            Told=T;
            j_iter=j_iter+1;
       
       end
                
                Tpre=T;
    end    
                time=toc;
    
    %plotting
    figure(1)
    contourf(x,y,T,'showtext','on');
    colorbar;
    colormap(jet)
    set(gca,'YDIR','reverse')
    xlabel('x')
    ylabel('y')
  
    a2=sprintf('\niteration=%d',j_iter);  
    title(['jacobi',a2]);
    
end


%Gauss Seidal
if iteration==2
        tic
        Gs_iter=1;
    %space loop
    for k=1:nt
        tic
        error=9e9;
        
     
        while error>tol
          for i=2:nx-1
              for j=2:ny-1
                       const=((1+(2*k1)+(2*k2)))^-1;
                       term3=k1*((T(i+1,j)+T(i-1,j)));
                       term4=k2*((T(i,j+1)+T(i,j-1)));
                      
                       T(i,j)=const*(T(i,j)+term3+term4);
                       
              end
          end
        error=max(max(abs(Told-T)));
        Told=T;
        Gs_iter=Gs_iter+1;
        end
         
        Tpre=T;
        
    end
        time=toc;
    %plottong
    figure(1)
    contourf(x,y,T,'showtext','on');
    colorbar;
    colormap(jet)
    set(gca,'YDIR','reverse')
    xlabel('x')
    ylabel('y')
    
    a2=sprintf('\niteration=%d',Gs_iter);  
    title(['Gauss seidal',a2]);
end
        

OUTPUT

Conclusion 

It was expected that

Jacobi iterations >Gauss Seidal Iterations >SOR it is observed the same for Steady State.

i.e.206>112>26

Jacobi iterations >Gauss Seidal Iterations >SOR it is observed the same for UnSteady State.

i.e. 1986>1903>1692

thus can be concluded that SOR has the least iterations. 

By applying the appropriate value of relaxation factor the least iterations can be obtained.