Week 5.1 - Mid term project - Solving the steady and unsteady 2D heat conduction problem
AIM:
- To solve the 2D heat conduction problem for the steady and unsteady state using Jacobi, gauss seidel, and successive over-relaxation iterative techniques.
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).
- Difference equations for the steady-state and unsteady-state 2D conduction equation are derived to obtain the numerical solution.
| Steady State 2D heat conduction |  |
| Transient 2D heat conduction (explicit) |  |
| Transient 2D heat conduction (implicit) |  |
- Here iterative methods are used to solve the difference equation, Iterative methods are numerical techniques used to converge the solution(accurately) of a problem with some initial guess.
- Here Jacobi, gauss-seidel, and successive over-relaxation techniques are used.
MATLAB CODE:
- There is 2 Matlab code written separately for steady and unsteady conditions.
- But for both conditions inputs are the same
- 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 | 600 K |
| Bottom Boundary temperature | 900 K |
| Left Boundary temperature | 400 K |
| Right Boundary temperature | 800 K |
| 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.
Code:
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.
Code:
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 run 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
Steady-State analysis
- Jacobi
- Gauss-Siedel
- Successive over relaxation
Transient analysis
- Out of a number of iterations,1500 iterations is due to time step
- Jacobi
- Gauss-Siedel
- Successive over-relaxation
- Explicit
- 1500 iterations are due to time step
- According to CFL stability criteria, K1+K2<0.5.
