## Solving 2-D 2nd Order Linear Transient &Steady  Heat Conduction Equation Using Both Explicit And Implicit Numerical Schemes In MATLAB.

Solving 2-D 2nd Order Linear Transient &Steady  Heat Conduction Equation Using Both Explicit And Implicit Numerical Schemes In MATLAB.

• Explicit Scheme : We say that a scheme is explit when the information at time level "n+1" depends on previous time levels i.e., "n , n-1" e.g. the upwind difference schemes for the convection equation is an explicit schemes.
• Implicit Scheme : In this schemes the information at time level "n+1" is not only dependent on previous time steps "n , n-1" but also on current time step "n+1" e.g.      Crank-Nicolson method, Euler Implicit methos etc.

The 2-D Transient Heat Conduction Equation is given by:  The numerical difference equation in implicit form is given as:-   where,  The iterative methods implemented to solve in the implicit way are:

• Gauss-Jacobi method.
• Gauss-Seidel method.
• Successive Over Relaxation method For gauss-Seidel method.

The Matlab Codes are:

% Implicit solver method for 2-D Heat Equation.

clear all
close all
clc

% Grid points in x-direction
nx = 30;

% Grid Points in y-direction
ny = nx;
x = linspace(0,1,nx);
y = linspace(0,1,ny);
% Time Step
dt = 1e-1;

% length of the domain
L_x = 1; % length in x - direction
L_y = 1; % Lengthin y- direction

% Mesh size
dx = x(2)-x(1);
dy = y(2)-y(1);

% Grid
T = ones(nx,ny);

% % Initial Temperature over the grid = 100K
T = 300*ones(nx,ny);
% T(1,:) = linspace(100,800,nx);
% T(:,1) = 100;
%
% for k = 2:nx;
%
%     for l = 2:ny;
%
%         T(k,l) = T(k-1,l-1)+100;
%
%     end
%
% end
%
% T_initial_before_boundary_condition = T;

% Applying Boundary condition
% Top_Edge = 600K
% Left_Edge = 400K
% Right_Edge = 800k
% Bottom_Edge = 900K

T(1,:) = 900;
T(end,:) = 600;

T(:,1) = 400;
T(:,end) =800;

T_initial = T;
T_old = T;

iteration  = input('enter input 1.Jacobi 2.Gauss-Seidel 3. SOR Gauss-Seidel')
tic;

% Defining constants
% Thermal Diffusivity of the material (aluminium)

alpha = 1e-5;
K_1 = ((alpha*dt)/(dx^2));
K_2 = ((alpha*dt)/(dy^2));
Term_1 = ((1+(2*K_1)+(2*K_2))^-1);
p = K_1*Term_1;
q = K_2*Term_1;

tolerance = 1e-4;
error = 9e9;

jacob_iter = 1;
% Jacobi Iteration

if iteration ==1;
for a = 1:20000

error_1 = error;

while(error_1 > tolerance)

for  i = 2:nx-1;

for j = 2:ny-1;

T(i,j) =(T_old(i,j)*Term_1) + ((p*(T_initial(i-1,j)+T_initial(i+1,j)))+(q*(T_initial(i,j-1)+T_initial(i,j+1))));

end

end

error_1 = max(max(abs(T_initial-T)))
T_initial = T;
jacob_iter = jacob_iter + 1
end

T_old = T;

end

toc;
F = toc;

%Plotting the curve
figure(1)
[m n] = contourf(x,y,T);
clabel(m,n);
colormap(jet);
colorbar;
zlabel('x-axis');
ylabel('y-axis');
[X,Y] = meshgrid(x,y);
title(sprintf('Transient State Heat Conduction(Gauss-Jacobi):Total Iteration= %d,Total Time = %fs',jacob_iter,F));

end

% gauss-Seidel Iteration
if iteration ==2;
Gauss_Seidel = 2;
tic;

for b = 1:20000

error = 9e9;
error_2 = error;

while(error_2 > tolerance)

for k = 2:nx-1;

for l = 2:ny-1;

T(k,l) = (T_old(k,l)*Term_1) + (p*(T(k-1,l)+T_initial(k+1,l))+(q*(T(k,l-1)+T_initial(k,l+1))));

end

end

error_2 = max(max(abs(T_initial - T)))
T_initial = T;
Gauss_Seidel = Gauss_Seidel+1

end

T_old = T;

