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

AIM: Solving the steady and unsteady 2D heat conduction problem

OBJECTIVE: 

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 are 1e-4.

Part 1 - 

Steady-state heat conduction

Input conditions to solve the 2D heat conduction at a steady-state using an implicit technique 

Length = 1m

Grid points along with x-axis = 100

Grid points along with y-axis = 100

Temperature at left = 400K 

Temperature at Right = 500K 

Temperature at Top = 600K 

Temperature at Bottom = 900K 

$\omega$ = 1.25

2D steady state heat conduction equation and its solver 

implicit iterative solver 

$\frac{{\partial }^{2}T}{\partial x^2}+\frac{{\partial }^{2}T}{\partial y^2}=0$

$\frac{{\partial }^{2}T}{\partial x^2}=\frac{T(i-1,j)+2T(i,j)+T(i+1,j)}{\delta x^2}$

$\frac{{\partial }^{2}T}{\partial y^2}=\frac{T(i,j-1)+2T(i,j)+T(i,j+1)}{\delta y^2}$

1 jacobi 

$T(i,j)={(\frac{T(i-1,j)+T(i+1,j)}{k\cdot {}^{}}+\frac{T(i,j-1)+T(i,j-1)}{k\cdot {}^{}})}^{old}$

Where $k=\frac{2(\partial x^2+\partial y^2)}{\partial x^2\cdot \partial y^2}$

2. Gauss Seidel

$T(i,j)=(\frac{T{(i-1,j)}^{old}+T(i+1,j)}{k\cdot {dx}^2}+\frac{T{(i,j-1)}^{old}+T(i,j-1)}{k\cdot {dy}^2})$

$k=\frac{2(\partial x^2+\partial y^2)}{\partial x^2\cdot \partial y^2}$

3. Sucessive Over Relexation

$T(i,j)=(1-\omega )T{(i,j)}^{old}+\omega (\frac{T{(i-1,j)}^{old}+T(i+1,j)}{k\cdot {dx}^2}+\frac{T{(i,j-1)}^{old}+T(i,j-1)}{k\cdot {dy}^2})$

$k=\frac{2(\partial x^2+\partial y^2)}{\partial x^2\cdot \partial y^2}$

MATLAB CODE 

% 2D heat conduction at steady state using implicit technique 

clear all
close all
clc
L=1;%Length of the domain
omega=1.25;% SOR factor
nx=100;%grid size
ny=nx;
x=linspace(0,L,nx);
dx=x(2)-x(1);
y=linspace(0,L,ny);
dy=y(2)-y(1);
T=300*ones(nx,ny);
T(:,1)=400;
T(1,:)=900;
T(:,end)=800;
T(end,:)=600;
T(1,1)=(900+400)/2;
T(1,end)=(900+800)/2;
T(end,1)=(400+600)/2;
T(end,end)=(800+600)/2;
Told=T;
k1=2*(dx^2+dy^2)/(dx^2*dy^2);
tol=1e-4;
error=9e9;
solver=input('Enter the number for the solver method: n Jacobi=1 n Gauss Seidel=2 n SOR=3 n Solver Type: ');
if solver == 1
   tic
   jacobi_iter=1;
   while (error>tol)
    for j=2:ny-1
      for i=2:nx-1
        H=((Told(i-1,j)+Told(i+1,j))/(k1*dx^2));
        V=((Told(i,j-1)+Told(i,j+1))/(k1*dy^2));
        T(i,j)=H+V;
      end
    end
    error=max(max(abs(Told-T)))
    Told=T;
    jacobi_iter=jacobi_iter+1
    timecounter=toc;
    figure(1)
    [C,h]=contourf(x,y,T,'ShowText','on');
    colorbar
    colormap(jet)
    xlabel('xx');
    ylabel('yy');
    clabel(C,h)
    title(sprintf('2d Steady State Heat Conduction  n Method:Jacobi n No of iterations=%d n Time:%f s',jacobi_iter,timecounter))
   end
