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

Aim: Steady state analysis & Transient State Analysis

Solve the 2D heat conduction equation by using the point iterative techniques that were taught in the class. The Boundary conditions for the problem are as follows;

Top Boundary = 600 K

Bottom Boundary = 900 K

Left Boundary = 400 K

Right Boundary = 800 K

You will implement the following methods for solving implicit equations.

1. Jacobi

2. Gauss-seidel

3. Successive over-relaxation

Your absolute error criteria is 1e-4

Theory:

• Conduction is a mode of heat transfer, generally takes place in solids due to temperature differences associated with free electron and lattice vibration.
• A generalized conduction equation is given

• Thermal diffusivity is a thermophysical property of a material defined as the ratio of thermal conductivity of the material and its thermal capacity.
• Since the problem is 2D and there is no heat generation is considered, therefore those terms are neglected. So the equation represents the unsteady 2D heat conduction (parabolic).

For steady-state, time derivative is not considered, their fore below equation represents the 2D steady-state Heat conduction (elliptic).

MATLAB Code :

• while programming it is assumed that, the domain is the surface of area (1mX1m). Discretized into 10 nodes both horizontally and vertically.
• Constant values like alpha and time step 'data, by keeping the explicit condition in mind.
• relaxation factor omega is also assumed while programming only.

Inputs

Top boundary temperature = 600K

Bottom Boundary temperature = 900K

Left Boundary temperature = 400K

Right Boundary temperature = 800K

error = 1e-4

Steady-State 2D Heat conduction

• First inputs are defined into variables and the temperature array is created.
• Boundary conditions are introduced into the temperature array, and at corners, boundary values are averaged.
• Constants like alpha, omega, tolerance are defined to run the iteration.
• Input must be provided by the user, to choose the iterative method.
• The if-else loop is used to solve the equation using the selected(input from the user) iterative technique.
• while loop is used to run the iterative technique till the required accuracy is achieved.
• Two for loop is used to calculate the temperature values at every node at each iterative step.
• Finally, the temperature distribution graph is obtained using 'conturf' command.

clear all
close all
clc

% Domain and dicretization length
number_node = 10;
domain_length = 1;

%value of x and y at nodes in a discretization plane
x = linspace(0,domain_length,number_node);
y = linspace(0,domain_length,number_node);
dx = domain_length/(number_node-1);
dy = domain_length/(number_node-1);

% Let define a temperature value at all nodes with initial guess 300k
T = 300*ones(number_node,number_node);

% Initia boundary conditions at left and right
T(1:end,1) = 400;
T(1:end,end) = 800;

% Initia boundary conditions at bottom and top
T(end,1:end) = 900;
T(1,1:end) = 600;

%average values of T at corner boundaries
T(1,1) = 500;
T(end,1) = 650;
T(1,end) = 700;
T(end,end) = 850;

%copy of Temperature to compute to compute
T_old = T;

%constant
K = (2*(dx^2+dy^2))/(dx^2*dy^2);
omegha = 1.2;

% to enter while loop initialize error value randomly
error = 10;
tol = 1e-4;
no_iteration = 0;

method = input('Enter(1) for jacobi,(2) for gauss siediel,(3) for sos');
if method == 1
while(error>tol)
for i = 2:number_node-1
for j = 2:number_node-1
term1 = (T_old(i+1,j)+T_old(i-1,j))/dx^2;
term2 = (T_old(i,j+1)+T_old(i,j-1))/dy^2;
T(i,j) = (term1+term2)/K;
end
end
error = max(max(abs(T_old-T)));

% to upadte the previous velcity with new velocity
%got from the iteration
T_old = T;
no_iteration = no_iteration+1;

end
elseif method == 2
while(error>tol)
for i = 2:number_node-1
for j = 2:number_node-1
term1 = (T_old(i+1,j)+T(i-1,j))/dx^2;
term2 = (T_old(i,j+1)+T(i,j-1))/dy^2;
T(i,j) = (term1+term2)/K;
end
end
error = max(max(abs(T_old-T)));

% to upadte the previous velcity with new velocity
%got from the iteration
T_old = T;
no_iteration = no_iteration+1;
end

else
while(error>tol)
for i = 2:number_node-1
for j = 2:number_node-1
term1 = (T_old(i+1,j)+T(i-1,j))/dx^2;
term2 = (T_old(i,j+1)+T(i,j-1))/dy^2;
T(i,j) = ((1-omegha)*T_old(i,j))+(omegha*((term1+term2)/K));
end
end
error = max(max(abs(T_old-T)));

% to upadte the previous velcity with new velocity
%got from the iteration
T_old = T;
no_iteration = no_iteration+1;

end
end

% to visualize the result

[X,Y] = meshgrid(x,y);
contourf(X,Y,T)
colorbar
set(gca,'ydir','reverse')

Unsteady-State 2D Heat conduction

