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

AIM:

To write a code to solve the 2D heat conduction equation by using the point iterative techniques. 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

To implement for solving implicit equations for both(i.e for steady and transient analysis) and explicit schemes (for Transient Analysis only).

To implement the following methods for solving implicit equations

1. Jacobi

2. Gauss-seidel

3. Successive over-relaxation

Absolute error criteria is 1e-4

DESCRIPTION:

With the given temperature conditions for boundaries we find the temperatures at corners as follows and use them to denote the temperatures at corners in the code

STEADY STATE 2D HEAT CONDUCTION EQUATION

The steady state 2D heat conduction equation is `(del^2T)/(delx^2)+(del^2T)/(dely^2)=0rArrnabla^2T=0`

Discretizing the above equation by central differencing scheme we get

`(Tl-2Tp+Tr)/(Deltax^2)+(Tt-2Tp+Tb)/(Deltay^2)=0`

Here `Tl=Tunderset(i-1,j)`

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

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

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

`Tp=Tunderset(i,j)`

simpifying above equation

`(Tl+Tr)/(Deltax^2)+(Tt+Tb)/(Deltay^2)=(2Tp)/(Deltax^2)+(2Tp)/(Deltay^2)`

`(Tl+Tr)/(Deltax^2)+(Tt+Tb)/(Deltay^2)=(2Tp)(1/(Deltax^2)+1/(Deltay^2))`

`(Tl+Tr)/(Deltax^2)+(Tt+Tb)/(Deltay^2)=(2Tp)((Deltax^2+Deltay^2)/(Deltax^2Deltay^2))`

`(Tl+Tr)/(Deltax^2)+(Tt+Tb)/(Deltay^2)=kTp`

Here `k=2((Deltax^2+Deltay^2)/(Deltax^2Deltay^2))`

So finally the equation looks as `Tp=(1/k)((Tl+Tr)/(Deltax^2)+(Tt+Tb)/(Deltay^2))`

For Jacobi iterative Method

`Tp^(n+1)=(1/k)((Tl+Tr)/(Deltax^2)+(Tt+Tb)/(Deltay^2))^n`

For Gauss Seidel iterative method

`Tp^(n+1)=(1/k)((Tl+Tr)/(Deltax^2)+(Tt+Tb)/(Deltay^2))^(n+1)`

For Successive over relaxation iterative method

`Tp^(n+1)=Tp^n(1-omega)+omega(1/k)((Tl+Tr)/(Deltax^2)+(Tt+Tb)/(Deltay^2))^(n+1)`

`omega>1`Over relaxation

`omega<1` Under relaxation

TRANSIENT 2D HEAT CONDUCTION EQUATION

The transient 2D heat conduction equation is `(delT)/(delt)=alpha((del^2T)/(delx^2)+(del^2T)/(dely^2))rArr(delT)/(delt)=alphanabla^2T`

EXPLICIT METHOD

Discretizing for time derivative with forward differencing and with central differencing for space derivative we get `(Tp^(n+1)-Tp^n)/(Deltat)=alpha((Tl-2Tp+Tr)/(Deltax^2)+(Tt-2Tp+Tb)/(Deltay^2))^n`

Simplifying above equation we get

`Tp^(n+1)=Tp^n+DeltaTalpha((Tl-2Tp+Tr)/(Deltax^2)+(Tt-2Tp+Tb)/(Deltay^2))^n`

`Tp^(n+1)=Tp^n+(DeltaTalpha)/(Deltax^2)(Tl-2Tp+Tr)^n+(DeltaTalpha)/(Deltay^2)(Tt-2Tp+Tb)^n`

Here let `k1=(alphaDeltaT)/(delx^2)` and `k2=(alphaDeltaT)/(dely^2)`

Now the equation becomes `Tp^(n+1)=Tp^n+k1(Tl-2Tp+Tr)^n+k2(Tt-2Tp+Tb)^n`

