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

In this challenge, we analyze the 2D heat conduction equation for both steady state and transient state.

At first we take up the 'steady state' analysis of the 2D heat conduction equation

1. 2D heat conduction equation  =  d^2T/dx^2  + d^2T/dy^2 = 0

This is an exclusive example of 'elliptic PDE' , here 2D implies that this equation has 2 independent variable, steady state implies that it doesn't change with respect to time, hence we won't we 'time varible' in this equation only space varibles.

since this PDE is elliptic its physical behaviour changes correspondingly, it needs boundary conditions at all boundaries

To solve this equation commanly we go for matrix operations that is 'inverse sparse system' which is an analytical approach but as we know this type of approach is costly and inefficient, hence we use a 'numerical tecnique' called as ' point iterative' methods (jacobi, gauss seidel, SOR)

we use central space scheme in numerical discretization  

d^2T/dx^2  +  d^2T/dy^2 = 0

to find the number of nodes N = O + A -1 , 2 + 2 -1 , N = 3 

• (Ti+1,j -2Ti,j + Ti-1,j )/ dx^2  + (Ti,j+1 -2Ti,j + Ti,j-1)/dy^2 = 0

upon simplifying we get   Ti,j = 1/4[Ti+1,j + Ti-1,j + Ti,j+1 + Ti,j-1] this changes for jacobi, gauss seidel,SOR

 for jacobi =  Ti,j = 1/4[T^old i+1,j + T^old i-1,j + T^old i,j+1 + T^old i,j-1]

 for  gauss seidel = Ti,j = 1/4[T^new i+1,j + T ^newi-1,j + T^new i,j+1 + T^new i,j-1]

for SOR(succession over relaxation) = Ti,j = (1-r)*T^old i,j + r*T_GS

The program of this 2D heat diffussion equation is shared below using insert code option and the explanation is given here

First we give inputs such as number of grid points nx, ny and T(temparature) to form matrix

then  apply boundary conditions TL,TR,TT,TB. and populate it suitable values , also Tinitial,Told ,tolrence and error values

then we go for iterative solvers accordingly one by one to do that first open while loop i:e run our program till (error > tolrenace value)

then since its a two dimension open  2for loops for i an j and write jacobi method solver equation close 2for loops and then caluculate a absolute error for Told and Tinitial

to go for next iteration it needs updated T values hence we give Told = T and the close a while loop then run a program and calculate similary with gauss seidel and SOR 

when we run program we observe how absolute error values decreases and approaches 1 , also number iterations decreases from jacobi to SOR  jacobi take more number of iterations than gauss seidel and SOR, because jacobi uses old values as a reference for succesive iterations but gauss seidel uses the newly updated values for further iterations  

% 2D heat conduction on a flat plate
clear all
close all
clc

nx = 30;
ny = 30;
T = 300*ones(nx,ny)

% Boundary conditions
TL = 400;
TR = 800;
TT = 600;
TB = 900;

T(:,1) = TL;
T(:,end) = TR;
T(1,:) = TT;
T(end,:) = TB;

T(1,1) = 0.5*(TL+TT);
T(1,nx) = 0.5*(TR+TT);
T(ny,1) = 0.5*(TL+TB);
T(nx,ny) = 0.5*(TR+TB);

T_intial = T;
T_old = T;

tolerance = 1e-4;
error = 1;

% jacobi method
jacobi_iter = 0;
while(error>tolerance)
     for i = 2:nx-1
    for j = 2:ny-1
TL = T_old(i,j-1); TR = T_old(i,j+1); TT = T_old(i-1,j); TB = T_old(i+1,j);
T(i,j) = 0.25*(TL +TR +TT +TB);
    end
    end
    error = max(max(abs(T-T_old)))
    T_old = T;
    jacobi_iter = jacobi_iter+1;
end

% reinitialise

T_GS = T_intial;
T_old = T_intial;

tolerance = 1e-4;
error = 1;

% gauss siedel method
 iter_GS = 0;
 while(error>tolerance)
     for i = 2:nx-1
    for j = 2:ny-1
TL = T_GS(i,j-1); TR = T_GS(i,j+1); TT = T_GS(i-1,j); TB = T_GS(i+1,j);
T_GS(i,j) = 0.25*(TL +TR +TT +TB);
    end
    end
    error = max(max(abs(T_GS-T_old)))
    T_old = T_GS;
    iter_GS = iter_GS+1;
end

% Re-initialize fields

T_SOR = T_intial;
T_old = T_intial;

iter_SOR = 0;
tolerance = 1e-4;
error = 1;

r = 0.5

%SOR method

while(error>tolerance)
    for j =2:ny-1
    for i =2:nx-1
   TL = T_GS(i,j-1); TR = T_GS(i,j+1); TT = T_GS(i-1,j); TB = T_GS(i+1,j);
   T_GS(i,j) = 0.25*(TL +TR +TT +TB);
   T_SOR(i,j) = (1-r)*T_old(i,j) + r*T_GS(i,j);
    end
    end
    error = max(max(abs(T_SOR - T_old)));
    T_old = T_SOR;
    iter_SOR = iter_SOR + 1;
end

Now we shall see the analysis of the 2D heat conduction of unsteady state or transient state by both explicit and implicit methods

here 2D implies that it has two independent variables and transient state implies that this equation changes with respect to time also

2D heat conduction equation of an unsteady state = pCp dT/dt = k [ d^2T/dx^2  +  d^2T/dy^2 ]

 p = density ,Cp = heat capacity, k = thermal conductivity, alpha = k/pCp  = 'thermal diffusivity'

Heat capacity = amount of heat supplied to an object to produce a unit change in its temparature.

Thermal conductivity = a measure of its ability to conduct heat in an object 