• Inputs are defined into variables and the temperature array is created.
• Boundary conditions are introduced into the temperature array, and at corners, boundary values are averaged.
• Constants like alpha, omega, tolerance are defined to run the iteration.
• Here a user must provide, one or two inputs.
• The first input is to whether solve the problem using an explicit approach or an implicit approach.
• If an implicit approach is selected, then the user must provide another input, to select the iterative technique.
• The if-else loop is used to solve the equation using the selected(input from the user) approach and iterative technique.
• For loop is used to run the iteration for a large number of times to reach steady-state conditions.
• while loop is used to run the iterative technique till the required accuracy is achieved.
• Two for loop is used to calculate the temperature values at every node at each iterative step.
• Finally, a temperature distribution graph is obtained using 'conturf' command.

clear all
close all
clc

% Domain and dicretization length
number_node = 10;
domain_length = 1;

%value of x and y at nodes in a discretization plane
x = linspace(0,domain_length,number_node);
y = linspace(0,domain_length,number_node);
dx = (domain_length/number_node-1);
dy = (domain_length/number_node-1);
dt = 1/10;

% Let define a temperature value at all nodes with initial guess 300k
T = 300*ones(number_node,number_node);

% Initia boundary conditions at left and right
T(1:end,1) = 400;
T(1:end,end) = 800;
% Initia boundary conditions at bottom and top
T(end,1:end) = 900;
T(1,1:end) = 600;
%average values of T at corner boundaries
T(1,1) = 500;
T(end,1) = 650;
T(1,end) = 700;
T(end,end) = 850;

%copy of Temperature to compute for previous space terms
T_old = T;

% copy of T for previous time_step
T_prev = T;

%constants
alpha = 1.5;
omegha = 1.15;
K1 = (alpha*dt)/(dx^2);
K2 = (alpha*dt)/(dy^2);
iteration = 0;

approach = input('Enter 1 for implicit or Enter 2 for explicit')

if approach == 1
method = input('Enter 1 for jacobi or 2 for gauss siedel or 3 for sos')
if method == 1
for nt = 1:1500
error = 10;
tol = 1e-4;
while(error>tol)
for i = 2:number_node-1
for j = 2:number_node-1
term1 = 1/(1+(2*K1)+(2*K2));
term2 = K1*(T_old(i+1,j)+T_old(i-1,j));
term3 = K2 *(T_old(i,j+1)+T_old(i,j-1));
T(i,j) = term1*(T_prev(i,j)+term2+term3);
end
end
iteration = iteration+1;
error = max(max((abs(T_old-T))));
T_old = T;
end
T_prev = T;
end
elseif method ==2
for nt = 1:1500
error = 10;
tol = 1e-4;
while(error>tol)
for i = 2:number_node-1
for j = 2:number_node-1
term1 = 1/(1+(2*K1)+(2*K2));
term2 = K1*(T_old(i+1,j)+T(i-1,j));
term3 = K2 *(T_old(i,j+1)+T(i,j-1));
T(i,j) = term1*(T_prev(i,j)+term2+term3);
end
end
iteration = iteration+1;
error = max(max(abs(T_old-T)));
T_old = T;

end
T_prev = T;
end
else
for nt = 1:1500
error = 10;
tol = 1e-4;
while(error>tol)
for i = 2:number_node-1
for j = 2:number_node-1
term1 = 1/(1+(2*K1)+(2*K2));
term2 = K1*(T_old(i+1,j)+T(i-1,j));
term3 = K2 *(T_old(i,j+1)+T(i,j-1));
term4 = term1*(T_prev(i,j)+term2+term3);
T(i,j)=(1-omegha)*T_old(i,j)+(omegha*term4);
end
end
iteration = iteration+1;
error = max(max(abs(T_old-T)));
T_old = T;
end
T_prev = T;
end
end
else
for nt = 1:1500
for i = 2:number_node-1
for j = 2:number_node-1
term1 = K1*(T_prev(i+1,j)-2*T_prev(i,j)+T_prev(i-1,j));
term2 = K2*(T_prev(i,j+1)-2*T_prev(i,j)+T_prev(i,j-1));
T(i,j) = T_prev(i,j)+term1+term2;
end
end
% to upadte the previous timestep value with new timestep value
%got from the iteration
T_prev = T;
iteration = iteration+1;
end
end

% to visualize the result

[X,Y] = meshgrid(x,y);
figure(1)
contourf(X,Y,T)
colorbar
set(gca,'ydir','reverse')

Result:

• The program ran successfully.
• The temperature distribution plot obtained for both the steady and unsteady conditions is the same.
• Because in Transient analysis also, the solution approaches to steady-state.
• The top, bottom, left, and right side of the graphs representing their respective boundary conditions.

• The number of iterations taken to converge to the solution is compared among the 3 iterative techniques.
• In both steady and transient analysis, the number of iteration by the Successive over-relaxation method is less,
• The number of iteration taken by Successive over relaxaton<Number of iteration taken by gauss-siedel<Number of iteration taken by Jacobi.

Steady-State analysis:

1. Jacobi:

2. Gauss Siedel:

3. Successive over relaxation:

Transient analysis

• Out of a number of iterations,1500 iterations is due to time step

1. Jacobi:

2. Gauss Siedel:

3. Successive over-relaxation: