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 

ω">ω

= 1.25

2D steady state heat conduction equation and its solver 

implicit iterative solver 

∂2T∂x2+∂2T∂y2=0">∂2T∂x2+∂2T∂y2=0

∂2T∂x2=T(i-1,j)+2T(i,j)+T(i+1,j)δx2">∂2T∂x2=T(i−1,j)+2T(i,j)+T(i+1,j)δx2

∂2T∂y2=T(i,j-1)+2T(i,j)+T(i,j+1)δy2">∂2T∂y2=T(i,j−1)+2T(i,j)+T(i,j+1)δy2

1 jacobi 

T(i,j)=(T(i-1,j)+T(i+1,j)k⋅dx2+T(i,j-1)+T(i,j-1)k⋅dy2)old">T(i,j)=(T(i−1,j)+T(i+1,j)k⋅dx2+T(i,j−1)+T(i,j−1)k⋅dy2)old

Where k=2(∂x2+∂y2)∂x2⋅∂y2">k=2(∂x2+∂y2)∂x2⋅∂y2

2. Gauss Seidel

T(i,j)=(T(i-1,j)old+T(i+1,j)k⋅dx2+T(i,j-1)old+T(i,j-1)k⋅dy2)">T(i,j)=(T(i−1,j)old+T(i+1,j)k⋅dx2+T(i,j−1)old+T(i,j−1)k⋅dy2)

k=2(∂x2+∂y2)∂x2⋅∂y2">k=2(∂x2+∂y2)∂x2⋅∂y2

3. Sucessive Over Relexation

T(i,j)=(1-ω)T(i,j)old+ω(T(i-1,j)old+T(i+1,j)k⋅dx2+T(i,j-1)old+T(i,j-1)k⋅dy2)">T(i,j)=(1−ω)T(i,j)old+ω(T(i−1,j)old+T(i+1,j)k⋅dx2+T(i,j−1)old+T(i,j−1)k⋅dy2)

k=2(∂x2+∂y2)∂x2⋅∂y2">k=2(∂x2+∂y2)∂x2⋅∂y2

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 

ω">ω

= 1.3

α">α

= 4

2D heat conduction transient - state equation and its solver

∂T∂t=α(∂2T∂x2+∂2T∂y2)">∂T∂t=α(∂2T∂x2+∂2T∂y2)

∂2T∂x2=T(i-1,j)+2T(i,j)+T(i+1,j)δx2">∂2T∂x2=T(i−1,j)+2T(i,j)+T(i+1,j)δx2

∂2T∂y2=T(i,j-1)+2T(i,j)+T(i,j+1)δy2">∂2T∂y2=T(i,j−1)+2T(i,j)+T(i,j+1)δy2

1. Explicit 

T(i,j)n=(1-2k1-2k2)⋅T(i,j)n-1+k1⋅(T(i-1,j)+T(i.j-1))n-1+k2⋅(T(i,j-1)+T(i.j+1))n-1">T(i,j)n=(1−2k1−2k2)⋅T(i,j)n−1+k1⋅(T(i−1,j)+T(i.j−1))n−1+k2⋅(T(i,j−1)+T(i.j+1))n−1

Where k1=α⋅d2tdx2">k1=α⋅d2tdx2

            k2=α⋅d2tdy2">k2=α⋅d2tdy2

2. Implicit - jacobi

term1=11+2k1+2k2">term1=11+2k1+2k2

term2=k1⋅term1">term2=k1⋅term1

term3=k2⋅term1">term3=k2⋅term1

T(i,j)n=term1⋅T(i,j)n-1+term2⋅(T(i-1,j)+T(i+1,j))old+term3⋅(T(i,j-1)+T(i,j+1))old">T(i,j)n=term1⋅T(i,j)n−1+term2⋅(T(i−1,j)+T(i+1,j))old+term3⋅(T(i,j−1)+T(i,j+1))old

3. Implicit - gauss seidel

term1=11+2k1+2k2">term1=11+2k1+2k2

term2=k1⋅term1">term2=k1⋅term1

term3=k2⋅term1">term3=k2⋅term1

T(i,j)n=term1⋅T(i,j)n-1+term2⋅(T(i-1,j)+T(i+1,j)old)+term3⋅(T(i,j-1)+T(i,j+1)old)">T(i,j)n=term1⋅T(i,j)n−1+term2⋅(T(i−1,j)+T(i+1,j)old)+term3⋅(T(i,j−1)+T(i,j+1)old)

3. Implicit - SOR 

term1=11+2k1+2k2">term1=11+2k1+2k2

term2=k1⋅term1">term2=k1⋅term1

term3=k2⋅term1">term3=k2⋅term1

T(i,j)n=(1-ω)T(i,j)old+ω(term1⋅T(i,j)n-1+term2⋅(T(i-1,j)+T(i+1,j)old)+term3⋅(T(i,j-1)+T(i,j+1)old))">T(i,j)n=(1−ω)T(i,j)old+ω(term1⋅T(i,j)n−1+term2⋅(T(i−1,j)+T(i+1,j)old)+term3⋅(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

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

Unsteady state problems can be solved by both implicit and Explicit method 

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