end
toc;
F  = toc;
%Plotting the curve
figure(1)
[m n] = contourf(x,y,T);
clabel(m,n);
colormap(jet);
colorbar;
zlabel('x-axis');
ylabel('y-axis');
title(sprintf('Transient State Heat Conduction(Gauss-seidel):Total Iteration= %d,Total Time = %fs',Gauss_Seidel,F));
end

% Successive Over Relaxation Gauss-Sidel Iteration
SOR_Gauss_Seidel = 3;
tic;
if iteration ==3;
omega = 0.9;

for b = 1:20000
error = 9e9;
error_3 = error;

while(error_3 > tolerance)

for k = 2:nx-1;

for l = 2:ny-1;

T(k,l) = (T_old(k,l)*(1-omega)) + omega*((T_old(k,l)*Term_1) + (p*(T(k-1,l)+T_initial(k+1,l))+(q*(T(k,l-1)+T_initial(k,l+1)))));

end

end

error_3 = max(max(abs(T_initial - T)))
T_initial = T;
SOR_Gauss_Seidel = SOR_Gauss_Seidel+1

end

T_old = T;

end
toc;
F  = toc;
%Plotting the curve
figure(2)
[C,h] = contourf(x,y,T);
clabel(C,h);
colormap(jet);
colorbar;
zlabel('x-axis');
ylabel('y-axis');
title(sprintf('Transient State Heat Conduction(SOR Gauss-sidel):Total Iteration= %d n,Total Time = %f s',SOR_Gauss_Seidel,F));
end



Results(Plots):   Observation:

1. Gauss-jacobi methods takes total iterations of "41720" to converge the results.
2. Whereas, Gauss-Seidel methods takes total iterations of "40570" which is lesser than gauss-Jacobi iteration method.
3. However, using SOR method for Gauss-Seidel method we introduce a relaxation parameter "omega = 0.2" for faster convergence rate "omega ∈ (0,2)". the result converges at 40003 iteration which faster than the above iteration methods.
4. The method is stable  due to implicit approach of the solution.

The Explicit form of the 2-D Transient Heat equation: % Explicit Solver for 2-D Heat Equation.

clear all
close all
clc

% Making grid of square matrix

dt = 1e-3;

nx = 10;

ny = nx;

L = 1;
dx =L/nx-1;
dy = L/ny-1;
x = linspace(0,1,nx);
y = linspace(0,1,ny);

dy = dx;

T = 300* ones(nx,ny);

% T = linspace(200,1000,n)
%
% for a = 2:10
%
%     T(a,:) = T(a-1)+100
%
% end

% Applying Boundary condition
% Top_Edge = 600K
% Left_Edge = 400K
% Right_Edge = 800k
% Bottom_Edge = 900K

T(1,:) = 900;
T(end,:) = 600;
T(:,1) = 400;
T(:,end) =800;

T_initial = T;

alpha = 1e-5;
K_1 = ((alpha*dt)/(dx^2));
K_2 = ((alpha*dt)/(dy^2));

% The stability criteria of the method.
CFL = alpha*dt*((1/dx^2) + (1/dy^2))

% Time loop

Exp_iteration = 1;

tic;
for b = 1:10000

% Starting Time Loop
Exp_iteration = Exp_iteration;

for i = 2:nx-1

for j = 2:ny-1   % NESTED LOOPING

T(i,j) =(T_initial(i,j)) + ((K_1*(T_initial(i-1,j)-(2*T_initial(i,j))+(+T_initial(i+1,j)))+(K_2*(T_initial(i,j-1)-(2*T_initial(i,j))+T_initial(i,j+1)))))

end

end

T_initial = T;
Exp_iteration = Exp_iteration +1

end
toc;
F= toc;
figure(2)
[m n] = contourf(x,y,T);
clabel(m,n);
colormap(jet);
colorbar;
xlabel('x-axis');
ylabel('y-axis');
title('Steady State Heat conduction Explicit Method');


Results(plot): EFFECT OF TIME STEP ON RESULTS FOR EXPLICIT SOLVER:  OBSERVATIONS:

• Explicit approach is usually not as stable as that of Implicit approach as a result in explicit appraoch we need to use smaller time step in order to maintain the stability whereas we can use higher time step and maintain stability in implicit approach.
• The above approach is a time marching of temperature over a grid(plate) via explicit method in which new value of temperature is assigned to a node of time"n+1" from a previous time "n , n-1..".
• At time step = 1e+1, we see that the truncation error obtained = 0.0181 and in the plots the temperature ranges are not uniformly distributed where as for time step = 1e-1, the truncation error is =0.0013 and temperature ranges are more closer and more uniformly distributed.Hence, smaller the time step generates smaller truncation error and more accurate results are obtained.
• Also if the time step is increased beyond the optimum limit given as CFL (Courent Friedrichs-Lewy Number) then the solution gets unstable.