end
   if solver == 2
   tic
   gs_iter=1;
   while (error>tol)
    for j=2:ny-1
      for i=2:nx-1
        H=((T(i-1,j)+Told(i+1,j))/(k1*dx^2));
        V=((T(i,j-1)+Told(i,j+1))/(k1*dy^2));
        T(i,j)=H+V;
      end
    end
    error=max(max(abs(Told-T)))
    Told=T;
    gs_iter=gs_iter+1
    timecounter_1=toc;
    figure(1)
    [C,h]=contourf(x,y,T,'ShowText','on');
    colorbar
    colormap(jet)
    xlabel('xx');
    ylabel('yy');
    clabel(C,h)
    title(sprintf('2d Steady State Heat Conduction  n Method:Gauss Seidel n No of iterations=%d n Time:%f s',gs_iter,timecounter_1))
   end
   end
      if solver == 3
   tic
   sor_iter=1;
   while (error>tol)
    for j=2:ny-1
      for i=2:nx-1
        H=((T(i-1,j)+Told(i+1,j))/(k1*dx^2));
        V=((T(i,j-1)+Told(i,j+1))/(k1*dy^2));
        T(i,j)=(Told(i,j)*(1-omega))+(omega*(H+V));
      end
    end
    error=max(max(abs(Told-T)))
    Told=T;
    sor_iter=sor_iter+1
    timecounter_2=toc;
    figure(1)
    [C,h]=contourf(x,y,T,'ShowText','on');
    colorbar
    colormap(jet)
    xlabel('xx');
    ylabel('yy');
    clabel(C,h)
    title(sprintf('2d Steady State Heat Conduction  Method:SOR  No of iterations=%d  Time:%f s',sor_iter,timecounter_2))
   end
   end

RESULTS

Part 2 - Transient state heat conduction 

Input condition to solve 2D heat conduction under transient state using implicit & explicit technique 

Length = 1m

Grid points along with x-axis = 51

Grid points along with y-axis = 51

Temperature at left = 400K 

Temperature at Right = 500K 

Temperature at Top = 600K 

Temperature at Bottom = 900K 

$\omega$ = 1.3

$\alpha$ = 4

2D heat conduction transient - state equation and its solver

$\frac{\partial T}{\partial t}=\alpha (\frac{{\partial }^{2}T}{\partial x^2}+\frac{{\partial }^{2}T}{\partial y^2})$

$\frac{{\partial }^{2}T}{\partial x^2}=\frac{T(i-1,j)+2T(i,j)+T(i+1,j)}{\delta x^2}$

$\frac{{\partial }^{2}T}{\partial y^2}=\frac{T(i,j-1)+2T(i,j)+T(i,j+1)}{\delta y^2}$

1. Explicit 

$T{(i,j)}^n=(1-2k\underset{1}{}-2k\underset{2}{})\cdot T{(i,j)}^{n-1}+k\underset{1}{}\cdot {(T(i-1,j)+T(i.j-1))}^{n-1}+k\underset{2}{}\cdot {(T(i,j-1)+T(i.j+1))}^{n-1}$

Where $k\underset{1}{}=\frac{\alpha \cdot d^{2}t}{{dx}^2}$

            $k\underset{2}{}=\frac{\alpha \cdot d^{2}t}{{dy}^2}$

2. Implicit - jacobi

$term1=\frac{1}{1+2k\underset{1}{}+2k\underset{2}{}}$

$term2=k1\cdot term1$

$term3=k2\cdot term1$

$T{(i,j)}^n=term1\cdot T{(i,j)}^{n-1}+term2\cdot {(T(i-1,j)+T(i+1,j))}^{old}+term3\cdot {(T(i,j-1)+T(i,j+1))}^{old}$

3. Implicit - gauss seidel

$term1=\frac{1}{1+2k\underset{1}{}+2k\underset{2}{}}$

$term2=k1\cdot term1$

$term3=k2\cdot term1$

