Week 5 - Mid term project - Solving the steady and unsteady 2D heat conduction problem
AIM:-
To solve the 2D heat conduction equation
- Case 1: In steady state by the following iterative methods
- Gauss - Jacobi
- Gauss - Seidel
- Successive over relaxation (SOR)
- Case 2: In transient state by
- Explicit scheme
- Implicit scheme by
- Gauss - Jacobi
- Gauss - Siedel
- Successive over relaxation (SOR)
ASSUMPTIONS:-
- The equations will be solved for a Cartesian unit square domain
- The number of grid points in x and y directions will be equal
- The boundary conditions are same for all the above cases.
- The material properties are isotropic in nature
DOMAIN:-
BOUNDARY CONDITIONS:-
All the boundaries of the square domain are provided with Dirichlet’s B.C as shown in the above figure
CASE 1:
Steady state 2D heat conduction
The steady state 2D heat conduction equation is defined by the following linear, 2nd order partial differential equation.
`(∂^2 T)/(∂x^2 )+ (∂^2 T)/(∂y^2 )=0`
To solve this equation numerically we have to discretize using central difference scheme, which results in:-
`(T_(i-1,j )-2T_(i,j)+T_(i+1,j))/(Δx^2 ) + (T_(i,j-1 )-2T_(i,j)+T_(i,j+1))/(Δy^2 )=0`
Further simplification of this by taking on one side leads to
`T_(i,j) =1/K [ (T_(i-1,j )+ T_(i+1,j) )/(Δx^2 )+(T_(i,j-1 )+T_(i,j+1))/(Δy^2 )]` Eq. (1)
We have to solve the above equation at all the node points except the boundaries of the domain. The use of iterative solvers is done to accomplish this task. A MATLAB code is used to solve this equation by three methods Gauss Jacobi, Gauss- Siedel and SOR method.
GAUSS JACOBI METHOD
For Gauss Jacobi method Eq. (1) can be rewritten as
`T_p^(k+1)=1/k [(T_L+T_R)/(Δx^2 )+ (T_T+T_B)/(Δy^2 )]^k`
Where the subscripts , p = current node
K = current iteration
L = node left to the current node
R = node right to the current node
B = node above the current node
T = node below the current node
In Gauss Jacobi method , the RHS term contain values of temp at the current iteration.
GAUSS SEIDEL METHOD
`T_p^(k+1)=1/k [ (T_L^(K+1)+T_R)/(Δx^2 )+ (T_T+T_B^(K+1))/(Δy^2 ) ]`
In this method, the terms which are already calculated in the current iteration are used for the calculation of the next iteration.
SUCCESSIVE OVER RELAXATION METHOD
`T_p^(k+1)=α*(1/k [(T_L^(K+1)+T_R)/(Δx^2 )+ (T_T+T_B^(K+1))/(Δy^2 )])+ (1-α)T_p^k`
In this method, we multiply the Temperature values got from either Gauss Jacobi or Gauss Seidel method by a factor and then the value of temperature for previous iteration by . Therefore if
`α=1` SOR = Jacobi/Gauss Seidel
`α> 1` Convergence rate increases
`α< 1` Convergence rate decreases
MATLAB CODE:-
The code consists of a single scripts consisting of if and else loop which gives the results of GAUSS JACOBI if value of variable , 'method' = 1
GAUSS SEIDEL if value of variable , 'method' = 2
SOR if value of variable , 'method' = 3
%2D steady state heat conduction
close all
clear all
clc
%%INPUT PARAMETERS
L = 1; %Lenght of domain
B = 1; %Breadth of domain
nx = 20; %Grid points in x direction
ny = 20; %Grid points in y direction
error = 1; % Starting value of error
limit = 1e-4; %Allowable error between iterations
iteration = 0; %Starting value of iteration
method = 2; %METHOD = 1,2 and 3 for Jacobi,Siedel and SOR techniques respectively
alpha =1.7; %value for SOR relaxation factor
%%DISRETIZATION
x =linspace(0,L,nx); %Discretizing x domain
y= linspace(0,B,ny); %Discretizing y domain
dx = L/(nx-1); %grid spacing in x axis
dy = B/(ny-1); %grid spacing in y axis
k = 2*((1/dx^2) + (1/dy^2)); %value of constants in Eq.(1)
%%INITIALIZATON
T_old = 298*ones(nx,ny); %Initializing with ambient temperature
%%BOUNDARY CONDITIONS
T_old(2:nx-1,1) =400; %LEFT FACE = 400°C
T_old(2:nx-1,ny) = 800; %RIGHT FACE = 800°C
T_old(nx,2:ny-1) =900; %BOTTOM FACE = 900°C
T_old(1,2:ny-1)=600; %TOP FACE = 600°C
T_old(1,1) = (T_old(2,1) + T_old(1,2))/2; %Top left corner average value
T_old(1,ny) = (T_old(1,ny-1) + T_old(2,ny))/2; %Top right corner average value
T_old(nx,1) = (T_old(nx-1,1) + T_old(nx,2))/2; %Bottom left average value
T_old(nx,ny) = (T_old(nx,ny-1) + T_old(nx-1,ny))/2; %Bottom right average value
contourf(x,y,T_old);
%% GAUSS JACOBI METHOD
if method == 1
T_new = T_old; % initializing temp values
while error > limit %loop will run till error between iteration is below 1e-4
for i = 2:nx-1 %counts for x dimension leaving the boundary nodes
for j = 2:ny-1 %counts for y dimension leaving the boundary nodes
term1 = (1/(k*dx^2))*(T_old(i-1,j) + T_old(i+1,j));
term2 = (1/(k*dy^2))*(T_old(i,j-1) + T_old(i,j+1));
T_new(i,j) = term1+term2; % Eq. (1)
end
end
error = max(max(abs(T_old - T_new))); % finds max value of error from two matricies
iteration = iteration +1;
T_old =T_new; % prepares matrix for next iteration
end
a = contourf(x,y,T_new);
clabel(a);
set(gca,'YDIR','reverse');
title(['Gauss- Jacobi method , total iterations =',num2str(iteration)]);
colorbar;
saveas(gcf,'Gauss Jacobi method.png');
elseif method == 2
%% GAUSS SEIDEL METHOD
T_new = T_old;
while error > limit
for i = 2:nx-1
for j = 2:ny-1
term1 = (1/(k*dx^2))*(T_new(i-1,j) + T_old(i+1,j));
term2 = (1/(k*dy^2))*(T_new(i,j-1) + T_old(i,j+1));
T_new(i,j) = term1+term2;
end
end
error = max(max(abs(T_old - T_new)));
iteration = iteration +1;
T_old =T_new;
end
a = contourf(x,y,T_new);
clabel(a);
set(gca,'YDIR','reverse');
title(['Gauss- Siedel method , total iterations =',num2str(iteration)]);
colorbar;
saveas(gcf,'Gauss Seidel method.png');
elseif method == 3
%% SOR METHOD
T_new = T_old;
while error > limit
for i = 2:nx-1
for j = 2:ny-1
term1 = (1/(k*dx^2))*(T_new(i-1,j) + T_old(i+1,j));
term2 = (1/(k*dy^2))*(T_new(i,j-1) + T_old(i,j+1));
T_new(i,j) = alpha *(term1+term2) + (1- alpha)*T_old(i,j);
end
end
error = max(max(abs(T_old - T_new)));
iteration = iteration +1;
T_old =T_new;
end
a = contourf(x,y,T_new);
clabel(a);
set(gca,'YDIR','reverse');
title(['SOR Method, total iterations =',num2str(iteration)]);
colorbar;
saveas(gcf,'SOR method.png');
endRESULTS :-
From the above results it can be concluded that the SOR method gives the fastest result in terms of iterations required to converge. Whereas the Jacobi method requires largest number of iterations to converge. The number of iterations required for SOR to converge depend upon the value of its relaxation factor. From trial and error it was found to be 1.7.
CASE 2: Un-Steady 2D heat conduction
The unsteady state 2D heat conduction equation is defined by the following linear, 2nd order partial differential equation.
`(∂^2 T)/(∂x^2 )+ (∂^2 T)/(∂y^2 )=1/α ( ∂T)/(∂t)` Eq (2)
Here the constant term `α` is called thermal diffusivity, it determines how fast (in terms of time) the heat will flow from hot to cold surface.
EXPLICIT METHOD
We can write the above equation in numerical form by using central difference scheme for spatial terms and forward difference scheme for the time derivative.
`(T_(i,j)^(n+1)-T_(i,j)^n)/(Δt)=α((T_(i+1,j)^n-2T_(i,j)^n+T_(i-1,j)^n)/(Δx^2 )+(T_(i,j+1)^n-2T_(i,j)^n+T_(i,j-1)^n)/(Δy^2 ) )`
Where the subscripts i,j represent the x and y co-ordinate of the discretized domain. The superscript n defines the time-step. Since the RHS of the above equation only contains values of the current time step, which are known to us, this method is called explicit method of time integration.
We can simplify the above equation for coding as follows,
`T_p^(n+1)=T_p^n+αΔt((T_L-2T_p+T_R)/(Δx^2 )+(T_B-2T_p+T_L)/(Δy^2 ) )^n`
Since this is an explicit method there is a condition of stability associated with this which is as follows,
Therefore the terms here are so chosen that they sum up to be lower than 0.5.
`α` =0.001040
dt = 0.5
dx = dy = 0.0526
MATLAB CODE :-
%%Main program for WEEK 5 challenge
%2D transient state heat conduction explicit form
close all
clear all
clc
%%INPUT PARAMETERS
L = 1; %Length
B = 1; %Breadth
nx = 20; %nodes in x
ny = 20; %nodes in y
limit = 1e-4; %limit for error
alpha =0.001040; %Value of thermal diffusivity
T = 500; %Total time in secs
dt = 0.5; %Timestep
nt = T/dt; %no. of timesteps
%%DISRETIZATION
x =linspace(0,L,nx); %discretizing length
y= linspace(0,B,ny); %discretizing breadth
dx = x(2)-x(1); %x grid spacing
dy = dx; %y grid spacing
t = linspace(0,T,nt); %discretizing time
%constant k
k = (alpha*dt)/(dx^2);
CFL = 2*k;
%%INITIALIZATON
T_old = 298*ones(nx,ny);
%%BOUNDARY CONDITIONS
T_old(2:nx-1,1) =400; %LEFT FACE 400°C
T_old(2:nx-1,ny) = 800; %RIGHT FACE 800°C
T_old(nx,2:ny-1) =900; %BOTTOM FACE 900°C
T_old(1,2:ny-1)=600; %TOP FACE 600°C
T_old(1,1) = (T_old(2,1) + T_old(1,2))/2; %Top left corner
T_old(1,ny) = (T_old(1,ny-1) + T_old(2,ny))/2; %Top right corner
T_old(nx,1) = (T_old(nx-1,1) + T_old(nx,2))/2; %Bottom left
T_old(nx,ny) = (T_old(nx,ny-1) + T_old(nx-1,ny))/2; %Bottom right
%Boundary CONDITIONS
T_new = T_old;
% displaying initial conditions
for m = 1: nt %TIME LOOP
%SPACE LOOP
for j = 2:nx-1
for i = 2:ny-1
term1 = T_old(i,j);
term2 = T_old(i-1,j) - 2*T_old(i,j) + T_old(i+1,j);
term3 = T_old(i,j+1) - 2*T_old(i,j) + T_old(i,j-1);
T_new(i,j) = term1 + k*term2 +k*term3;
end
end
T_old = T_new;
end
a = figure(2);
contourf(x,y,T_new,'showtext','on');
axis square;
title(['Unsteady 2D heat conduction - Explicit method ','Time =', num2str(t(m))]);
xlabel('x axis');
ylabel('y axis');
set(gca,'YDIR','reverse');
colorbar;
saveas(a,'Explicit Method for Transient Heat conduction.png')
RESULTS
The results for unsteady explicit 2D heat conduction are run till 500 secs so that the contour plots near about match the steady state results. In unsteady state simulation, the temperature changes at different locations as time advances can be visualised .Moreover the steady state results can also be observed if the simulation is run for a longer amount of time.
IMPLICIT METHOD
If in Eq (2) the RHS term has values of temperature computed at n+1 timestep, then the method of time integration is called Implicit method.
`T_p^(n+1) ( 1+2k_(1 )+2k_2 )= T_p^n+k_1 ( T_L^(n+1)+T_R^(n+1) )+k_2 (T_T^(n+1)+T_B^(n+1) )`
Where , `k_1=(α Δt)/(Δx^2 )` & `k_2=(α Δt)/(Δy^2 )`
For solving the above equations we can use the above mentioned iterative methods whose formulas are as follows
Jacobi method
`(T_p^(n+1) ( 1+2k_(1 )+2k_2 ))^(k+1)= (T_p^n+k_1 ( T_L^(n+1)+T_R^(n+1) )+k_2 (T_T^(n+1)+T_B^(n+1) ) )^k`
Gauss Seidel Method
`(T_p^(n+1) ( 1+2k_(1 )+2k_2 ))^(k+1) = T_p^n+k_1 ( (T_L^(n+1) )^(k+1)+T_R^(n+1) )+k_2 (T_T^(n+1)+(T_B^(n+1) )^(k+1) )`
SOR Method
`(T_p^(n+1) )^(k+1) =α(( T_p^n+k_1 ( (T_L^(n+1) )^(k+1)+T_R^(n+1) )+k_2 (T_T^(n+1)+(T_B^(n+1) )^(k+1) ))/(1+2k_1+2k_2 ))+(α-1)*(T_p^(n+1) )^k`
Where subscripts and superscripts have the same meaning as discussed earlier
MATLAB CODE
The code consists of a single scripts consisting of if and else loop which gives the results of
GAUSS JACOBI if value of variable , 'method' = 1
GAUSS SEIDEL if value of variable , 'method' = 2
SOR if value of variable , 'method' = 3
%%Main program for WEEK 5 challenge
%2D transient state heat conduction
close all
clear all
clc
%%INPUT PARAMETERS
L = 1;
B = 1;
nx =20;
ny = 20;
limit = 1e-4;
alpha =0.001040;
T = 500;
dt = 0.5;
nt = (T/dt)+1 ;
omega = 1.1;
method = 3;
iteration = 0;
%%DISRETIZATION
x =linspace(0,L,nx);
y= linspace(0,B,ny);
dx = L/(nx-1);
dy = B/(ny-1);
t = linspace(0,T,nt);
%k1 and k2 values
k1 = (alpha*dt)/(dx^2);
k2 = k1;
k3 = (1 + 2*k1 + 2*k2);
%%INITIALIZATON
T_old = 300*ones(nx,ny);
%%BOUNDARY CONDITIONS
T_old(2:nx-1,1) =400; %LEFT FACE
T_old(2:nx-1,ny) = 800; %RIGHT FACE
T_old(nx,2:ny-1) =900; %BOTTOM FACE
T_old(1,2:ny-1)=600; %TOP FACE
T_old(1,1) = (T_old(2,1) + T_old(1,2))/2; %Top left corner
T_old(1,ny) = (T_old(1,ny-1) + T_old(2,ny))/2; %Top right corner
T_old(nx,1) = (T_old(nx-1,1) + T_old(nx,2))/2; %Bottom left
T_old(nx,ny) = (T_old(nx,ny-1) + T_old(nx-1,ny))/2; %Bottom right
T_new = T_old;
T_prev = T_old;
%% GAUSS JACOBI METHOD
if method == 1
for m = 1: nt
error = 1; %% TO RESET THE ERROR VALUE SO THAT WHILE LOOP GETS EXECUTED FOR THE CURRENT TIMESTEP!!!
while error > limit
for i = 2:nx-1
for j = 2:ny-1
left = T_old(i-1,j);
right = T_old(i+1,j);
top = T_old(i,j+1);
bottom = T_old(i,j-1);
term1 = T_prev(i,j)/k3;
term2 = (k1/k3)*(left + right);
term3 = (k2/k3)*(top + bottom);
T_new(i,j) = term1 +term2 +term3;
end
end
error = max(max(abs(T_old - T_new)));
iteration = iteration +1;
T_old = T_new;
end
T_prev = T_new;
end
a = figure(2);
contourf(x,y,T_new,'showtext','on');
axis square;
title(['Transient Implicit - Jacobi method ','Iteration =', num2str(iteration)]);
xlabel('x axis');
ylabel('y axis');
set(gca,'YDIR','reverse');
colorbar;
saveas(a,'Jacobi Implicit Method for Transient Heat conduction.png')
%% GAUSS SEIDEL METHOD
elseif method == 2
for m = 1: nt
error = 1; %% TO RESET THE ERROR VALUE SO THAT WHILE LOOP GETS EXECUTED FOR THE CURRENT TIMESTEP!!!
while error > limit
for i = 2:nx-1
for j = 2:ny-1
left = T_new(i-1,j);
right = T_old(i+1,j);
top = T_old(i,j+1);
bottom = T_new(i,j-1);
term1 = T_prev(i,j)/k3;
term2 = (k1/k3)*(left + right);
term3 = (k2/k3)*(top + bottom);
T_new(i,j) = term1 +term2 +term3;
end
end
error = max(max(abs(T_old - T_new)));
iteration = iteration +1;
T_old = T_new;
end
T_prev = T_new;
end
a = figure(2);
contourf(x,y,T_new,'showtext','on');
axis square;
title(['Transient Implicit - Seidel method ','Iteration =', num2str(iteration)]);
xlabel('x axis');
ylabel('y axis');
set(gca,'YDIR','reverse');
colorbar;
saveas(a,'Seidel Implicit Method for Transient Heat conduction.png')
%% SOR METHOD
else
for m = 1: nt
error = 1; %% TO RESET THE ERROR VALUE SO THAT WHILE LOOP GETS EXECUTED FOR THE CURRENT TIMESTEP!!!
while error > limit
for i = 2:nx-1
for j = 2:ny-1
left = T_new(i-1,j);
right = T_old(i+1,j);
top = T_old(i,j+1);
bottom = T_new(i,j-1);
term1 = T_prev(i,j)/k3;
term2 = (k1/k3)*(left + right);
term3 = (k2/k3)*(top + bottom);
T_new(i,j) = omega*(term1 +term2 +term3) + (1-omega)*T_new(i,j);
end
end
error = max(max(abs(T_old - T_new)));
iteration = iteration +1;
T_old = T_new;
end
T_prev = T_new;
end
a = figure(2);
contourf(x,y,T_new,'showtext','on');
axis square;
title(['Transient Implicit - SOR method ','Iteration =', num2str(iteration)]);
xlabel('x axis');
ylabel('y axis');
set(gca,'YDIR','reverse');
colorbar;
saveas(a,'SOR Implicit Method for Transient Heat conduction.png')
end
RESULTS
The results show similar results to explicit method and steady state results. The total number of iterations required for SOR method is lowest when the value of relaxation factor is 1.1. The Jacobi method takes lots of iteration compared to other methods for giving the same results.
CONLCUSION
The 2D heat conduction equation was solved for both steady and unsteady state. The steady state method and implicit transient method is solved with Jacobi,Seidel and SOR methods . The number of iterations for SOR method is very less compared to other methods and the rate of convergence depends on the under relaxation factor. The SOR method gives fastest convergence only at some particular value . It is also observed that wrong values can also lead to the slowest convergence.
The Explicit method on the other hand is to be used only when CFL condition is below 0.5 otherwise the solution becomes unstable. This is not seen with implicit methods which are fairly stable even for values exceeding the CFL condition.