Steady and Unsteady 2D heat conduction Problem

Aim: To solve the steady and unsteady/Transient 2D Heat conduction problem using iterative methods(Jacobi,Gauss seidal, Successive Over Relaxation).

1. Assume that the domain is a unit square

2. Assume nx =ny[number of points along x direction is equal to the number of points along y direction]

3. Boundary conditions for steady and transient case

(i) Left : 400 K

(ii) Top : 600 K

(iii) Right : 800 K

(iv) Bottom : 900 K

Introduction:

Heat conduction is a mode of heat transfer, which transfers the heat from one substance to other substance when they are in physical contact or one end of the substance to other end of the same substance.

In solids,

Conduction is accomplished by lattice vibrations and drift of electrons.

In liquids,

Conduction is accomplished by collision of molecules and diffusion.

Convection and Radiation are the other two modes of heat transfer.

Generalized Heat conduction through 3 dimensional coordinates is

α({ ∂^2 T}/{∂x^2} + {∂^2 T}/{∂y^2} + {∂^2 T}/{∂z^2} ) +{dot q}/k = {∂T}/{∂t}

By eliminating the internal heat generation, finally we get

α({ ∂^2 T}/{∂x^2} + {∂^2 T}/{∂y^2} + {∂^2 T}/{∂z^2} ) = {∂T}/{∂t}

From the above 3D governing equation it is clear that the conduction is a time dependent, the transient nature arises from 3 dimensions and eventually move towards steady state(the temperature doesnot vary with time).

Derivation of numerical scheme:

From the generalized Energy equation in 3 dimensions we make the following assumptions

(i) convection and Radiation effects doesnot occurs

(ii) No internal heat generation or dissipation

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

Hence for 2 dimensional  condition we have:

α({ ∂^2 T}/{∂x^2} + {∂^2 T}/{∂y^2} ) = {∂T}/{∂t}

The second order of spatial terms is discretized by central differencing

{∂^2 T}/{∂x^2} = (T_(i-1,j) -2*(T)_(i,j)+T_(i+1,j) )/{Δx^2}

{∂^2 T}/{∂y^2} = (T_(i,j-1) -2*(T)_(i,j)+T_(i,j+1) )/(Δy^2)

The temporal terms is discretized by forward differencing

{∂T}/{∂t} = ((T^(n+1))_(i,j) - (T^n)_(i,j))/{Δt}

Hence we have:

α((T_(i-1,j) - 2*(T)_(i,j)+T_(i+1,j) )/{Δx^2} + (T_(i,j-1) - 2*(T)_(i,j)+T_(i,j+1) )/{Δy^2}) = {(T^(n+1))_(i,j) - (T^n)_(i,j)}/{Δt}

we define two constants K1 and K2 such that:

K_1 =α {∆t}/{(Δx)^2} and  K_2 =α{∆t}/{(Δy)^2}

Here  Δx = Δy

Hence the equation:

(T^(n+1))_(i,j) = (T^n)_(i,j) + K_1((T^(n+1))_(i-1,j) - 2(T^(n+1))_(i,j) + (T^(n+1))_(i+1,j)) + K_2((T^(n+1))_(i,j-1) - 2(T^(n+1))_(i,j) + (T^(n+1))_(i,j+1))

(T^(n+1))_(i,j) - K_1(T^(n+1))_(i-1,j) + 2 K_1(T^(n+1))_(i,j) - K_1 (T^(n+1)) _(i+1,j) - K_2(T^(n+1))_(i,j-1) + 2 K_2(T^(n+1))_(i,j) - K_2(T^(n+1))_(i,j+1) = (T^n)_(i,j)

(T^(n+1))_(i,j)(1+ 2 K_1 + 2 K_2) - K_1(T^(n+1))_(i-1,j) - K_1(T^(n+1))_(i+1,j) - K_2(T^(n+1)) _(i,j-1) - K_2(T^(n+1))_(i,j+1) = (T^n)_(i,j)

(T^(n+1))_(i,j) = ((T^n)_(i,j) + K_1(T^(n+1))_(i-1,j )+ K_1 (T^(n+1))_(i+1,j) + K_2(T^(n+1))_(i,j-1) + K_2(T^(n+1))_(i,j+1))/((1+ 2 K_1 + 2 K_2))

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

The other parameters are assumed as:

Thermal diffusivity,α = 1xx10^-1

Time step, dt = 0.01 seconds

number of grid points, nx =ny = 41

size of the grid, dx=dy= 0.0250

tolerance value = 10^-4

omega = 1.1