$T{(i,j)}^n=term1\cdot T{(i,j)}^{n-1}+term2\cdot (T(i-1,j)+T{(i+1,j)}^{old})+term3\cdot (T(i,j-1)+T{(i,j+1)}^{old})$

3. Implicit - SOR 

$term1=\frac{1}{1+2k\underset{1}{}+2k\underset{2}{}}$

$term2=k1\cdot term1$

$term3=k2\cdot term1$

$T{(i,j)}^n=(1-\omega )T{(i,j)}^{old}+\omega (term1\cdot T{(i,j)}^{n-1}+term2\cdot (T(i-1,j)+T{(i+1,j)}^{old})+term3\cdot (T(i,j-1)+T{(i,j+1)}^{old}))$

%% Transient state 2D heat conduction using Explicit technique 
clear all
close all
clc

%% intializing variables
Nx = 51;
Ny = 51;
x = linspace (0,1,Nx);
y = linspace (0,1,Ny);
dx = 1/(Nx-1);
dy = 1/(Ny-1);
alpha = 0.0001;
dt = 1;

tol = 1e-4;

%% Temperature gride 
T = 298*ones(Nx,Ny);
T(1,:) = 900; %bottom temp
T(end,:) = 600; %Top temp
T(:,1) = 400; %Left temp
T(:,end) = 800; %right temp
%edge refining
T(1,1) = (900+400)/2;
T(1,end) = (900+800)/2;
T(end,1) = (400+600)/2;
T(end,end) = (800+600)/2;

%% CFL number calculation and adaptive Time-step control
CFL_number = (alpha*dt)*((1/dx^2)+(1/dy^2));

%% Calculation of temperature distribution explicitly 

k1 = (alpha*dt)/(dx^2);
k2 = (alpha*dt)/(dy^2);
[xx,yy] = meshgrid(x,y);
Told = T;
T_prev_dt = T;
tic;
counter = 1;
for z = 1:2500
    % Explicit Iteration Scheme
    for i = 2:Nx-1
        for j = 2:Ny-1
            term1 = (1-2*k1-2*k2)*Told(i,j);
            term2 = k1*(Told(i-1,j)+Told(i+1,j));
            term3 = k2*(Told(i,j-1)+Told(i,j+1));
            T(i,j) = term1+term2+term3;
        end
    end
    Told = T;
    counter = counter+1;
    time_counter = toc;
    figure (1)
    [C,h] = contourf(xx,yy,T);
    colorbar;
    colormap(jet);
    clabel(C,h);
    title(sprintf('2D Unstedy state heat conduction  Method:Explicit Method  Number of iterations: %d(Time: %0.2fs)  Computation Time: %0.4f  CFL Number: %0.4f', counter-1,z*dt,time_counter,CFL_number));
    xlabel('xx');
    ylabel('yy');
end 

Result

Transient state heat conduction using implicit method

%% 2D transient state heat conduction equation solved using implicit method
clear all
close all
clc

%% Initializing Variable
Nx = 51;
Ny = 51;
x = linspace (0,1,Nx);
y = linspace (0,1,Ny);
dx = 1/(Nx-1);
dy = 1/(Ny-1);
alpha = 0.0001;
dt = 1;

tol = 1e-4;

%% Temperature gride 
T = 298*ones(Nx,Ny);
T(1,:) = 900; %bottom temp
T(end,:) = 600; %Top temp
T(:,1) = 400; %Left temp
T(:,end) = 800; %right temp
%edge refining
T(1,1) = (900+400)/2;
T(1,end) = (900+800)/2;
T(end,1) = (400+600)/2;
T(end,end) = (800+600)/2;

%% CFL number calculation and adaptive Time-step control
CFL_number = (alpha*dt)*((1/dx^2)+(1/dy^2));

%% Calculation of temperature distribution using implicit method

k1 = (alpha*dt)/(dx^2);
k2 = (alpha*dt)/(dy^2);
[xx,yy] = meshgrid(x,y);
Told = T;
T_prev_dt = T;

