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

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