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

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

GIVEN PROBLEM:

1. Steady state analysis & Transient State Analysis

Solve the 2D heat conduction equation by using the point iterative techniques that were taught in the class. The Boundary conditions for the problem are as follows;

Top Boundary = 600 K

Bottom Boundary = 900 K

Left Boundary = 400 K

Right Boundary = 800 K

You will implement the following methods for solving implicit equations.

1. Jacobi

2. Gauss-seidel

3. Successive over-relaxation

Your absolute error criteria is 1e-4

Make sure to use "INSERT CODE SAMPLE" option to paste the code here. Post screenshots of the plots and use the projects.skill-lync.com page to Publish your findings.

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

3. Write a detailed explanation for your results, plots or codes without explanation will awarded zero marks.

Theory:

conduction is a mode of heat transfer generally takes place in solids due to temperature difference associated with

free electron and lattice vibration. The equation is given by:

Thermal diffusivity is a property of a material defined as thermal conductivity and thermal capacity

since it a 2D heat conduction problem and there is no internal heat generation,so the equation is parabolic.

For steady state time derivative is not constant,so the equation is eliptical

To obtain numeric solution of 2D steady and unsteady heat conduction equation difference equations are

steady state 2-D heat conduction:

Transient 2-D heat conduction equation(explicit):

Transient 2-D heat conduction (implicit):

CODE FOR STEADY-STATE EQUATION:

FUNCTION CODE:

JACOBI:

function out=jacobi(T,T_old,i,j)
nx=10;
ny=10;
x=linspace(0,1,nx);
dx=x(2)-x(1);
y=linspace(0,1,ny);
dy=y(2)-y(1);
k=(2*dx^2+2*dy^2)/(dx^2*dy^2);
term1=(T_old(i+1,j)+T_old(i-1,j))/dx^2;
term2=(T_old(i,j+1)+T_old(i,j-1))/dy^2;
out=(term1+term2)/k;

GAUSS:

function out=GS(T,T_old,i,j)
nx=10;
ny=10;
x=linspace(0,1,nx);
dx=x(2)-x(1);
y=linspace(0,1,ny);
dy=y(2)-y(1);
k=(2*dx^2+2*dy^2)/(dx^2*dy^2);
term1=(T_old(i+1,j)+T(i-1,j))/dx^2;
term2=(T_old(i,j+1)+T(i,j-1))/dy^2;
out=(term1+term2)/k;

SOR:

unction out=sor(T,T_old,i,j,omega)
nx=10;
ny=10;
x=linspace(0,1,nx);
dx=x(2)-x(1);
y=linspace(0,1,ny);
dy=y(2)-y(1);
k=(2*dx^2+2*dy^2)/(dx^2*dy^2);
term1=(T_old(i+1,j)+T(i-1,j))/dx^2;
term2=(T_old(i,j+1)+T(i,j-1))/dy^2;
out=((1-omega)*T_old(i,j))+(omega*((term1+term2)/k));

MAIN CODE:

clear 
close all
clc

%inputs
%grid points along x and y axis
nx=10;
ny=10;
x=linspace(0,1,nx);
dx=x(2)-x(1);
y=linspace(0,1,ny);
dy=y(2)-y(1);


%initial temperature values for all nodes
T=300*ones(nx,ny);

%initial boundary condition at the top and botttom
T(1,:)=600;
T(end,:)=900;
%initial boundary condition at the left and right
T(:,1)=400;
T(:,end)=800;