Thermal diffusivity = measures the rate of transfer of heat of a material from hot end to cold end

This is an exclusive example of parabolic PDE[partial differential equation] it requires both initial and boundary conditions 

For Explicit method numerical discretization scheme used is [FTCS] i:e forward time central space scheme

 Final equation would be =  T^n+1 (i,j) = (1-2*k1 -2*k2)T^n(i,j) + k1*(T^n(i+1,j) + T^n(i-1,j)) + k2*(T^n(i,j+1) + T^n(i,j-1))

here k1 = aplha*dt/dx^2 ,k2 = alpha*dt/dy^2 which is CFL number therefore we can say if CFL <0.5 then the equation is stable if CFL>0.5 it is unstable

in explicit time marching we take two matrices for converging we consider previous space and time values and approximate.

For implicit method numerical discretization scheme used is [BTCS] backward time cenral space scheme

Final equation would be = T^n+1 (i,j) = 1/(1+2k1+2k2) * T^n (i,j) + k1*(T^n+1(i+1,j) + T^n+1(i-1,j)) + k2*(T^n+1(i,j+1) + T^n+1(i,j-1))

In implicit time marching we take three marices for approxomation or to converge, here first we converge space values by taking previous matrix space values, then we converge time values by taking present matrix into consideration  

The program for both implicit and explicit methods are solved and given below seperately

plots for temprature have been shared as image and attached it below

%% 2D HEAT CONDUCTION ON A FLAT PLATE TRANSIENT STATE
clear all
close all
clc

%% EXPLICIT METHOD


nx = 10;
ny = 10;
Lx = 1;
Ly = 1;

dx = Lx/nx;
dy = Ly/ny;

% x = linspace(0,Lx,nx);
% y = linspace(0,Lx,ny);
% T = meshgrid(x,y);
T = 300*ones(nx,ny);

TL = 400;
TR = 800;
TT = 600;
TB = 900;

T(:,1) = TL;
T(:,end) = TR;
T(1,:) = TT;
T(end,:) = TB;

T(1,1) = 0.5*(TL + TT);
T(1,nx) = 0.5*(TR + TT);
T(ny,1) = 0.5*(TL + TB);
T(nx,ny) = 0.5*(TR + TB);

error = 1;
tol = 1e-4;

% CFL

dt = 1e-2;
CFL = 0.1;
alpha = (CFL*dx^2)/dt;

k1 = alpha*dt/dx^2;
k2 = alpha*dt/dy^2;

%copy
% Previous time matrix
T_prev = T;
% Current time matrix 
T_initial = T;

% EXPLICIT SOLVER

iter = 0;
tspan = 0;

% CONVERGENCE CRITERION FOR TIME
while(error>tol) % FOR LOOP (REPEATS UNTIL nt)

    for j = 2:ny-1
    for i = 2:nx-1
    % Previous time level neighbours
    TL = T_prev(i,j-1); TR = T_prev(i,j+1); TT = T_prev(i-1,j); TB = T_prev(i+1,j);
    H = TL + TR; V = TT + TB;
    k1k2 = 1-2*k1-2*k2;
    T(i,j) = T_prev(i,j)*k1k2 + H*k1 + V*k2;
    end
    end

  iter = iter + 1;
  tspan = tspan + (dt*iter);
  error(iter) = max(max(abs(T-T_prev)));
  T_prev = T;

end

%% 2D HEAT CONDUCTION ON A FLATE PLATE TRANSIENT PLATE
clear all
close all
clc

%% IMPLICIT METHOD

nx = 10;
ny = 10;
Lx = 1;
Ly = 1;

dx = Lx/nx;
dy = Ly/ny;

% x = linspace(0,Lx,nx);
% y = linspace(0,Lx,ny);
% T = meshgrid(x,y);
% Current time level matrix
T = 300*ones(nx,ny);

TL = 400;
TR = 800;
TT = 600;
TB = 900;

T(:,1) = TL;
T(:,end) = TR;
T(1,:) = TT;
T(end,:) = TB;

T(1,1) = 0.5*(TL + TT);
T(1,nx) = 0.5*(TR + TT);
T(ny,1) = 0.5*(TL + TB);
T(nx,ny) = 0.5*(TR + TB);

tol = 1e-4;
Terror =1;

% CFL

dt = 1;
alpha = 1e-5;

k1 = alpha*dt/dx^2;
k2 = alpha*dt/dy^2;

CFL = k1 + k2;

% copy 
% Previous time matrix
T_prev = T;
% Previous iteration current time matrix
T_old = T;
% Initial matrix
T_initial = T;

% IMPLICIT SOLVER

Titer = 0;
Siter = 0;
tspan = 0;
Serror = 1;

% Convergence criterion for time
% TIME LOOP
while(Terror>tol)
    Titer = Titer + 1;
    Serror = 1;
% SPACE LOOP - JACOBIAN
while(Serror>tol)

    for j = 2:ny-1
    for i = 2:nx-1
% Current time level but old iteration neighbours
TL = T_old(i,j-1); TR = T_old(i,j+1); TT = T_old(i-1,j); TB = T_old(i+1,j);
    H = TL + TR; V = TT + TB;
    k1k2 = 1+2*k1+2*k2;
    T(i,j) = (T_prev(i,j) + H*k1 + V*k2)/k1k2;
    end
    end

Siter = Siter + 1;
Serror(Siter) = max(max(abs(T-T_old)));
T_old = T;

% disp('space iter \n')
% disp(Siter)

end

tspan = tspan + (dt*Titer);
Terror(Titer) = max(max(abs(T-T_prev)));
T_prev = T;

% disp('Time iter \n')
%disp(Titer)

end