ambient temperature = 30^oC

SOLVING THE UNSTEADY STATE PROBLEM EXPLICITLY:

MATLAB CODE:

% Matlab code to solve 2D transient state Heat conduction equation explicitly

close all
clear all
clc

% Inputs
% number of grid points in x direction, nx
nx = 41
% number of grid points in y direction, ny
ny = nx
% thermal diffusivity, alpha
alpha = 1e-4
% Initializing time to start in seconds
t_start = 0
% Initializing time to end in seconds
t_end = 10000
% time step in seconds
dt = 0.01
% calculation
% x array along the length of the domain in x direction
x = linspace(0,1,nx)
% y array along the length of the domain in y direction
y = linspace(0,1,ny)
% Grid spacing in x direction, dx
dx = 1/(nx-1)
% Grid spacing in y direction, dx = 1/(nx-1)
dy = 1/(ny-1)

% Initializing the temperatures in x and y direction
T = 303*ones(nx,ny)
% Assigning Boundary conditions
T(1,:) = 600; % Top temperature
T(end,:) = 900; % Bottom temperature
T(:,1) = 400; % Left temperature
T(:,end) = 800;  % Right temperature
% Edge Refining
T(1,1) = (600+400)/2 ;
T(1,end) = (600+800)/2 ;
T(end,1) = (900+400)/2 ;
T(end,end) = (900+800)/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
error = 1e9
tol = 1e-4
for t = 1:10000
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);
xlabel('xx');
ylabel('yy');
title(sprintf('2D Unsteady state Heat transfer (Explicit Method) n Number of iterations:%d  (Time: %0.2fs)', counter-1,time_counter));

RESULTS & CONCLUSION FOR UNSTEADY STATE 2D HEAT CONDUCTION(EXPLICITLY): Here the timetaken for the iterations is very small because of the temperature at the time n+1 depends explicitly on the temperature at time n(i.e., Told). Here this is relatively simple and fast to compute the temperature but main drawback is stable solutions are obtained only if 0 < k{Delta t}/{Delta x^2} + k{Delta t}/{Delta y^2} le 0.5 . From the above plot it is clear that the temperature is a function of time.

SOLVING THE  UNSTEADY STATE PROBLEM IMPLICITLY:

MATLAB CODE:

%% Matlab program to solve 2D Unsteady state Heat conduction equation implicitly

close all
clear all
clc

% Inputs
% number of grid points in x direction, nx
nx = 41
% number of grid points in y direction, ny
ny = nx
% thermal diffusivity, alpha
alpha = 1e-4
% Initializing time to start in seconds
t_start = 0
% Initializing time to end in seconds
t_end = 10000
% time step in seconds
dt = 0.01

% calculation
% x array along the length of the domain in x direction
x = linspace(0,1,nx)
% y array along the length of the domain in y direction
y = linspace(0,1,ny)
% Grid spacing in x direction, dx
dx = 1/(nx-1)
% Grid spacing in y direction, dx = 1/(nx-1)
dy = 1/(ny-1)
% No. of time steps
nt = [t_start:dt:t_end]
% Initializing the temperatures in x and y direction
T = 303*ones(nx,ny)
% Assigning Boundary conditions
T(1,:) = 600; % Top temperature
T(end,:) = 900; % Bottom temperature
T(:,1) = 400; % Left temperature
T(:,end) = 800;  % Right temperature
% Edge Refining
T(1,1) = (600+400)/2 ;
T(1,end) = (600+800)/2 ;
T(end,1) = (900+400)/2 ;
T(end,end) = (900+800)/2;

% optimal value of over relaxation parameter (omega);
omega = 1.1
% calculation of temperature distribution
k1 =  (alpha*dt)/(dx^2);
k2 =  (alpha*dt)/(dy^2);
[xx, yy] = meshgrid(x,y);
Told = T;
T_prev_dt = Told;

counter = 1;
iterative_solver = input('Enter the number for the dedicated solver type: n 1 for jacobi n 2 for gauss seidel n 3 for SOR n')
error = 1e9
tol = 1e-4
tic
for k = 1:10000
% Jacobi Iteration scheme
if iterative_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)=((T_prev_dt(i,j)*term1)+(H*term2)+(V*term3));
end
end
error= max(max(abs(Told-T)))
Told = T
counter = counter + 1
end
end

% Gauss seidal scheme
if iterative_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)=((T_prev_dt(i,j)*term1)+(H*term2)+(V*term3));
end
end
error= max(max(abs(Told-T)))
Told = T
counter = counter + 1
end
end