`Tp^(n+1)=Tp^n(1-2k1-2k2)+k1(Tl+Tr)^n+k2(Tt+Tb)^n`

Where `CFLx=k1`and `CFLy=k2`

`CFL=CFLx+CFLy`

For explicit diffusion problems `CFL<=0.5`(stable i.e solution doesn't blows)

`CFL>0.5`(unstable i.e sloution blows)

IMPLICIT METHOD

Discretizing the transient equation time derivative with forward differencing and the space derivative with central differencing the equation is `(Tp^(n+1)-Tp^n)/(Deltat)=alpha((Tl-2Tp+Tr)/(Deltax^2)+(Tt-2Tp+Tb)/(Deltay^2))^(n+1)`

Simplifying above equation we get

`Tp^(n+1)=Tp^n+DeltaTalpha((Tl-2Tp+Tr)/(Deltax^2)+(Tt-2Tp+Tb)/(Deltay^2))^(n+1)`

`Tp^(n+1)=Tp^n+(DeltaTalpha)/(Deltax^2)(Tl-2Tp+Tr)^(n+1)+(DeltaTalpha)/(Deltay^2)(Tt-2Tp+Tb)^(n+1)`

Here let `k1=(alphaDeltaT)/(delx^2)` and `k2=(alphaDeltaT)/(dely^2)`

Now the equation becomes `Tp^(n+1)=Tp^n+k1(Tl-2Tp+Tr)^(n+1)+k2(Tt-2Tp+Tb)^(n+1)`

`Tp^(n+1)=Tp^n+k1(Tl+Tr)^(n+1)-2k1Tp^(n+1)+k2(Tt+Tb)^(n+1)-2k2Tp^(n+1)`

`Tp^(n+1)+2k1Tp^(n+1)+-2k2Tp^(n+1)=Tp^n+k1(Tl+Tr)^(n+1)+k2(Tt+Tb)^(n+1)`

`Tp^(n+1)(1+2k1+2k2)=Tp^n+k1(Tl+Tr)^(n+1)+k2(Tt+Tb)^(n+1)`

`Tp^(n+1)=(1/(1+2k1+2k2))(Tp^n+k1(Tl+Tr)^(n+1)+k2(Tt+Tb)^(n+1))`

For Jacobi iteration method

`Tp^(n+1)=(1/(1+2k1+2k2))(Tp^n+k1(Tl^n+Tr^n)+k2(Tt^n+Tb^n))`

For Gauss Seidel iteration method

`Tp^(n+1)=(1/(1+2k1+2k2))(Tp^n+k1(Tl^(n+1)+Tr^n)+k2(Tt^n+Tb^(n+1))`

For successive over relaxation method

`Tp^(n+1)=Tp^n(1-omega)+omega((1/(1+2k1+2k2))(Tp^n+k1(Tl^(n+1)+Tr^n)+k2(Tt^n+Tb^(n+1)))`

MAIN CODE FOR THE STEADY AND UNSTEADY 2D CONDUCTION EQUATION

%code for 2nd order PDE steady state heat conduction
%d^2u/dx^2+d^2u/dy^2=0

clear all
close all
clc

%inputs defining the no of nodes for grid size 10m along x and y direction
nx=10;
ny=nx;

%defining temperature array 1xnx and boundary conditions
T=300*ones(nx,ny); % assuming the inside temperature of grid
T(1,:) =600; %top boundary temperature
T(end,:) = 900; %bottom boundary temperature
T(:,1) = 400; %left boundary temperature
T(:,end) = 800; %right boundary temperature

%calculating temperature at edges taking average
T(1,1)= (600+400)/2; %top left node
T(1,end)= (600+800)/2; %top right node
T(end,1)= (900+400)/2; %bottom left node
T(end,end)= (900+800)/2; %bottom right node

%defining the grid size and mesh size along x and y axis with mesh size
x=linspace(0,1,nx); %along x no of nodes
dx=x(2)-x(1);%mesh size
y=linspace(0,1,ny);
dy=y(2)-y(1);
dt=1e-3;

%defining coefficient of conduction 
k = 2*((dx^2+dy^2)/((dx^2)*(dy^2))); %for steady state
alpha = 1.4;
k1 = (alpha*dt)/(dx^2); %CFL along x for transient state
k2 = (alpha*dt)/(dy^2); %CFL along y for transient state
cfl = k1+k2;

for method=linspace(1,7,7)
    if method ==1
        error1 = steadystate_jacobi(x,y,nx,ny,dx,dy,T,k);
        fprintf('Error using jacobi method for steady state=%d\n',error1);
    elseif method ==2
        error2 = steadystate_Gaussseidel(x,y,nx,ny,dx,dy,T,k);
        fprintf('Error using Gauss seidel method for steady state=%d\n',error2);
    elseif method ==3
        error3 = steadystate_sor(x,y,nx,ny,dx,dy,T,k);
        fprintf('Error using SOR method for steady state=%d\n',error3);
     elseif method ==4
         trans_explicit = transientstate_jacobiexplicit(x,y,nx,ny,T,k1,k2);
         fprintf('Error using jacobi method for explicit transient state=%d\n',trans_explicit);
     elseif method ==5
         error5 = transientstate_jacobiimplicit(x,y,nx,ny,T,k1,k2);
         fprintf('Error using jacobi method for implicit transient state=%d\n',error5);
     elseif method ==6
         error6 = transientstate_Gaussseidelimplicit(x,y,nx,ny,T,k1,k2);
         fprintf('Error using Gauss seidel method for implicit transient state=%d\n',error6);
     elseif method==7
         error7 = transientstate_sorimplicit(x,y,nx,ny,T,k1,k2);
         fprintf('Error using sor method for implicit transient state=%d\n',error7);
    end
end

FUNCTION CODE FOR STEADY IMPLICIT JACOBI METHOD

function [error]=steadystate_jacobi(x,y,nx,ny,dx,dy,T,k)

error = 9e9; %taken as default for comparision
tolerance = 1e-4; %given tolerance to run for
Told = T;
jacob_iteration =1;%defining as 1st iteration

while(error>tolerance)
    for i=2:nx-1
        for j=2:ny-1
            term1 = (Told(i-1,j)+Told(i+1,j))/(k*dx^2);
            term2 = (Told(i,j-1)+Told(i,j+1))/(k*dy^2);
            T(i,j)=term1+term2;
        end
    end
    error=max(max(abs(Told-T)));
    Told=T; %updating Told for next iteration
    jacob_iteration=jacob_iteration+1; %to run for next iteration
end

%plotting the figure
figure(1)
meshgrid(x,y);
[a,b]=contourf(x,y,T);
clabel(a,b)
colorbar %to get colors according to temperature strip
colormap("jet") %to get jet colouring
set(gca,'ydir','reverse')
xlabel('Grid length in x direction')
ylabel('Grid length in y direction')
title_text=sprintf('No of iterations for steady state jacobi method =%d',jacob_iteration);
title(title_text)

FUNCTION CODE FOR STEADY IMPLICIT GAUSS SEIDEL METHOD

function [error]=steadystate_Gaussseidel(x,y,nx,ny,dx,dy,T,k)

error = 9e9; %taken as default for comparision
tolerance = 1e-4; %given tolerance to run for
Told = T;
gauss_iteration =1;%defining as 1st iteration

while(error>tolerance)
    for i=2:nx-1
        for j=2:ny-1
            term1 = (T(i-1,j)+Told(i+1,j))/(k*dx^2);
            term2 = (T(i,j-1)+Told(i,j+1))/(k*dy^2);
            T(i,j)=term1+term2;
        end
    end
    error=max(max(abs(Told-T)));
    Told=T; %updating Told for next iteration
    gauss_iteration=gauss_iteration+1; %to run for next iteration
end

%plotting the figure
figure(2)
meshgrid(x,y);
[a,b]=contourf(x,y,T);
clabel(a,b)
colorbar %to get colors according to temperature strip
colormap("jet") %to get jet colouring
set(gca,'ydir','reverse')
xlabel('Grid length in x direction')
ylabel('Grid length in y direction')
title_text=sprintf('No of iterations for steady state Gauss seidel method =%d',gauss_iteration);
title(title_text)

FUNCTION CODE FOR STEADY IMPLICIT SUCCESSIVE OVER RELAXATION METHOD

function [error]=steadystate_sor(x,y,nx,ny,dx,dy,T,k)

error = 9e9; %taken as default for comparision
tolerance = 1e-4; %given tolerance to run for
omega=1.2; %over relaxation>1 under relaxation<1  
Told = T;
sor_iteration =1;%defining as 1st iteration

while(error>tolerance)
    for i=2:nx-1
        for j=2:ny-1
            term1 = (T(i-1,j)+Told(i+1,j))/(k*dx^2);
            term2 = (T(i,j-1)+Told(i,j+1))/(k*dy^2);
            T(i,j)=((1-omega)*Told(i,j))+(omega*(term1+term2));
        end
    end
    error=max(max(abs(Told-T)));
    Told=T; %updating Told for next iteration
    sor_iteration=sor_iteration+1; %to run for next iteration
end

%plotting the figure
figure(3)
meshgrid(x,y);
[a,b]=contourf(x,y,T);
clabel(a,b)
colorbar %to get colors according to temperature strip
colormap("jet") %to get jet colouring
set(gca,'ydir','reverse')
xlabel('Grid length in x direction')
ylabel('Grid length in y direction')
title_text=sprintf('No of iterations for steady state SOR method =%d',sor_iteration);
title(title_text)

FUNCTION CODE FOR TRANSIENT EXPLICIT

function [error]= transientstate_jacobiexplicit(x,y,nx,ny,T,k1,k2)

Told =T;
for nt=1:1400
    for i=2:nx-1
        for j=2:ny-1
            term1 = (Told(i-1,j)+Told(i+1,j));
            term2 = (Told(i,j-1)+Told(i,j+1));
           T(i,j)=(Told(i,j)*(1-(2*k1)-(2*k2)))+(k1*term1)+(k2*term2);
        end
    end
   
    error = max(max(abs(Told-T)))
    Told=T;
end

%plotting the figure
figure(4)
meshgrid(x,y);
[a,b]=contourf(x,y,T);
clabel(a,b)
colorbar %to get colors according to temperature strip
colormap("jet") %to get jet colouring
set(gca,'ydir','reverse')
xlabel('Grid length in x direction')
ylabel('Grid length in y direction')
title('Transient state Explicit Jacobi method');

FUNCTION CODE FOR TRANSIENT IMPLICIT JACOBI METHOD

function [error]= transientstate_jacobiimplicit(x,y,nx,ny,T,k1,k2)

Told=T;
Tpre=T;
jacobi_iteration =1;

for nt=1:1400
    error =9e9;
    tolerance=1e-4;
    while(error>tolerance)
        for i=2:nx-1
            for j=2:ny-1
                term1=1/(1+(2*k1)+(2*k2));
                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)=(Tpre(i,j)*term1)+(H*term2)+(V*term3);
            end
        end
       error =max(max(abs(Told-T)));
       Told=T;
       jacobi_iteration=jacobi_iteration+1;
    end
    Tpre=T;