%copy temperature to compute
T_old=T;
T_prev=T;
% initializing error value
error=9e9;
tolerance=1e-4;
omega=1.2;
for iterative_solver=1:3
   
    if  iterative_solver==1
         iter=1;
    while error>tolerance
        for i=2:nx-1
         for   j=2:ny-1             
         T(i,j)=jacobi(T,T_old,i,j);
         end
        end
    
         
          error=max(max(abs(T_old-T)));
    %updating value of T
         T_old=T;
         figure(1)
        
      [u,v]=meshgrid(x,y);
     [c,h]= contourf(u,v,T);
       colorbar
       colormap("jet")
       clabel(c,h)
        set(gca,'ydir','reverse')
 
       iter_method=sprintf('jacobi iteration is:%d',iter);
       title(iter_method);
        iter=iter+1;
    end
   



    elseif iterative_solver==2
        T_old=T_prev;
     iter1=1;
    while error>tolerance
        for i=2:nx-1
         for   j=2:ny-1             
         T(i,j)=GS(T,T_old,i,j);
         end
        end
    
         
          error=max(max(abs(T_old-T)));
    %updating value of T
         T_old=T;
         figure(2)
        
      [u,v]=meshgrid(x,y);
      [c,h]=contourf(u,v,T);
       colorbar
       colormap("jet")
       set(gca,'ydir','reverse')
       clabel(c,h)
 
       iter_method=sprintf('GS iteration is:%d',iter1);
       title(iter_method);
       iter1=iter1+1;
    end

  



    elseif iterative_solver==3
        iter2=1;
        T_old=T_prev;
    while error>tolerance
        for i=2:nx-1
         for   j=2:ny-1             
         T(i,j)=sor(T,T_old,i,j,omega);
         end
        end
    
         
          error=max(max(abs(T_old-T)));
    %updating value of T
         T_old=T;
         figure(3)
        
      [u,v]=meshgrid(x,y);
     [c,h]= contourf(u,v,T);
       colorbar
       colormap("jet")
        set(gca,'ydir','reverse')
        clabel(c,h);
 
       iter_method=sprintf('sor iteration is:%d',iter2);
       title(iter_method);
       iter2=iter2+1;
    end
    end
    
    
error=9e9;
T_old=T;
end

output:

for jacobi:

for gs:

for sor:

EXPLANATION:

1.First of all we take n number of nodes along x and y axis as nx and ny . Distance between two nodes along x and y axis are

dx and dy. we take temperature at each point is 300k and then take boundary conditions like top and bottom temperature are

600k and 900k and left and right are 400k and 800k.

2.We take two copy of initial conditions of temperature, T as T_old and T_prev.Take a random value of error is 9e^9 and tolerance is 1e^-4 so as error should be greater than tolrance as we have to reduce error value up to tolerance so as to achieve perfection from imperfection.Also taking the value of omega as 1.2 , if we incraese omega more so we mayface more iteration thus incorporating wrong result.

3.we have to adopt a iterative_solver (a variable)to check this method for three formulas as jacobi,gauss-seidal and sor.

Vary the iterative_solver from 1 to 3 and check each formula by if-else condition.

4.In previous its discussed that we should decrease error upto level of tolerance so we adopt a while loop where error>tolerance.Then as its a 2D problem so we have to take two for loop as i & j and vary it from 2 to (n-1) as boundary

conditions are known.

5.After for loops we are updating error value as maximum difference of T_old and newly derived T. we have to update

T_old=T as the newly derived T will be T_old in the loop for the next iteration. In each loop we have to take iteration value

before while loop and increment it by "itreration=iteration+1" before while loop finishes.These all are done at each iterative_solver.

6.We are taking x,y command as 2d form by meshgrid command and plot it by contourf command.We are taking a column like entity to show how the colors are distributed by colorbar command , we can make different types of shades of color by

colormap command.we are pointing sizes of each plot by clabel command.

7.After declaring title_method by sprintf we are taking it as title of the plot.

8.At iterative_solver=2 before while loop we are taking T_old=T_prev, as the T_old should be first value before iterative_solver=1 ,so that it can start iteration from the first T_old value otherwise T_old value will be taken after it had been

iterated at the end of the iterative_solver=1. same is done for the iterative_solver=3.

CODE FOR NON-STEADY STATE EQUATION:

clear all
close all
clc
%inputs
%number of grid points along x and y axis
nx=10;
ny=10;
x=linspace(0,1,nx);
y=linspace(0,1,ny);
dx=x(2)-x(1);
dy=y(2)-y(1);
nt=1400;
alpha=1.4;
omega=1.2;
dt=0.001;
k1=alpha*(dt/dx^2);
k2=alpha*(dt/dy^2);
t=300*ones(nx,ny);

