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

`(∂^2T)/(∂x^2) + (∂^2T)/(∂y^2) = (1/alpha)(∂T)/(∂t)` unsteady state`

`(∂^2T)/(∂x^2) + (∂^2T)/(∂y^2) = 0` steady state condition`

Descretization

descretization of a second order partial differential equation is performed by central differencing methods and is the order of second degree accuracy.

`(∂^2T)/(∂x^2) = (T(x-h) -2T(x)+T(x+h))/h^2`

`(∂^2T)/(∂y^2) = (T(y-h) -2T(y)+T(y+h))/h^2`

`(∂T)/(∂t) = (T^(n)-T)/(△t)`

Steady state conditions

Sine there is no temperature variation after reaching steady state conditions the convergence of the temperature is obtained bt ierative techniques like Jacobi, Gauss siedel and successive over relaxation. The descretization of the steady state condition is written as

`(T(x-h) -2T(x)+T(x+h))/h^2 + (T(y-h) -2T(y)+T(y+h))/h^2 = 0`

Temperature at each location is analyzed by

`T = (T(x-h)+T(x+h)+T(y-h)+T(y+h))/4`

ans the ieration is carried out by Jaccobi, Guass and succesive over relaxation methods.

Transient conduction

Temperature changes with time and reaches a steady state. Based on the discretization the conduction equation (diffusion) is manipulated as

`(T(x-h) -2T(x)+T(x+h))/h^2 + (T(y-h) -2T(y)+T(y+h))/h^2 = (1/alpha)(T^n-T)/(△t)`

This equation is iterated by both explicit as well as implicit methods. In explicit methods, the temperature at next time level at a point is calculated by based on the neighbouring temperature values. Increasing the time step value more than the maximum allowable limit, the result will oscillate unpredictably and diverges from the true results. The maximum possible value of '△t' calculated for the given problem is

`1-4*(alpha*(△t)/(△x^2)) >=0`

where △x = 0.1111, alpha = 0.001 From the above expression the maximum value of △t is approximated as 3 seconds.

Implicit methods used to avoid the unstable conditions during the iterations.

`((T(x-h) -2T(x)+T(x+h))/h^2 + (T(y-h) -2T(y)+T(y+h))/h^2)^n= (1/alpha)(T^n-T)/(△t)`

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.

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.

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.

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

Matlab codes for steady state

************************************************************************

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

******************************************************************

Matlab codes for transient explicit

******************************************************************

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

******************************************************************

Matlab codes for transient implicit

******************************************************************

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

*******************************************************************