Week 5 - Mid term project - Solving the steady and unsteady 2D heat conduction problem
The 2D heat conduction equation in steady state form is given as:
`(del^2T)/(delx^2)+(del^2T)/(dely^2)= 0`
from above equation we can say that the temperature is the function of both x co-ordinate and y co-ordinate. To solve the equation we have will use all the three iterative techniques which are:
1) Jacobi.
2) Gauss Siedal.
3) SOR (Successive over-relaxation)
In the jacobi method we are calculating the next space derivating term from the boundary condition value and initial value. For each iteration we are just using the value which is given at boundary.
In gauss siedal method we are calculating the next space derivative term from the boundary condition and the previous value which is calculated in further space.
In Successive over-relaxation we can use both jacobi method and gauss seidal method to calculate the next space derivating in which we are allowing one constant whose value can be changed to get the convergence more quickly. The value is ranged from 0 to 2. If the value is less than 1 then it is called under-relaxed method and if the value is between 1 to 2 then it is called over-relaxed method.
There is an optimum value of that constant (w) for which we will have at least number of iteration steps needed to reduce the residual by the factor of 10. It will provide the higher convergence rate but it is very rarely used and in some occasions. The optimum value is also depend on the spectral radius, if spectral radius is smaller the rate of convergence is faster.
To solve the equation numerically we will just take some notation
current point = `T_P`
left point in space = `T_L`
right point in space = `T_R`
top point in space = `T_T`
bottom point in space = `T_B`
Therefore, for the steady state function the final equation is given as:
`T_P = (1/k)*((T_L+T_R)/(Deltax^2))+(1/k)*((T_T+T_B)/(Deltay^2))`
Can be further written as:
`T_P = (1/k)*((T_L+T_R)/(Deltax^2)+(T_T+T_B)/(Deltay^2))`
Where `k = 2*((Deltax^2+Deltay^2)/(Deltax^2*Deltay^2))`
Assuming grid point notation for x-direction as 'i' and y-direction for 'j' the above equation can be written in the programming format for jacobi iterative method as:
`T(i,j) = (1/k){((Told(i-1,j)+Told(i+1,j))/dx^2)+((Told(i,j-1)+Told(i,j+1))/dy^2)}`
The assumption we have made before the programing the code:
1) The unit grid is consider which is nx = ny.
2) Boundary conditions are:
a) Temperature at left boundary = 400k
b) Temperature at right boundary = 800k
c) Temperature at Top boundary = 600k
d) Temperature at bottom boundary = 900k.
3) Tolerance = 1e-5.
4) Error = 9e9.
5) Thermal diffusivity = 1.4
6) Lenght of Domain = 1m
7) Stability Criteria = `alpha*(dt/(dx^2))`
***********************************************************************************
The program code for steady state function can be given as:
%steady state equation by jacobi, guass siedal and SOR
close all
clear all
clc
%inputs for the simulation
nx = 10;
ny = nx;
x = linspace(0,1,nx);
y = linspace(0,1,ny);
dx = x(2)-x(1);
dy = y(2)-y(1);
[xx,yy] = meshgrid(x,y);
%Boundary condition defination
T = ones(nx,ny);
T(:,1) = 400; %Left boundary Temperature value
T(:,end) = 800; %Right boundary Temperature value
T(1,:) = 900; %Bottom boundary Temperature value
T(end,:) = 600; %Top boundary Temperature value
Told = T; %Value assinging to old temperature format.
%we will creat a variable to call the specific iterative solvers
iterative_solver = str2double(input('Please enter the No. between 1 to 3 for selection of solver: ', 's'));
while isnan(iterative_solver) || fix(iterative_solver) ~= iterative_solver
iterative_solver = str2double(input('Please enter and INTEGER only to select solver: ', 's'));
end
error_SOR = 9e9;
error_GS = 9e9;
error_jacobi = 9e9;
tol = 1e-4; %if error is less than 1e-4 we can say the solution converge
k = 2*((dx^2+dy^2)/(dx^2*dy^2));
%Jacobi Iteration
tic
if iterative_solver ==1;
jacobi_iteration = 1;
while (error_jacobi>tol)
for i = 2:nx-1;
for j = 2:ny-1;
T(i,j) = (((Told(i-1,j)+Told(i+1,j))/dx^2)+((Told(i,j-1)+Told(i,j+1))/dy^2))/k;
end
end
error_jacobi = max(max(abs(Told-T)));
Told = T;
jacobi_iteration = jacobi_iteration + 1;
figure(1)
[c,d]= contourf(xx,yy,T);
colormap(jet)
colorbar
title_text = sprintf('jacobi Iteration No. = %d @Steady state',jacobi_iteration);
title(title_text)
xlabel('X-axis')
ylabel('Y-axis')
clabel(c,d)
pause(0.03)
end
fprintf('Number of Iteration in jacobi = %d',jacobi_iteration);
end
time_elapsed_Jacobi = toc;
%Gauss Siedal Iteration
tic
if iterative_solver ==2;
GS_iteration = 1;
while (error_GS>tol)
for i = 2:nx-1;
for j = 2:ny-1;
T(i,j) = (((T(i-1,j)+Told(i+1,j))/dx^2)+((T(i,j-1)+Told(i,j+1))/dy^2))/k;
end
end
error_GS = max(max(abs(Told-T)));
Told = T;
GS_iteration = GS_iteration + 1;
figure(1)
[c,d]= contourf(xx,yy,T);
colormap(jet)
colorbar
title_text = sprintf('jacobi Iteration No. = %d @Steady state',GS_iteration);
title(title_text)
xlabel('X-axis')
ylabel('Y-axis')
clabel(c,d)
pause(0.03)
end
fprintf('Number of Iteration in jacobi = %d',GS_iteration);
end
time_elapsed_GS = toc;
%SOR Iteration
tic
if iterative_solver ==3;
w = input('Enter Value of Relaxation factor between 1.5 and 1.8:');
SOR_iteration = 1;
while (error_SOR>tol)
for i = 2:nx-1;
for j = 2:ny-1;
T_GS = (((T(i-1,j)+Told(i+1,j))/dx^2)+((T(i,j-1)+Told(i,j+1))/dy^2))/k;
T(i,j) = w*(T_GS)+(1-w)*(Told(i,j));
end
end
error_SOR = max(max(abs(Told-T)));
Told = T;
SOR_iteration = SOR_iteration + 1;
figure(1)
[c,d]= contourf(xx,yy,T);
colormap(jet)
colorbar
title_text = sprintf('jacobi Iteration No. = %d @Steady state',SOR_iteration);
title(title_text)
xlabel('X-axis')
ylabel('Y-axis')
clabel(c,d)
pause(0.03)
end
fprintf('Number of Iteration in jacobi = %d',SOR_iteration);
end
time_elapsed_SOR = toc;
Explanation of Code structure:
1) The code start with the assumption we have made during the initial process of problem defination.
2) We have defined the boundary condition as per the given problem statements.
3) There is use of while loop for the entering the number of specific solver. If user enter any decimal value it will again show a input string which ask user to enter the 'INTEGER' number only to call the solver.
4) If user enter no. 3 as the solver he will directed to the SOR solver in which program will ask user to enter the value for w as omaga between 1.5 to 1.8 in which he will get less number of iteration.
5) The plot is created in the contour format to show the temperature filed varing in the graph.
6) For each iterative solver we are using tic toc function to calculate elapsed time for each program.
Results for Steady State Equation:
1) Jacobi Solver:
Time required to complete the program = 81.8664 sec.
2) Gauss Seidal:
Time required to complete the program = 43.7115 sec.
3) SOR Iterative method:
Considering Omega value as 1.5
Time required to complete the program = 19.3326 sec.
Conclusion:
From the above three contour graph we can see the difference in the number of iteration required for calculation of Temperature. The SOR method gives least amount of iteration as compare to the jacobi and gauss seidal approach. The program is run for different value of omega in the SOR solver and it found that the value of iteration are considerably lower in range of 1.5 to 1.8. Also from the time elapsed for each solver as we can see SOR provide faster results as compare to others.
***********************************************************************************
Unsteady State Equation
The 2D unsteady state equation is given by:
`(delT)/(delt) = alpha*((del^2T)/(delx^2)+(del^2T)/(dely^2))`
The above equation can be solved using explicitly and implicitly. In the explicit solution it does not required any iterative solution as the next time step and spacial derivative is solved on the basis of boundary conditions and initial condition. In the implicit approach as in the right hans side of equation we are getting some unknown terms which need to be evaluated hence it can be solved using the iterative solvers.
The assumption that we have to made for this kind of simulation are the time step size as the temperature is the function of time in left hand side of equation.
Solution of unsteady state equation with explicit scheme.
`((T_P^(n+1) - T_P^n)/(Deltat)) = alpha*((T_L-2T_P+T_R)/(Deltax^2)+(T_T-2T_P+T_B)/(Deltay^2))^n`
All the value in above equation are written at 'n' time level and left hand side only one 'n+1' time level.
Assuming grid point notation for x-direction as 'i' and y-direction for 'j' the above equation can be written in the programming format for jacobi iterative method as:
`T^(n+1)(i,j) = T^n+k1*(T^nold(i-1,j)-2T^n(i,j)+T^n(i+1,j))+k2*(T^n(i,j-1)-2T^n(i,j)+T^n(i,j+1))`
where `k1 = (alphaDeltat)/(Deltax^2)` & `k2 = (alphaDeltat)/(Deltay^2)`
The above constant are the CFL constant for the explicit scheme. So selection of this constant are very important as stability of solution is depend on them. For stabilize solution the value of CFL number should be less than and equal to 0.5.
The assumption will be prety much same for the explicit scheme only the difference is to provide number of time step size. Also we will increase the grid size to get the convergence and stability to the solution.
The code for Explicit scheme for Unsteady state equation:
%Unsteady state equaiton with Implicit scheme
clear all
close all
clc
%inputs for the simulation
nx = 51;
ny = nx;
x = linspace(0,1,nx);
y = linspace(0,1,ny);
dx = x(2)-x(1);
dy = y(2)-y(1);
dt = 0.01;
[xx,yy] = meshgrid(x,y);
alpha = 1.88e-5;
%Boundary condition defination
T = ones(nx,ny);
T(:,1) = 400; %Left boundary Temperature value
T(:,end) = 800; %Right boundary Temperature value
T(1,:) = 900; %Bottom boundary Temperature value
T(end,:) = 600; %Top boundary Temperature value
Told = T; %Value assinging to old temperature format.
%Defining corner value
T(1,1) = (400+900)/2;
T(1,ny) = (900+800)/2;
T(nx,1) = (400+600)/2;
T(nx,ny) = (600+800)/2;
error_jacobi = 9e9;
tol = 1e-4; %if error is less than 1e-4 we can say the solution converge
k1 = alpha*(dt/dx^2);
k2 = alpha*(dt/dy^2);
%we will creat a variable to call the specific iterative solvers
iterative_solver = 1;
tic;
if iterative_solver == 1
T_jacobi_explicit = Jacobi_iter_explicit(xx,yy,nx,ny,T,k1,k2,Told);
end
time_elapsed_jacobi = toc
The Function defined as:
function T = Jacobi_iter_explicit(xx,yy,nx,ny,T,k1,k2,Told)
jacobi_iteration = 1;
for k=1:140000
for i = 2:nx-1
for j = 2:ny-1
Term1 = Told(i,j);
Term2 = k1*(Told(i-1,j)-2*Told(i,j)+Told(i+1,j));
Term3 = k2*(Told(i,j+1)-2*Told(i,j)+Told(i,j-1));
T(i,j) = Term1+Term2+Term3;
end
end
error_jacobi = max(max(abs(Told-T)));
Told = T;
jacobi_iteration = jacobi_iteration + 1;
end
figure(1)
[c,d]= contourf(xx,yy,T);
colormap(jet)
colorbar
title_text = sprintf('jacobi Iteration No. = %d @UnSteady state',jacobi_iteration);
title(title_text)
xlabel('X-axis')
ylabel('Y-axis')
clabel(c,d)
pause(0.03)
fprintf('Number of Iteration in jacobi = %d',jacobi_iteration);
endExplanation of Code:
1) In this code we have used the function called Jacobi_iter_explicit to provide at main function for calculation of temperature.
2) In this code we not used the while loop as it mostly not required at this stage for the convergence.
Results:
As we can see in the above image how the temperature profile is changes with respect to the time.
***********************************************************************************
Solution of unsteady state equation with Implicit scheme.
The equation for unsteady state equation with implicit scheme can be given as:
`((T_P^(n+1) - T_P^n)/(Deltat)) = alpha*((T_L-2T_P+T_R)/(Deltax^2)+(T_T-2T_P+T_B)/(Deltay^2))^(n+1)`
Assuming grid point notation for x-direction as 'i' and y-direction for 'j' the above equation can be written in the programming format for jacobi iterative method as:
`T_P^(n+1) = (1/(1+2*k1+2*k2))*((T^n(i,j)+k1*(T^(n+1)(i-1,j)+T^(n+1)(i+1,j))+k2*((T^(n+1)(i,j-1)+T^(n+1)(i,j+1)))`
As the above equation contains (n+1) on the both side hence we have to use the iterative solver to get temperature value. The equation look bit tidy so the effort is taken to minimize the error while writing it by introducing some terms.
In the Jacobi iteration scheme as we have already discussed that it will calculate the each time step value and spacial value from the old temperature which is given as per the boundary condition. Also for the next time step it wont use the previous time step value as a input.
The most of input value are same the one improvement value that we have to use is T_prev i.e. temperature at previous time step because we have 'T^n' on the right hand side.
The code for Unsteady state equation with Implicit scheme is given as:
%Unsteady state equation solution by iterative method
close all
clear all
clc
%inputs for the simulation
nx = 41;
ny = nx;
L = 1;
x = linspace(0,1,nx);
y = linspace(0,1,ny);
dx = x(2)-x(1);
dy = y(2)-y(1);
dt = 0.01;
[xx,yy] = meshgrid(x,y);
alpha = 1.88e-5;
%Boundary condition defination
T = 298*ones(nx,ny);% Initial Value
T(:,1) = 400; %Left boundary Temperature value
T(:,end) = 800; %Right boundary Temperature value
T(1,:) = 900; %Bottom boundary Temperature value
T(end,:) = 600; %Top boundary Temperature value
%Calculation of average value at corner.
T(1,1) = (400+900)/2;
T(1,ny) = (900+800)/2;
T(nx,1) = (400+600)/2;
T(nx,ny) = (600+800)/2;
Told = T;
T_pre = Told;
%error_jacobi = 9e9;
error_GS = 9e9;
error_SOR = 9e9;
tol = 1e-4; %if error is less than 1e-4 we can say the solution converge
k1 = alpha*(dt/dx^2);
k2 = alpha*(dt/dy^2);
CFL = alpha*(dt/(dx^2));
%we will creat a variable to call the specific iterative solvers
iterative_solver = str2double(input('Please enter the No. between 1 to 3 for selection of solver: ', 's'));
while isnan(iterative_solver) || fix(iterative_solver) ~= iterative_solver
iterative_solver = str2double(input('Please enter and INTEGER only to select solver: ', 's'));
end
%Jacobi Iteration
tic;
if iterative_solver == 1
jacobi_iterative = 1;
for a = 1:120000
error_jacobi = 9e9;
while (error_jacobi>tol)
for i = 2:nx-1
for j = 2:ny-1
Term1 = (1+(2*k1)+(2*k2))^-1;
Term2 = k1*Term1;
Term3 = k2*Term1;
P = (Told(i-1,j)+Told(i+1,j));
Q = (Told(i,j-1)+Told(i,j+1));
T(i,j) = (T_pre(i,j)*Term1)+(Term2*P)+(Term3*Q);
end
end
error_jacobi = max(max(abs(Told-T)));
Told = T;
jacobi_iterative = jacobi_iterative + 1;
end
T_pre = T;
end
time_elapsed_jacobi = toc;
contourf(xx,yy,T,'ShowText','on')
colormap(jet)
colorbar
xlabel('X-axis')
ylabel('Y-axis')
colorbar
title_bar = sprintf('2D Implicit Unsteady State Jacobi iteration required = %d and Time taken %f sec:',jacobi_iterative,time_elapsed_jacobi);
title(title_bar)
fprintf('No. of Iteration required for jacobi = %d and Time taken in sec = %f',jacobi_iterative,time_elapsed_jacobi)
end
%Gauss Seidal
tic;
if iterative_solver == 2
GS_iterative = 1;
for a = 1:120000
error_GS = 9e9;
while (error_GS>tol)
for i = 2:nx-1
for j = 2:ny-1
Term1 = (1+(2*k1)+(2*k2))^-1;
Term2 = k1*Term1;
Term3 = k2*Term1;
P = (T(i-1,j)+Told(i+1,j));
Q = (T(i,j-1)+Told(i,j+1));
T(i,j) = (T_pre(i,j)*Term1)+(Term2*P)+(Term3*Q);
end
end
error_GS = max(max(abs(Told-T)));
Told = T;
GS_iterative = GS_iterative + 1;
end
T_pre = T;
end
time_elapsed_GS = toc;
contourf(xx,yy,T,'ShowText','on')
colormap(jet)
colorbar
xlabel('X-axis')
ylabel('Y-axis')
colorbar
title_bar = sprintf('2D Implicit Unsteady State GS iteration required = %d and Time taken %f sec:',GS_iterative,time_elapsed_GS);
title(title_bar)
fprintf('No. of Iteration required for GS = %d and Time taken in sec = %f',GS_iterative,time_elapsed_GS);
end
%SOR
tic;
if iterative_solver == 3
SOR_iterative = 1;
w = input('Enter Value of Relaxation factor between 1.5 and 1.8:');
for a = 1:120000
error_SOR = 9e9;
while (error_SOR>tol)
for i = 2:nx-1
for j = 2:ny-1
Term1 = (1+(2*k1)+(2*k2))^-1;
Term2 = k1*Term1;
Term3 = k2*Term1;
P = (T(i-1,j)+Told(i+1,j));
Q = (T(i,j-1)+Told(i,j+1));
T(i,j) = w*((T_pre(i,j)*Term1)+(Term2*P)+(Term3*Q))+(1-w)*(T_pre(i,j));
end
end
error_SOR = max(max(abs(Told-T)));
Told = T;
SOR_iterative = SOR_iterative + 1;
end
T_pre = T;
end
time_elapsed_SOR = toc;
contourf(xx,yy,T,'ShowText','on')
colormap(jet)
colorbar
xlabel('X-axis')
ylabel('Y-axis')
colorbar
title_bar = sprintf('2D Implicit Unsteady State SOR iteration required = %d and Time taken %f sec:',SOR_iterative,time_elapsed_SOR);
title(title_bar)
fprintf('No. of Iteration required for SOR = %d and Time taken in sec = %f',SOR_iterative,time_elapsed_SOR);
end
The whole code comprises of each type of Iterative methods.
Explanation of Code:
1) There is no such differece in this code with respect to steady state equation code and unsteady state explicit code.
2) The only difference is the time value and the updation of T_pre value with respect to orginal temperaure is most important.
3) If you put your plot after the calculation of T_pre value then it will take hell lot time to compute all the value, hence it is better it should be done after time loop.
Results:
1) Jacobi Iteration
2) Gauss Seidal Iteration
3) SOR Iterative method
for omega = 1.5
Conclusion:
1) From the above three method as we can see the time taken for implicit scheme is lower if you take well approach of your coding.
2) The time required for Jacobi method is less as compare to the other two method which is very opposite in steady state type of equation.
3) But as for comparision the SOR method is more stable and reliable for the calculation of Unsteady state equation.