%applying boundary conditions
t(1,:)=600;
t(end,:)=900;
t(:,1)=400;
t(:,end)=800;
told=t;
t_prev=t;
iterative_method=input('Enter 1 for implicit or enter 2 for explicit');
if iterative_method==1
for iterative_solver=1:3

    if iterative_solver==1
        iter=1;      
 for nt=1:1400
    error=9e9;
    tol=1e-4;
    while error>tol
        for i=2:nx-1
            for j=2:ny-1
                term1=1/(1+2*k1+2*k2);
                term2=k1*(told(i+1,j)+told(i-1,j));
                term3=k2*(told(i,j+1)+told(i,j-1));
                t(i,j)=term1*(t_prev(i,j)+term2+term3);
            end
        end
        error=max(max(abs(told-t)));
        told=t;
        iter=iter+1;
    end
   t_prev=t;
 end
 figure(1)
        
        [c,h]=contourf(x,y,t);
         clabel(c,h);
        colorbar
        colormap("jet")
        told=t;
      
        title_text=sprintf("jacobi iteration:%d",iter);
        title(title_text);
        set(gca,'ydir','reverse')


    elseif iterative_solver==2
        iter1=1;      
 for nt=1:1400
    error=9e9;
    tol=1e-4;
    while error>tol
        for i=2:nx-1
            for j=2:ny-1
                term1=1/(1+2*k1+2*k2);
                term2=k1*(told(i+1,j)+t(i-1,j));
                term3=k2*(told(i,j+1)+t(i,j-1));
                t(i,j)=term1*(t_prev(i,j)+term2+term3);
            end
        end
        error=max(max(abs(told-t)));
        told=t;
        iter1=iter1+1;
    end
   t_prev=t;
 end
 figure(2)
        
        [c,h]=contourf(x,y,t);
         clabel(c,h);
        colorbar
        colormap("jet")
        told=t;
      
        title_text=sprintf("GS iteration:%d",iter);
        title(title_text);
       set(gca,'ydir','reverse')


    elseif iterative_solver==3
        iter2=1;      
 for nt=1:1400
    error=9e9;
    tol=1e-4;
    while error>tol
        for i=2:nx-1
            for j=2:ny-1
                 term1=1/(1+2*k1+2*k2);
                term2=k1*(told(i+1,j)+t(i-1,j));
                term3=k2*(told(i,j+1)+t(i,j-1));
                term4=term1*(t_prev(i,j)+term2+term3);
                t(i,j)=(1-omega)*told(i,j)+(omega*term4);
            end
        end
        error=max(max(abs(told-t)));
        told=t;
        iter2=iter2+1;
    end
   t_prev=t;
 end
 figure(3)
        
        [c,h]=contourf(x,y,t);
        clabel(c,h);
        colorbar
        colormap("jet")
        told=t;
      set(gca,'ydir','reverse')
        title_text=sprintf("sor  iteration:%d",iter);
        title(title_text);
    end
end

elseif iterative_method==2
    
iter3=1;
for nt=1:1400
    error=9e9;
    tol=1e-4;
    while error>tol
        for i=2:nx-1
            for j=2:ny-1
                term1=k1*(t_prev(i+1,j)-2*t_prev(i,j)+t_prev(i-1,j));
                term2=k2*(t_prev(i,j+1)-2*t_prev(i,j)+t_prev(i,j-1));
                t(i,j)=t_prev(i,j)+term1+term2;
            end
        end
         error=max(max(abs(told-t)));
        told=t;
        iter3=iter3+1;
    end
    t_prev=t;
end
 figure(4)
        
        [c,h]=contourf(x,y,t);
         clabel(c,h);
        colorbar
        colormap("jet")
        told=t;
      set(gca,'ydir','reverse')
        title_text=sprintf("explicit  iteration:%d",iter3);
        title(title_text);
end

OUTPUT:

IMPLICIT:

JACOBI:

GAUSS:

SOR:

EXPLICIT:

EXPLANATION:

All are same like steady state explanation but we are taking an extra time loop nt=1400 for unsteady state.We are declaring unsteady-state by implicit and explicit by iterative_method( a variable) where iterative_method=1 for implicit and 2 for explicit.

CONCLUSION:

1.We can conclude that for steady or unsteady sor is faster than gauss and then gauss is faster than jacobi.

2.Then for unsteady-state explicit method is much faster than implicit.