Simulation of a 2D Heat Conduction problem in steady and unsteady/transient forms using iterative methods.

Project Objectives:

1. Solving the 2 Dimensional Heat conduction equation in the generalized form using various iterative techniques:

i. Explicit Solver (for unsteady state)

(a). Gauss Jacobi method

(b). Gauss Siedel method

(c). Successive Over Relaxation method

2. Analysis of steady state solution and unsteady state solution and comparison of them:

3. Stability Analysis of the unsteady state solution: ______________________________________________________________________

[Part: 1/3]: SOLVING THE 2D HEAT CONDUCTION EQUATION:

Introduction: Thermal conduction is the transfer of heat by molecular collisions of particles and movement of electrons with high kinetic energy within a body. It is an important mode of heat transfer that occurs between two bodies in contact. Conduction requires an temperature gradient between two bodies to take place, the other important modes of heat transfer are convection and radiation.

Conduction is essentially a time-dependant phenomena, the transient nature arises from the governing energy equation in 3 dimensions and transient conduction eventually moves towards steady state. (i.e. when the temperature gradients that have been established within a body do not change with time)

Problem specifications:

i. Assume that the domain is an unit square

ii. Assume Nx = Ny (Number of gridpoints is same in both x and y directions or the mesh is a square mesh)

iii. Boundary Conditions: (Drichlet type)

1. Left wall = 400K

2. Top wall = 600K

3. Right wall = 800K

4. Bottom wall = 900K

The Plate is subjected to the boundary conditions at time t = 0. The Temperature of the plate is considered to be the ambient temperature i.e. T = 298K

Derivation of the Numerical Scheme:

From the generalized Energy equation in 3 Dimensions we make the following assumptions:

i. No convection occurs

ii. No internal heat generation or dissipation

iii. Constant properties of system i.e. alpha != f(t,x,y)

Hence for a 2 Dimensional condition we have:

(del^2T)/(delx^2)+(del^2T)/(dely^2) = alpha(delT)/(delt)

We Discretize the spatial terms by central differencing (O(h^2)) i.e.

(del^2T)/(delx^2) = (T_(i-1,j)-2T_(i,j)+T_(i+1,j))/(Deltax^2)

(del^2T)/(dely^2) = (T_(i,j-1)-2T_(i,j)+T_(i,j+1))/(Deltay^2)

We discretize the temporal terms by forward differencing (O(h)) i.e.

(delT)/(delt) = ((T^(n+1))_(i,j)-(T^(n))_(i,j))/(Deltat)

Hence we have:

(T_(i-1,j)-2T_(i,j)+T_(i+1,j))^(n+1)/(Deltax^2) + (T_(i,j-1)-2T_(i,j)+T_(i,j+1))^(n+1)/(Deltay^2) = ((T^(n+1))_(i,j)-(T^(n))_(i,j))/(Deltat)

We define two constants K1 and K2 such that:

K1 = alpha(Deltat)/((Deltax)^2)

K2 = alpha(Deltat)/((Deltay)^2)

Here Deltax = Deltay

Hence the equation:

(T^(n+1))_(i,j)=(T^(n))_(i,j)+K1((T^(n+1))_(i-1,j)-2(T^(n+1))_(i,j)+(T^(n+1))_(i+1,j))+K2((T^(n+1))_(i,j-1)-2(T^(n+1))_(i,j)+(T^(n+1))_(i,j+1))

(T^(n+1))_(i,j)=((T^(n))_(i,j)+(K1(T_(i-1,j)+T_(i+1,j))+K2(T_(i,j-1)+T_(i,j+1)))^(n+1))/(1+2K1+2K2)

Here n denotes the time index and (i,j) denotes the space index

Other conditions are specified for the problem:

Thermal Diffusivity alpha = 10^-4

Time-step = 0.01 seconds

Number of grid-points = 51

dx = dy = 0.02

tolerance value = 10^-4

1. Explicit Iteration Scheme:

(A) MATLAB PROGRAM:

%% MATLAB Program to generate temperature profile of 2D Transient State Heat Conduction Equation explicitly:
clear all;
close all;
clc;

%% Initializing Variables:
Nx = 51;
Ny = Nx;
x = linspace(0,1,Nx);
y = linspace(0,1,Ny);
dx = 1/(Nx-1);
dy = 1/(Ny-1);
alpha = 1e-4;
t_start = 0;
t_end = 1/alpha;
dt = 0.01;
error = 6e9;
tol = 1e-4;