end
figure(5)
meshgrid(x,y);
[a,b]=contourf(x,y,T);
clabel(a,b)
colorbar %to get colors according to temperature strip
colormap("jet") %to get jet colouring
set(gca,'ydir','reverse')
xlabel('Grid length in x direction')
ylabel('Grid length in y direction')
title_text=sprintf('No of iterations for Transient state Implicit Jacobi method =%d',jacobi_iteration);
title(title_text)

FUNCTION CODE FOR TRANSIENT IMPLICIT GAUSS SEIDEL METHOD

function [error]= transientstate_Gaussseidelimplicit(x,y,nx,ny,T,k1,k2)

Told=T;
Tpre=T;
gauss_iteration =1;

for nt=1:1400
    error =9e9;
    tolerance=1e-4;
    while(error>tolerance)
        for i=2:nx-1
            for j=2:ny-1
                term1=1/(1+(2*k1)+(2*k2));
                term2=k1*term1;
                term3=k2*term1;
                H=(T(i-1,j)+Told(i+1,j));
                V=(T(i,j-1)+Told(i,j+1));
                T(i,j)=(Tpre(i,j)*term1)+(H*term2)+(V*term3);
            end
        end
       error =max(max(abs(Told-T)));
       Told=T;
       gauss_iteration=gauss_iteration+1;
    end
    Tpre=T;
end
figure(6)
meshgrid(x,y);
[a,b]=contourf(x,y,T);
clabel(a,b)
colorbar %to get colors according to temperature strip
colormap("jet") %to get jet colouring
set(gca,'ydir','reverse')
xlabel('Grid length in x direction')
ylabel('Grid length in y direction')
title_text=sprintf('No of iterations for Transient state Implicit gauss siedel method =%d',gauss_iteration);
title(title_text)

FUNCTION CODE FOR TRANSIENT IMPLICIT SUCCESSIVE OVER RELAXATION METHOD

function [error]= transientstate_sorimplicit(x,y,nx,ny,T,k1,k2)

Told=T;
Tpre=T;
sor_iteration =1;
omega=1.2;

for nt=1:1400
    error =9e9;
    tolerance=1e-4;
    while(error>tolerance)
        for i=2:nx-1
            for j=2:ny-1
                term1=1/(1+(2*k1)+(2*k2));
                term2=k1*term1;
                term3=k2*term1;
                H=(T(i-1,j)+Told(i+1,j));
                V=(T(i,j-1)+Told(i,j+1));
                T(i,j)=(1-omega)*(Told(i,j))+(omega*((Tpre(i,j)*term1)+(H*term2)+(V*term3)));
            end
        end
       error =max(max(abs(Told-T)));
       Told=T;
       sor_iteration=sor_iteration+1;
    end
    Tpre=T;