Steady State Heat Equation Analysis By Implicit Method: MATLAB CODES:

% Steady state heat equation
clear all
close all
clc

% defining the number of grid points in X,Y direction.
nx = 100;
ny = nx;

% Defining the domain length and mesh size.
L = 1;
dx = L/(nx-1);
dy = L/(ny-1);
x = linspace(0,L,nx);
y = linspace(0,L,ny);
[X Y] = meshgrid(x,y);

% Defining the Temperature distribition over the grid.
T = zeros(nx,ny);
T(1,:) = linspace(100,800,nx);
for a = 2:nx;

for b = 1:ny;

T(a,b) = T(a-1,b);

end
end

T_Before_boundary_condition = T;
T_initial = T;

% Applying Boundary Condition
% Top_Edge = 600K
% Left_Edge = 400K
% Right_Edge = 800k
% Bottom_Edge = 900K

T(1,:)   = 600;
T(end,:) = 900;
T(:,1)   = 400;
T(:,end) = 800;
iteration = input('1. Jacobi_Iteration  2. Gauss-sidel_iteration  3.SOR_Gauss-sidel_iteration');
k = (2*(dx^2 + dy^2)/(dx^2)*(dy^2));

% # Jacobi Iteration for the Steady State Heat Equation.
% # Initializing Convergence Loop
% # Tolerance for steaady state convergence.
Tolerance = 1e-10;

tic;
jacob_iter = 1;
if iteration == 1;

error   = 9e9;
error_1 = error;

while (error_1 > Tolerance);

for i = 2:nx-1;

for j = 2:ny-1;

T(i,j) =((T_initial(i-1,j) + T_initial(i+1,j))/k*(dx^2)) + ((T_initial(i,j-1) + T_initial(i,j+1))/k*(dy^2));

end

end

error_1 = max(max(abs(T_initial - T)))
T_initial = T;
jacob_iter = jacob_iter + 1

end
toc;
F = toc;

%Plotting the curve
figure(1)
[m n] = contourf(x,y,T);
clabel(m,n);
colormap(jet);
colorbar;
zlabel('x-axis');
ylabel('y-axis');
title(sprintf('Steady State heat Conduction(Jacobi):total iteration= %d, Time taken = %fs',jacob_iter,F));
end

% initializing Gauss-sidel Iteration
tic;
Gauss_Sidel_iter = 2;
if iteration == 2;

error   = 9e9;
error_2 = error;
while (error_2 > Tolerance);

for i = 2:nx-1;

for j = 2:ny-1;

T(i,j) = ((T(i-1,j) + T_initial(i+1,j))/k*(dx^2)) + ((T_initial(i,j+1) + T(i,j-1))/k*(dy^2));

end

end

error_2 = max(max(abs(T_initial-T)))
T_initial = T;
Gauss_Sidel_iter = Gauss_Sidel_iter + 1

end
toc;
F = toc;
%Plotting the curve
figure(2)
[m n] = contourf(x,y,T);
clabel(m,n);
colormap(jet);
colorbar;
zlabel('x-axis');
ylabel('y-axis');
title(sprintf('Stady State heat Condition: Total iteration = %d, Time Taken = %fs',Gauss_Sidel_iter,F));
end

% initializing  Successive Over Relaxation Gauss-sidel Iteration
% Defining Relaxation constant (omega)[0 < omega < 2].
omega = 1.5;
tic;
SOR_Gauss_Sidel_iter = 3;
if iteration ==3;

error   = 9e9;
error_2 = error;
while (error_2 > Tolerance);

for i = 2:nx-1;

for j = 2:ny-1;

T(i,j) = (1-omega)*T_initial(i,j) +omega*(((T(i-1,j) + T_initial(i+1,j))/k*(dx^2)) + ((T_initial(i,j+1) + T(i,j-1))/k*(dy^2)));

end

end

error_2 = max(max(abs(T_initial-T)))
T_initial = T;
SOR_Gauss_Sidel_iter = SOR_Gauss_Sidel_iter + 1

end
toc;
F=toc;
%Plotting the curve
figure(2)
[m n] = contourf(x,y,T);
clabel(m,n);
colormap(jet);
colorbar;
zlabel('x-axis');
ylabel('y-axis');
title(sprintf('Steady State Heat conduction:total iteration = %d,Total Time  = %fs',SOR_Gauss_Sidel_iter,F));

end