term1 = 1/(1+2*k1+2*k2);
term2 = k1*term1;
term3 = k2*term1;
tic;
counter = 1;
solver = input('Enter the number for the solver method: n Jacobi=1 n Gauss Seidal=2 n SOR=3 n Solver type: ');
% Jacobin Method
if solver == 1
    for k = 1:2500
        error = 9e9;
        while(tol<error)
            for j = 2:Ny-1
                for i = 2:Nx-1
                    H = (Told(i-1,j)+Told(i+1,j));
                    V = (Told(i,j-1)+Told(i,j+1));
                    T(i,j) = (T_prev_dt(i,j)*term1)+(H*term2)+(V*term3);
                end 
            end
            error = max(max(abs(Told-T)));
            Told = T;
            counter = counter+1;
        end
        T_prev_dt = T;
        figure (1)
         [C,h] = contourf(xx,yy,T);
         colorbar;
         colormap(jet);
         clabel(C,h);   
         time_counter = toc;
         title(sprintf('2D Unsteady state heat conduction  Method:Implicit Jacobi  Number of iterations: %d (Time: %0.2fs)  Computation Time: %0.4f  CFL Number: %0.4f',counter-1,k*dt,time_counter,CFL_number));
         xlabel('xx');
         ylabel('yy');
    end 
end 

% Gauss Seidel Method
if solver == 2 
    for k = 1:2500
        error = 9e9;
        while (tol<error)
            for j = 2:Ny-1
                for i = 2:Nx-1
                    H = (T(i-1,j)+Told(i+1,j));
                    V = (T(i,j-1)+Told(i,j+1));
                    T(i,j) = (T_prev_dt(i,j)*term1)+(H*term2)+(V*term3);
                end 
            end
            error = max(max(abs(Told-T)));
            Told = T;
            counter = counter+1;
        end 
        T_prev_dt = T;
        figure (1)
        [C,h] = contourf(xx,yy,T);
     colorbar;
     colormap(jet);
     clabel(C,h);   
     time_counter = toc;
     title(sprintf('2D Unsteady state heat conduction  Method:Implicit Gauss Seidel  Number of iterations: %d (Time: %0.2fs)  Computation Time: %0.4f  CFL Number: %0.4f',counter-1,k*dt,time_counter,CFL_number));
     xlabel('xx');
     ylabel('yy');
    end
end

%SOR method
if solver == 3
    for k=1:2500
    error = 9e9;
    while(tol<error)
     for j = 2:Ny-1
        for i = 2:Nx-1
            H= (T(i-1,j)+Told(i+1,j));
            V= (T(i,j-1)+Told(i,j+1));
            T(i,j) =(1-1.3)*(Told(i,j))+(1.3*((T_prev_dt(i,j)*term1)+(H*term2)+(V*term3)));
        end
      end
     error=max(max(abs(Told-T)));
     Told = T;
     counter = counter+1;
     end
     T_prev_dt = T;
     figure(1)
     [C,h] = contourf(xx,yy,T);
     colorbar;
     colormap(jet);
     clabel(C,h);   
     time_counter = toc;
     title(sprintf('2D Unsteady state heat conduction n Method:Implicit SOR n Number of iterations: %d (Time: %0.2fs) n Computation Time: %0.4f n CFL number: %0.4f',counter-1,k*dt,time_counter,CFL_number));
     xlabel('xx');
     ylabel('yy');
    end
end 

RESULT

Conclusion 

PART 1: Steady-state heat conduction

From the above result, the iteration for convergence of plot is very high for the Jacobi method when compare with gauss- Seidel and SOR. By increasing the grid size we can drastically reduce the iteration which as well gives a convergence plot with stability.

PART 2 - transient- state heat conduction

The above result shows the convergence of transient heat equation using an explicit and implicit technique

1. Unsteady state problems can be solved by both implicit and Explicit method 
2. For stiff equations required lesser timesteps and larger computational time to solve using the explicit method.

Among implicit solvers, SOR gives the result in less time compared to other techniques. It gives stable results irrespective of the CFL number