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

Abstract

Temperature plot on a 2D plate is performed by using MATLAB. Steady state condition, as well as the unsteady state, is described. The governing equations are implemented and Jacobi, Gauss-Seidel and SOR iteration process is used to converging the final results.

Problem description

A 2D pyrolytic graphite with thermal diffusivity of 0.001 m2/s is subjected to constant temperature surfaces at all four sides. Boundary conditions conditions are left side : 400K, Topside : 600K, Right side : 800K, Bottom side : 900K. The temperature developed in the plate based on the time and the steady-state condition of the plate has to be analyzed. The effect of implicit and explicit methods of describing the conduction equations needs to be used.

Governing equation

2D conduction equation in cartesion coordinates is shown below

`(∂T/∂t)=α((∂^2 T)/(∂x^2 )+(∂^2 T)/(∂y^2 ))` ...........1

`α((∂^2 T)/(∂x^2 )+(∂^2 T)/(∂y^2 ))=0` ...........2

`(∂T)/(∂t)=(T_(i,j)^(n+1)-T_(i,j)^n)/dt` .................3

`(∂^2 T)/(∂x^2) =(T_(i+1,j)-2*T_(i,j)+T_(i-1,j))/(dx)^2` .............4

`(∂^2 T)/(∂y^2) =(T_(i,j+1)-2*T_(i,j)+T_(i,j-1))/dy^2 ...........5

Substituting the above equations 3,4,5 in 1 we get

`(T_(i,j)^(n+1)-T_(i,j)^n)/dt=α((T_(i+1,j)-2*T_(i,j)+T_(i-1,j))/dx^2 +(T_(i,j+1)-2*T_(i,j)+T_(i,j-1))/dy^2 )` ...............................6

Rearranging

`T_(i,j)^(n+1)=T_(i,j)^n+α*dt((T_(i+1,j)-2*T_(i,j)+T_(i-1,j))/dx^2 +(T_(i,j+1)-2*T_(i,j)+T_(i,j-1))/dy^2 )………….7`

`K_x=(α*dt)/dx^2 ……….8`

`K_y=(α*dt)/dy^2`..............9

Substituting in the above equation

`T_(i,j)^(n+1)=T_(i,j)^n+K_x (T_(i+1,j)-2*T_(i,j)+T_(i-1,j) )+K_y (T_(i,j+1)-2*T_(i,j)+T_(i,j-1) )`

Steady state 2d conduction equation coding

clear all
close all
clc
%% steady state heat condition 
n = 11;
T = zeros(n,n); % initiai temperature conditions
T(:,1) = 400; % bondary condition left
T(:,n) = 800; % bondary condition right
T(1,:) = 600; % bondary condition top
T(n,:) = 900; % bondary condition botom
T_old1 = T;T_old2 = T;T_old3 = T;
T1 = T;T2 = T;T3 = T;
XX = linspace(0,1,n);
YY = XX;
[x,y] = meshgrid(XX,YY);
% iteration input jacobi = 1, Gauss =2, Successive = 3
for method = 1:3
%% jacobi iteration
if method ==1
error1 = 9e9;
jcobi_iter = 0;
while error1 > 1e-4
    jcobi_iter =jcobi_iter + 1;
for j = 2:(n-1)   
for i = 2:(n-1)
    T1(i,j) = (T_old1(i-1,j)+T_old1(i+1,j)+T_old1(i,j-1)+T_old1(i,j+1))/4;
end
end
error1 = max(max(abs(T1-T_old1)));
T_old1 = T1;
subplot(3,3,1)
surf(x,y,T_old1)
title(['jcobi no iterations = ' num2str(jcobi_iter)])
xlabel('location (m)')
ylabel('location (m)')
zlabel('Temperature (C)')
colorbar
subplot(3,3,2)
pcolor(x,y,T_old1)
xlabel('location (m)')
ylabel('location (m)')
zlabel('Temperature (C)') 
subplot(3,3,3)
contour(x,y,T_old1)
xlabel('location (m)')
ylabel('location (m)')
zlabel('Temperature (C)') 
pause(0.01)
end
%% Gauss iteration
elseif method ==2
error2 = 9e9;
gauss_iter = 0;
while error2 > 1e-4
    gauss_iter =gauss_iter + 1;
for j = 2:(n-1)  
for i = 2:(n-1)
    T2(i,j) =(T2(i-1,j)+T_old2(i+1,j)+T2(i,j-1)+T_old2(i,j+1))/4;
end
end
error2 = max(max(abs(T2-T_old2)));
T_old2 = T2;
subplot(3,3,4)
surf(x,y,T_old2)
title(['Gauss no iterations = ' num2str(gauss_iter)])
xlabel('location (m)')
ylabel('location (m)')
zlabel('Temperature (C)')
colorbar
subplot(3,3,5)
pcolor(x,y,T_old2)
xlabel('location (m)')
ylabel('location (m)')
zlabel('Temperature (C)')
subplot(3,3,6)
contour(x,y,T_old2)
xlabel('location (m)')
ylabel('location (m)')
zlabel('Temperature (C)')
pause(0.01)
end
%% succesive over relaxation
else
c = 2*cos(pi/n);    
w = 4/(2+sqrt(4-c^2));% relaxation parameter
succover_iter = 0;
error3 = 9e9;
while error3 > 1e-4
    succover_iter =succover_iter + 1;
for j = 2:(n-1)  
for i = 2:(n-1)
    T3(i,j) =(1-w)*T3(i,j)+w*(T3(i-1,j)+T_old3(i+1,j)+T3(i,j-1)+T_old3(i,j+1))/4;
end
end
error3 = max(max(abs(T3-T_old3)));
T_old3 = T3;
subplot(3,3,7)
surf(x,y,T_old3)
title(['succ no iterations = ' num2str(succover_iter)])
xlabel('location (m)')
ylabel('location (m)')
zlabel('Temperature (C)')
colorbar
subplot(3,3,8)
pcolor(x,y,T_old3)
xlabel('location (m)')
ylabel('location (m)')
zlabel('Temperature (C)')
subplot(3,3,9)
contour(x,y,T_old3)
xlabel('location (m)')
ylabel('location (m)')
zlabel('Temperature (C)')
pause(0.01)
end
end
end

Transient explicit coding

clear all
close all
clc
% steady state heat condition 
n = 11;
T = zeros(n,n); % initiai temperature conditions
T(:,1) = 400; % bondary condition left
T(:,n) = 800; % bondary condition right
T(1,:) = 600; % bondary condition top
T(n,:) = 900; % bondary condition botom
T_old1 = T;
T1 = T;
T_left = [];
T_right = [];
T_top = [];
T_bottom = [];
XX = linspace(0,1,n);
YY = XX;
[x,y] = meshgrid(XX,YY);
dt = 1;        % time step size 
alpha = 0.001;  % thermal diffusivity
dx = XX(2)-XX(1);
T_centre = [];
index = 0;
duration = [];
time = 820;  % total time for which the analysis is performed
for t = 1:time/dt
    index = index + 1;
    duration = [duration index];
for j = 2:(n-1)   
for i = 2:(n-1)
    T1(i,j) = T_old1(i,j)+(alpha*dt/(dx^2))*(T_old1(i,j-1)+T_old1(i,j+1)+T_old1(i-1,j)+T_old1(i+1,j)-4*T_old1(i,j));
end
end 
T_old1 = T1;
subplot(1,2,1)
surf(x,y,T_old1)
title([ 'time = ' num2str(index) ' seconds'])
xlabel('location (m)')
ylabel('location (m)')
zlabel('Temperature (C)')
pause(0.001)
T_centre = [T_centre T1(6,6)];
T_left = [T_left T1(6,3)];
T_right = [T_right T1(6,9)];
T_top = [T_top T1(3,6)];
T_bottom = [T_bottom T1(9,6)];
subplot(1,2,2)
plot(duration,T_centre,duration,T_left,duration,T_right,duration,T_top,duration,T_bottom)
xlabel('time (s)')
ylabel('Temperature (C)')
title([ 'temperature at center of the plate is = ' num2str(T1(6,6)) 'C'])
legend('T1(6,6)','T1(6,3)','T1(6,9)','T1(3,6)','T1(9,6)')
pause(0.001)
end

Transient implicit coding

clear all
close all
clc
% steady state heat condition 
n = 11;
T = zeros(n,n); % initiai temperature conditions
T(:,1) = 400; % bondary condition left
T(:,n) = 800; % bondary condition right
T(1,:) = 600; % bondary condition top
T(n,:) = 900; % bondary condition botom
T_old1 = T;T_old2 = T;T_old3 = T;
T1 = T;T2 = T;T3 = T;
T_1= T; T_2= T; T_3= T;
XX = linspace(0,1,n);
YY = XX;
[x,y] = meshgrid(XX,YY);
dt = 10;        % time step size 
alpha = 0.001;  % thermal diffusivity
dx = XX(2)-XX(1);
time = 200;  % total time for which the analysis is performed
k = alpha*dt/(dx^2);
%% iteration input jacobi = 1, Gauss = 2, Successive = 3
for method = 1:3
if method == 1
% jacobi teration method
jacobi_iter = 0;
for t = 1:time/dt
    error1 = 9e9;
while error1 > 1e-4
   jacobi_iter = jacobi_iter +1; 
for j = 2:(n-1)   
for i = 2:(n-1)
    T1(i,j) = T_1(i,j)/(1+4*k)+(k/(1+4*k))*(T_old1(i,j-1)+T_old1(i,j+1)+T_old1(i-1,j)+T_old1(i+1,j));
end
end
error1 = max(max(abs(T1-T_old1))); 
T_old1 = T1;
subplot(3,1,1)
surf(x,y,T_1)
title(['jacobi no iterations = ' num2str(jacobi_iter) ' & time = ' num2str(time) ' seconds'])
xlabel('location (m)')
ylabel('location (m)')
zlabel('Temperature (C)')
pause(0.03)
end
T_1 = T1;
end
elseif method == 2
% Gauss siedel iteration method
Gauss_iter = 0;
for t = 1:time/dt
    error2 = 9e9;
while error2 > 1e-4
    Gauss_iter = Gauss_iter +1;
for j = 2:(n-1)   
for i = 2:(n-1)
    T2(i,j) = T_2(i,j)/(1+4*k)+(k/(1+4*k))*(T2(i,j-1)+T_old2(i,j+1)+T2(i-1,j)+T_old2(i+1,j));
end
end
error2 = max(max(abs(T2-T_old2))); 
T_old2 = T2;
subplot(3,1,2)
surf(x,y,T_2)
title(['Gauss no iterations = ' num2str(Gauss_iter) ' & time = ' num2str(time) ' seconds'])
xlabel('location (m)')
ylabel('location (m)')
zlabel('Temperature (C)')
pause(0.03)
end
T_2 = T2;
end
else
% succesive over relaxation method
 w = 1.5;  % relaxation parameter
Succ_iter = 0;
for t = 1:time/dt
    error3 = 9e9;
while error3 > 1e-4
    Succ_iter = Succ_iter +1;
for j = 2:(n-1)   
for i = 2:(n-1)
    T3(i,j) =(1-w)*T3(i,j)+w*(T_3(i,j)/(1+4*k)+(k/(1+4*k))*(T3(i,j-1)+T_old3(i,j+1)+T3(i-1,j)+T_old3(i+1,j)));
end
end
error3 = max(max(abs(T3-T_old3))); 
T_old3 = T3;
subplot(3,1,3)
surf(x,y,T_3)
title([' Succ no iterations = ' num2str(Succ_iter) ' & time = ' num2str(time) ' seconds'])
xlabel('location (m)')
ylabel('location (m)')
zlabel('Temperature (C)')
pause(0.03)
end
T_3 = T3;
end    
end
end

Result and conclusion

Steady-state conditions are formed by using three iterative techniques. (Jacobi, Gauss Seidel and SOC)

Each node is located at 0.1m from both sides. Jacobi method converged after 264 iterations and gauss seidel method took 140 iterations. The successive overrelaxation runs faster and only took 35 iterations to get the convergence. The optimum value of relaxation constant in this case is 1.5 which gives quick convergence.

Implicit methods are found to be unconditionally stable and noted a very small change in temperature at nodes compared to the explicit method.

Explicit transient heat conduction results show the plate becomes the steady state in 810 seconds. After 810 seconds there is no change in temperature at any node is noticed. It is also found that the temperature profile at each node is different with respect to time. The magnitude of temperature at nodes increases and become steady in 810 seconds. The variation in temperature with time at five different locations is shown below. The steady state is achieved at all nodes at same time.

The contour plot of the graph shows the isothermal lines. The upper half the isotherms are more close and more heat is transferred because the isotherms are close. The bottom half of the plate involves less heat transfer. Isotherms are more close at plate corners.