RESULTS(PLOTS):   OBSERVATIONS:

From the results of steady state heat equation, we can figure out that using Gauss-jacobi method takes more iteration than any other method to determine the temperature distribution. Whereas, using successive over relaxation for Gauss-Seidel method using an optimum relaxation parameter "omega = 0.2" gives a huge margin in the rate of convergence from its predecessor. The final error obtained for following methods are:-

• Gauss-Jacobi (error)          =  9.9931e-11
• Gauss-Seidel (error)          =  9.9817e-11
• SOR Gauss -Seidel (error) =  9.9817e-11

### Meshing of an automotive hood for structural analysis Ankit Chakrabarty · 2019-12-27 04:07:31

Title:- Geometry cleanup and meshing of an automotive hood for structural analysis   Objective:  Perform geometry cleanup of automotive hood (both inner and outer panel). Perform mid-surface extraction of both the panels and Hem both the surfaces . Perf Read more

### Solving A Quasi-1 Dimensional Supersonic Incompressible Isentropic Nozzle Flow Using Maccormack Explict Method In MATLAB Ankit Chakrabarty · 2019-10-14 08:26:31

Abstract: This project presents the Modelling of the  behaviour of  fluid and its properties when it passes through a convergent-divergent nozzle using a Finite Difference Method. However, to simplify the understanding of the dynamics of the flow in the nozz Read more

### Case Study: Thermodynamic Performance Analysis Of a Twin-scroll Turbocharger vs Single-Scroll Turbocharger Using GT-Power. Ankit Chakrabarty · 2019-09-22 18:00:53

Turbocharger: Literature Review: A turbocharger is mechanical rotating device which forcebly induces air at high pressure into the internal combustion engine, in turn increases the power and effeciency (thermal and volumetric effeciency) of the engine. A turb Read more

### Effect Of Spectral Radius And Diagonal Magnification In The Convergence Rate Of An Iterative Solver For A Linear System Using MATLAB. Ankit Chakrabarty · 2019-09-20 17:47:52

(A) Problem defination:- Solve the following Linear equation using Gauss-Jacobi, Gauss-Seidal and SOR iterative method and understand the following parameters:- Spectral radius and the impact of spectral radius on the convergence of the solution. Effect of d Read more

### CASE STUDY:- A Thermodynamic Analysis Of A Tractor EngineFor Use In India With Given Engine Performance Targets. Ankit Chakrabarty · 2019-09-09 06:55:40

CASE STUDY:-  A Thermodynamic Analysis Of A Tractor Engine For Use In India With Given Engine Performance Targets. A Brief Literature Review: A tractor is an engineered vehicle specifically designed to deliver Power and traction/tractive effort to haul heavy equi Read more

### CASE STUDY: Effect of Fuel Burn Rates on Engine Performance Of A 6-Cyl-Turbocharged-DI Engine Using A Non-predictive Combustion model (DI WIEBE) Ankit Chakrabarty · 2019-08-23 20:03:11

CASE STUDY: Effect of Burn Rate on engine performance Of A 6-Cyl-Turbocharged-DI Engine Using A Non-predictive Combustion model (DI WIEBE) Abstract: Analytical functions approximating the burn rate in internal combustion engines are useful and cost-effective to Read more

### Solving 1-D Linear Convection Using First-Order Backward Difference And Forward Difference Numerical Schemes Using MATLAB Ankit Chakrabarty · 2019-08-12 10:58:54

Solving 1-D Linear Convection Using First-Order Backward Difference And Forward Difference method Using MATLAB LITERATURE REVIEW: WHAT IS CONVECTION? Convection is the sum of bulk transport of the  fluid  and brownian/ osmotic dispersion of fluid constituen Read more

### Computation Of 4th order Approximation For 2nd Order Derivative function Using Taylor Table In MATLAB Ankit Chakrabarty · 2019-08-07 17:51:14

Deriving the 4th order approximation for 2nd order derivative of a function using Taylor table for centered, skewed right-sided stencil and skewed left_sided stencil and compare the error results for the same using MATLAB.   Given function, f(x)    =& Read more

### Single Cylinder Spark Plug Engine Calibration at 3600 rpm for Higher Load Applications And Lower NOx And Other Emmisions using GT-Power Ankit Chakrabarty · 2019-07-25 10:55:20

Objective:- (1) Determine the engine parameters running at 1800 rpm and record the same. (2) Further increase the power output obtained at 3600 rpm by 10%.   Model:-  The GT model Of the 1-cylinder engine is shown below:- Observations (@1800rpm):-  Read more