end
figure(7)
meshgrid(x,y);
[a,b]=contourf(x,y,T);
clabel(a,b)
colorbar %to get colors according to temperature strip
colormap("jet") %to get jet colouring
set(gca,'ydir','reverse')
xlabel('Grid length in x direction')
ylabel('Grid length in y direction')
title_text=sprintf('No of iterations for Transient state Implicit SOR method =%d',sor_iteration);
title(title_text)

RESULTS:

For steady jacobi method,the number of iterations for covergence of solution it took is 208iterations

For steady Gauss seidel method,the number of iterations it took to converge the solution with tolerance is 112iterations

For steady successive over relaxation,the number of iterations it took is 76 to converge the solution

For Transient Explicit

if CFL<0.5, here it was 0.22 so the solution is coverged for dt=1e-3

if we change the dt to 0.1 the CFL changes to 22 which is >0.5 so the soltion will be unstable i.e blows up as below

For Transient Jacobi method,the number of iterations for jacobi method is 3666

For Transient Gausss siedel method,the number of iterations for gauss siedel is 3117

For Transient successive over relaxation method,the number of iterations for successive over relaxation is 2991

CONCLUSION:

Thus from the project we can conslude that the number of iterations run is in order Jacobi method>Gauss seidel method>successive over relaxation and if CFL<0.5 the results will be converged and stable if >0.5 then it blows up.