%% Initializing the Temperature grid:
T = 298*ones(Nx,Ny);
T(1,:) = 900; % Bottom Temperature
T(end,:) = 600; % Top Temperature
T(:,1) = 400; % Left Temperature
T(:,end) = 800; % Right Temperature
% 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/(dx^2);

%% Calculation of Temperature Distribution explicitly using iterative solvers:

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


(B) RESULTS:

2. Implicit Iteration Scheme:

(A) MATLAB CODE:

%% MATLAB Program to solve the unsteady 2D Heat Conduction equation Implictly:
clear all;
close all;
clc;

%% Initializing Variables:
Nx = 51;
Ny = Nx;
x = linspace(0,1,Nx);
y = linspace(0,1,Ny);
dx = 1/(Nx-1);
dy = 1/(Ny-1);
alpha = 1e-4;
t_start = 0;
t_end = 1/alpha;
dt = 0.01;
Nt = [t_start:dt*50:t_end];
error = 6e9;
tol = 1e-4;

%% Initializing the Temperature grid:
T = 298*ones(Nx,Ny);
T(1,:) = 900; % Bottom Temperature
T(end,:) = 600; % Top Temperature
T(:,1) = 400; % Left Temperature
T(:,end) = 800; % Right Temperature
% 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;

%% Calculation of optimal value of over-relaxation parameter (omega):
omega = 2/(1+sin(pi*dx));

%% Calculation of Temperature Distribution explicitly using iterative solvers:

k1 = (alpha*dt)/(dx^2);
k2 = (alpha*dt)/(dy^2);
[xx,yy] = meshgrid(x,y);
Told = T;
T_prev_dt = Told;
counter = 1;
solver = input('Enter the number for the dedicated solver type: n 1 for Gauss Jacobi n 2 for Gauss Siedel n 3 for SOR n');
tic;
for k = 1:20000
%% Gauss Jacobi implicit Iteration Scheme:
if solver == 1
term1 = (1+2*k1+2*k2)^(-1);
term2 = k1*term1;
term3 = k2*term1;
error = 1e9;
while (error > tol)
for i = 2:Nx-1
for j = 2:Ny-1
H = (Told(i-1,j)+Told(i+1,j));
V = (Told(i,j-1)+Told(i,j+1));
T(i,j) = term1*T_prev_dt(i,j)+term2*H+term3*V;
end
end
error = max(max(abs(Told-T)));
Told = T;
counter = counter + 1;
end
end
if solver == 2
term1 = (1+2*k1+2*k2)^(-1);
term2 = k1*term1;
term3 = k2*term1;
error = 1e9;
while (error > tol)
for i = 2:Nx-1
for j = 2:Ny-1
H = (T(i-1,j)+Told(i+1,j));
V = (T(i,j-1)+Told(i,j+1));
T(i,j) = term1*T_prev_dt(i,j)+term2*H+term3*V;
end
end
error = max(max(abs(Told-T)));
Told = T;
counter = counter + 1;
end
end
if solver == 3
term1 = (1+2*k1+2*k2)^(-1);
term2 = k1*term1;
term3 = k2*term1;
error = 1e9;
while (error > tol)
for i = 2:Nx-1
for j = 2:Ny-1
%H = (T(i-1,j)+Told(i+1,j));
%V = (T(i,j-1)+Told(i,j+1));
%T(i,j) = (omega)*(term1*T_prev_dt(i,j)+term2*H+term3*V) + (1-omega)*T_prev_dt(i,j);
H=(omega*T(i-1,j)+(1-omega)*Told(i-1,j)+Told(i+1,j));
V=(omega*T(i,j-1)+(1-omega)*Told(i,j-1)+Told(i,j+1));
T(i,j)=T_prev_dt(i,j)*term1+(term2*H+term3*V);
end
end
error = max(max(abs(Told-T)));
Told = T;
counter = counter + 1;
end
end
T_prev_dt = T;
end
execution_time = toc;
figure(1);
str_array = ["Gauss-Jacobi","Gauss-Siedel","SOR"];
formatSpec = "Unsteady state 2D Heat Transfer by implicit: (%s Method) n Time: %0.3fs n Iteration Counter: %d n Computation time: %0.3fs";

[C,h] = contourf(xx,yy,T);
colorbar;
colormap(jet);
clabel(C,h);
xlabel('xx');
ylabel('yy');
title(sprintf(formatSpec,str_array(solver),k*dt,counter,execution_time));


(B) RESULTS:

i. Gauss Jacobi Method:

ii. Gauss Siedel Method:

iii. Successive Over Relaxation Method:

Derivation of the numerical scheme:

From the generalized Energy equation in 3 Dimensions we make the following assumptions:

i. No convection occurs

ii. No internal heat generation or dissipation

iii. (delT)/(delt) = 0

Hence the equation becomes:

(del^2T)/(delx^2)+(del^2T)/(dely^2) = 0

We Discretize the spatial terms by central differencing (O(h^2)) i.e.

(del^2T)/(delx^2) = (T_(i-1,j)-2T_(i,j)+T_(i+1,j))/(Deltax^2)

(del^2T)/(dely^2) = (T_(i,j-1)-2T_(i,j)+T_(i,j+1))/(Deltay^2)

hence:(T_(i-1,j)-2T_(i,j)+T_(i+1,j))/(Deltax^2) + (T_(i,j-1)-2T_(i,j)+T_(i,j+1))/(Deltay^2) = 0

Also Deltax = Deltay

T_i = 1/4(T_(i-1,j)+T_(i+1,j)+T_(i,j-1)+T_(i,j+1))

Number of gridpoints Nx = Ny = 81

other problem specifications are same as before.

SOLVING THE SYSTEM IMPLICTLY:

(A) MATLAB PROGRAM:

%% MATLAB Program to generate the 2D temperature distribution at steady state:
clear all;
close all;
clc;

%% Initial Conditions:
Nx = 81;
Ny = Nx;
x = linspace(0,1,Nx);
y = linspace(0,1,Ny);
dx = 1/(Nx-1);
dy = 1/(Ny-1);

k = 2*(1/(dx^2)+1/(dy^2));
error = 1e8;
tol = 1e-4;

%% Initializing the temperature matrix:

T = 298*ones(Nx,Ny);
T(1,:) = 900; % Bottom Temperature
T(end,:) = 600; % Top Temperature
T(:,1) = 400; % Left Temperature
T(:,end) = 800; % Right Temperature
% 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;

%% Calculation of over-relaxation parameter Omega:

omega = 2/(1+sin(pi*dx));
%% Solving the 2D steady state system iteratively:
[xx,yy] = meshgrid(x,y);
Told = T;

solver = input('Enter the desired number: ')

tic;
%% Gauss Jacobi Iteration Scheme:
if solver == 1
iter_count_jac = 1;
while (error > tol)
for i = 2:Nx-1
for j = 2:Ny-1
term1 = (Told(i-1,j)+Told(i+1,j))/(k*(dx^2));
term2 = (Told(i,j-1)+Told(i,j+1))/(k*(dy^2));
T(i,j) = term1+term2;
end
end
error = max(max(abs(Told-T)));
Told = T;
iter_count_jac = iter_count_jac + 1;
end
end
%% Gauss Siedel Iteration Scheme:
if solver == 2
iter_count_gs = 1;
while (error > tol)
for i = 2:Nx-1
for j = 2:Ny-1
term1 = (T(i-1,j)+T(i+1,j))/(k*(dx^2));
term2 = (T(i,j-1)+T(i,j+1))/(k*(dy^2));
T(i,j) = term1+term2;
end
end
error = max(max(abs(Told-T)));
Told = T;
iter_count_gs = iter_count_gs + 1;
end
end
%% Successive Over Relaxation (SOR) Iteration Scheme:
if solver == 3
iter_count_sor = 1;
while (error > tol)
for i = 2:Nx-1
for j = 2:Ny-1
term1 = (T(i-1,j)+Told(i+1,j))/(k*(dx^2));
term2 = (T(i,j-1)+Told(i,j+1))/(k*(dy^2));
T(i,j) = omega*(term1+term2)-(omega-1)*Told(i,j);
end
end
error = max(max(abs(Told-T)));
Told = T;
iter_count_sor = iter_count_sor + 1;
end
end
time_required = toc;
figure(1)
[C,h] = contourf(xx,yy,T);
colorbar;
colormap(jet);
clabel(C,h);
if solver == 1
title(sprintf('Temperature profile of steady state heat conduction using Gauss Jacobi method n No. of Iterations: %d n Computation Time: %0.3f s',iter_count_jac,time_required));
else if solver == 2
title(sprintf('Temperature profile of steady state heat conduction using Gauss Siedel method n No. of Iterations: %d n Computation Time: %0.3f s',iter_count_gs,time_required));
else if solver == 3
title(sprintf('Temperature profile of steady state heat conduction using SOR method n No. of Iterations: %d n Computation Time: %0.3f s n Relaxation Parameter: %f',iter_count_sor,time_required,real(omega)));
end
end
end
xlabel('xx');
ylabel('yy');


(B) RESULTS:

i. Gauss Jacobi Method:

ii. Gauss Siedel Method:

iii. Successive Over Relaxation Method:

The Optimal Value of omega is given by the formula:

omega_(opt)=2/(1+sin(pidx)

### Pipe flow simulation in OpenFOAM (Part 2/2: Symmetry boundary conditions) Priyotosh Bairagya · 2018-10-08 11:27:26

LAMINAR INCOMPRESSIBLE FLOW SIMULATION THROUGH A PIPE IN OPENFOAM: [PART 2/2]: SYMMETRY BOUNDARY CONDITION: OBJECTIVES:1. Creating the Mesh-Script (blockMeshDict file) with symmetry boundary condition.2. Simulation Results and post processing the velocity profile at Read more

### Pipe flow simulation in OpenFOAM (Part 1/2: Wedge boundary conditions) Priyotosh Bairagya · 2018-10-07 05:05:54

LAMINAR INCOMPRESSIBLE FLOW SIMULATION THROUGH A PIPE IN OPENFOAM: [PART 1/2]: WEDGE BOUNDARY CONDITION: OBJECTIVES:1. Calculation of Pre-eliminary quantinites related to the flow.2. Creating the Mesh script (blockMeshDict file) for specifying the geometry and bounda Read more

### Flux Limiters and Interpolation Schemes in Finite Volume Method Priyotosh Bairagya · 2018-10-05 12:18:40

I. Interpolation Schemes in Finite Volume Method: The approximation of surface and volume integrals may require values of the variable at locations other than the computational nodes of the CV. Values at these locations are obtained using interpolation formulae. Some o Read more

### Analysis of Solution Stability in a 2D Heat conduction problem Priyotosh Bairagya · 2018-09-20 22:56:03

ANALYSIS OF NUMERICAL STABILITY OF VARIOUS ITERATIVE SOLVERS FOR TRANSIENT 2D HEAT CONDUCTION: [Part: 3/3] INTRODUCTION: The criterion of stability of a numerical scheme is determined by the way the errors propagate while the solution moves from one time-step to th Read more

### BlockMesh Analysis of a Backward Facing Step Priyotosh Bairagya · 2018-09-20 17:08:11

MESH GENERATION AND ANALYSIS USING BLOCKMESH FOR FLOW OVER A BACKWARD FACING STEP: The purpose of the following project is to generate the geometry for a variation of the incompressible cavity flow problem in OpenFOAM. For this purpose we have modified the lid-driven c Read more

### Analysis of Steady and Unsteady State solutions of a 2D Heat conduction problem Priyotosh Bairagya · 2018-09-16 10:37:07

ANALYSIS OF VARIOUS ITERATIVE SCHEMES FOR THE SOLUTION OF A 2D HEAT CONDUCTION PROBLEM: [PART: 2/3] In the previous part we had explained the problem statement and the MATLAB Program in detail. In this Part we are going to explain the outputs from the 2D Heat Conduct Read more

### Flow simulation of a 1D Super-Sonic nozzle using the Mac-Cormack method Priyotosh Bairagya · 2018-09-02 12:55:50

NUMERICAL SOLUTION OF 1D SUPERSONIC NOZZLE FLOW SIMULATION BY MAC-CORMACK METHOD PROJECT OBJECTIVES: i. Numerical solution of the governing equations in both conservative and non-conservative forms. ii. Creating user defined functions for calculating the flow quan Read more

### Iterative solution of a system of linear equations and an analysis of spectral radius of a matrix Priyotosh Bairagya · 2018-08-31 10:13:18

UNDERSTANDING LINEAR SYSTEMS(ANALYSIS OF VARIOUS ITERATIVE SCHEMES TO SOLVE A SYSTEM OF LINEAR EQUATIONS TO FIND THE EIGEN VALUES AND SPECTRAL RADIUS) (A) PROBLEM STATEMENT: Given coefficient matrix: A = [[5,1,2],[-3,9,4],[1,2,-7]] Given Solution Matrix: X = [[x Read more

### MATLAB Program to solve the 1D linear wave equation Priyotosh Bairagya · 2018-08-18 22:51:02

WEEK 4: (Effect of Grid-Size on output for the solution of 1D linear wave equation) 1. Problem Setup: Given Partial Differential Equation: (delu)/(delt) + c(delu)/(delx) = 0 Numerical Discretization:  u_(i,n+1) = u_(i,n) + (cDeltat)/(Deltax)*(u_(i-1,n+1) Read more