% Successive over-relaxation scheme
if iterative_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)*((T_prev_dt(i,j)*term1)+(H*term2)+(V*term3)) + (1-omega)*(T_prev_dt(i,j))
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 = ["Jacobi","Gauss Seidel","Successful Over Relaxation"];
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(iterative_solver), k*dt, counter, execution_time));

RESULTS & CONCLUSION FOR UNSTEADY STATE 2D HEAT CONDUCTION(IMPLICITLY):

1. Jacobi Method 2. Gauss seidel Method 3. Successive over Relaxation Here the timetaken for the iterations is very large because of the temperature at the time n+1 depends implicitly on the temperature at new time step n+1 at different nodes. Here this is relatively complicated and slow to compute the temperature From the above plots it is clear that the temperature is a function of time.

Derivation of numerical scheme:

From the generalized Energy equation in 3 dimensions we make the following assumptions

(i) convection and Radiation effects doesnot occurs

(ii) No internal heat generation or dissipation

(iii) For steady state, the Tne f(time) then

{∂T}/{∂t} =0

Hence for 2D the equation becomes

({ ∂^2 T}/{∂x^2} + {∂^2 T}/{∂y^2} ) = 0

The second order of spatial terms is discretized by central differencing

{∂^2 T}/{∂x^2} = (T_(i-1,j) -2*(T)_(i,j)+T_(i+1,j) )/{Δx^2}

{∂^2 T}/{∂y^2} = (T_(i,j-1) -2*(T)_(i,j)+T_(i,j+1) )/(Δy^2)

Hence we have

((T_(i-1,j) - 2*(T)_(i,j)+T_(i+1,j) )/{Δx^2} + (T_(i,j-1) - 2*(T)_(i,j)+T_(i,j+1) )/{Δy^2}) = 0

Also Δx = Δy

The other parameters are assumed as:

Time step, dt = 0.01 seconds

number of grid points, nx =ny = 41

size of the grid, dx=dy= 0.0250

tolerance value = 10^-4

SOLVING THE STEADY STATE PROBLEM IMPLICITLY:

MATLAB CODE:

%% Matlab program to solve 2D steady state Heat conduction equation implicitly

close all
clear all
clc

% Inputs
% number of grid points in x direction, nx
nx = 41
% number of grid points in y direction, ny
ny = nx
% thermal diffusivity, alpha
alpha = 1e-4
% Initializing time to start in seconds
t_start = 0
% Initializing time to end in seconds
t_end = 10000
% time step in seconds
dt = 0.01

% calculation
% x array along the length of the domain in x direction
x = linspace(0,1,nx)
% y array along the length of the domain in y direction
y = linspace(0,1,ny)
% Grid spacing in x direction, dx
dx = 1/(nx-1)
% Grid spacing in y direction, dx = 1/(nx-1)
dy = 1/(ny-1)
% No. of time steps
nt = [t_start:dt:t_end]
% Initializing the temperatures in x and y direction
T = 303*ones(nx,ny)
% Assigning Boundary conditions
T(1,:) = 600; % Top temperature
T(end,:) = 900; % Bottom temperature
T(:,1) = 400; % Left temperature
T(:,end) = 800;  % Right temperature
% Edge Refining
T(1,1) = (600+400)/2 ;
T(1,end) = (600+800)/2 ;
T(end,1) = (900+400)/2 ;
T(end,end) = (900+800)/2;

% optimal value of over relaxation parameter (omega);
omega = 1.1
% calculation of temperature distribution
k = 2*(1/(dx^2)+ 1/(dx^2));
[xx, yy] = meshgrid(x,y);
Told = T;
T_prev_dt = Told;

iterative_solver = input('Enter the number for the dedicated solver type: n 1 for jacobi n 2 for gauss seidel n 3 for SOR n')
error = 1e9
tol = 1e-4
tic
% Jacobi Iteration scheme
if iterative_solver == 1
jacobi_iter = 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
jacobi_iter = jacobi_iter+1;
end
end

% Gauss seidal scheme
if iterative_solver == 2
gs_iter = 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)= term1+term2;
end
end
error= max(max(abs(Told-T)))
Told = T
gs_iter = gs_iter+1;
end
end

% Successive Over Relaxation scheme
if iterative_solver == 3
sor_iter = 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
sor_iter = sor_iter+1;
end
end
time_required = toc;

figure(1)
[C,h]=contourf(xx,yy,T);
colorbar;
colormap(jet);
clabel(C,h);
if iterative_solver ==1
title(sprintf('Temperature profile for steady state Heat conduction using Jacobi Method n No. of Iterations:%d n computation time:%f s',jacobi_iter, time_required));
else
if iterative_solver == 2
title(sprintf('Temperature profile for steady state Heat conduction using Gauss siedel Method n No. of Iterations:%d n computation time:%f s',gs_iter, time_required));
else
if iterative_solver == 3
title(sprintf('Temperature profile for steady state Heat conduction using Successive over Relaxation Method n No. of Iterations:%d n computation time:%f s',sor_iter, time_required));
end
end
end
xlabel('xx');
ylabel('yy');

RESULTS AND CONCLUSION FOR STEADY STATE 2D HEAT CONDUCTION:

1. Jacobi Method 2. Gauss seidel Method 3. Successive Over Relaxation Method Jacobi method converged after 3187 iterations and Gauss-seidel method took 1711 iterations to converge where as successive over relaxation runs faster and took 1427 iterations to converge. The optimum value of overrelaxation constant omega is taken as 1.1 for Successive over Relaxation method, which gives quick convergence.

The contour plot show  the isothermal lines, If the isotherms are more close and more heat is transferred, if the isotherms are less close then less heat is transferred.

Parametric Study of Gate Valve Shravankumar Nagapuri · 2019-10-22 03:14:10

Objectives: To perform a parametric study on the gate valve simulation. To obtain the mass flow rates at the outlet for 5 design points. To show the cut section view for different design points, and also to show the gate disc lift and fluid volume extraction. To s Read more

Symmetry vs Wedge vs HP equation Shravankumar Nagapuri · 2019-10-09 22:50:52

Objectives: To Write a Matlab program that can generate the computational blockMesh file  automatically for a wedge angle of 3 degrees and compare the results obtained for Symmetry BC and Wedge BC. To Write a Matlab program that takes an angle as input and gene Read more

Combustion simulation on the combustor model Shravankumar Nagapuri · 2019-10-05 16:23:22

1Q. Briefly explain about the possible types of combustion simulations in FLUENT. A:  Combustion or burning is a high temperature exothermic chemical reaction between a fuel and an oxidant accompanied by the production of heat, light and unburnt gases in the fo Read more

Cyclone separator challenge Shravankumar Nagapuri · 2019-10-01 17:50:02

Cyclone Separator: Principle and Working: A high speed rotating (air)flow is established within a cylindrical or conical container called a cyclone. Air flows in a spiral pattern, beginning at the top (wide end) of the cyclone and ending at the bottom (narrow) end bef Read more

Gearbox sloshing effect Shravankumar Nagapuri · 2019-09-28 15:37:58

Objective: To analyse the flow pattern of the fluid inside the gearbox, for two different clearances of the same geometry and each geometry flow is analyzed by two different fluids (Water and Oil) as lubricants. Gearbox sloshing effect:  Slosh refers Read more

Simulation of Flow through a pipe in OpenFoam Shravankumar Nagapuri · 2019-09-26 19:10:15

Objectives: To write a program in Matlab that can generate the computational mesh automatically for any wedge angle and Grading scheme. To calculate, length of the pipe using the entry length formula for laminar flow through a pipe To show that entry length is suf Read more

Conjugate Heat Transfer CHT Analysis on a graphics card Shravankumar Nagapuri · 2019-09-20 16:08:56

Objectives: To find out the maximum temperature attained by the processor. To Prove that the simulation has achieved convergence with appropriate images and plots. To Find out the heat transfer coefficient at appropriate areas of the model. To Identify potential h Read more

Finite Volume Method FVM Literature Review Shravankumar Nagapuri · 2019-09-17 13:14:28

Finite volume Method: The finite volume method (FVM) is a method for representing and evaluating partial differential equations in the form of algebraic equations. Similar to the finite difference method (FDM) or finite element method (FEM), values are calculated at Read more

Rayleigh Taylor instability Shravankumar Nagapuri · 2019-09-14 17:48:30

Q1. What are some practical CFD models that have been based on the mathematical analysis of Rayleigh Taylor waves? In your own words, explain how these mathematical models have been adapted for CFD calculations  Kelvin–Helmholtz instability: The Kelvin Read more

BlockMesh Drill down challenge Shravankumar Nagapuri · 2019-09-11 18:17:33

Objectives: To Generate the BlockMesh file for the given Geometry using OpenFOAM and use the icoFoam solver to simulate the flow through a backward-facing step. To create multiple meshes and will be comparing the results obtained from each mesh. To show the velocit Read more