Week 5.1 - Mid term project - Solving the steady and unsteady 2D heat conduction problem
Week 5.1 - Mid term project - Solving the steady and unsteady 2D heat conduction problem
1. Steady state analysis & Transient State Analysis
Solve the 2D heat conduction equation by using the point iterative techniques that were taught in the class. The Boundary conditions for the problem are as follows;
Top Boundary = 600 K
Bottom Boundary = 900 K
Left Boundary = 400 K
Right Boundary = 800 K
You will implement the following methods for solving implicit equations.
- Jacobi
- Gauss-seidel
- Successive over-relaxation
Your absolute error criteria is 1e-4
Note: Make sure to use "INSERT CODE SAMPLE" option to paste the code here. Post screenshots of the plots and use the projects.skill-lync.com page to publish your findings.
2. You will write code for both implicit and explicit schemes (This is for the transient analysis only)
3. Compare the convergence rate [Number of iterations] of the steady state and transient [explicit and implicit] simulations and justify the numerics.
Hint:
The old values should be updated inside the convergence (error criteria) loop, while the previous values should be updated after each time loop.
Steady state analysis & Transient State Analysis
Solve the 2D heat conduction equation by using the point iterative techniques that were taught in the class. The Boundary conditions for the problem are as follows;
Top Boundary = 600 K
Bottom Boundary = 900 K
Left Boundary = 400 K
Right Boundary = 800 K
You will implement the following methods for solving implicit equations.
In order to solve for dependent variable iterative solvers are used. Iterative solvers used are
1. Jacobi iterative solver
2. Gauss seidal iterative solver
3. Successive over relaxation iterative solver.
2D Heat conduction Equation
The time derivative term
Discretization of Time derivative using forward difference scheme since the time derivative is varying in time direction which can be solved by marching method.
Discretization of spacial derivatives using central differencing scheme since the spacial derivatives are elliptical in nature, which uses information from neighbouring points to obtain a solution. Below are the discretization for x and y spacial derivatives.
From the PDE heat conduction equation,
where n is current time step and n+1 is next time step.
i,i+1 are the neighbouring points in x direction, j,j+1 are the neighbouring points in y direction.
Iterative Solver:
1. Jacobi iterative solver:
In this method the temperature at n+1th iteration is found using old values of the temperature defined at nth iterations. Using the approximation obtained iterations procedure it continued untill the desired accuracy of the solution is obtained.
2. Gauss Seidal iterative solver:
In this type iterative solver, temperature at n+1 iterations uses values of n+1 iterations that have already been computed. By this way old approximated values are overwritten with new approximate values for the iterations performed. Unlike jacobi, Gauss seidal require one storage vector which can be overwritten with each computation. This helps for the solution to coverge faster than the jacobi method.
3. Successive over relaxation iterative solver.
In this type of iterative solver relaxation factor is used which helps in increasing the convergence rate. It goes by the formula
where
when
This condition is used when the temperature gradient is high and due to this condition the time step does not over shoot and helps to converge at a faster rate.
when
This condition is used when the temperature gradient is low and due to this condition larger time steps are used to obtain the steady state i.e., for the values to converge at faster rate.
T(GS) or T(jacobi) is the temperature of value for n+1th iteration using gauss seidal or jacobi iterative solver.
Steady State 2D Heat Conduction Equation :
Steady state is where the fluid properties does not change with respect to time.
From the 2D heat conduction equation since the system is in steady state time derivative is zero and diffusivity is constant.
Hence we the equation for steady state as:-
Discretization for Steady State:
Solving for T(i,j), we get
, 
Solving the 2D Steady Heat Conduction equation implicitly.
Formula for Jacobi iterative solver:-
Formula for Gauss Seidal Iterative Solver:-
Formula for Successive over relaxation iterative solver:-
Transient state:
Unsteady state is where the fluid properties change with respect to time.
The formula for 2D transient heat conduction equation is :
Discretization for Unsteady State heat conduction equation:-
Solving 2D transient heat conduction equation explicitly we get,
Stability Criteria:
Solving for 2D Transient Heat Conduction equation implicitly we get,
here n+1 is the next time step and n is the preceeding time step.
Formula for jacobi iterative solver:-
Formula for Gauss seidal
Formula for Successive over relaxation
Boundary Conditions:
Boundary conditions are defined at the boundaries of independent variable of the function. In this problem mesh size of 10x10 is considered for which the boundary conditions are defined at the extreme rows and columns of the mesh.
For the mesh size 10x10 shown above the corner boxes are average of the two boundaries intersecting at that point. The values other than the boundaries are taken as guess values which is in this case in 300K (initial condition for transient heat conduction equation).
Matlab code for solving the 2D Steady state heat conduction equation by three iteration techniques.
clear all
close all
clc
% Steady state analysis
n = 50; % number of grid points
x = linspace(0,1,n);
dx = x(2)-x(1);
dy = dx;
y = linspace(0,1,n);
T = 298*ones(n,n); % initial temperature condition
T(:,1) = 400; % K ,temperature left boundary condition
T(end,:) = 600; % K ,temperature top boundary condition
T(:,end) = 800; % K ,temperature right boundary condition
T(1,:) = 900; % K ,temperature bottom boundary condition
tol = 1e-4; % absolute error criteria
Told1 = T; Told2 = T; Told3 = T;
T1 = T; T2 = T; T3 = T;
% method: 1.Jacobi Method, 2.Gauss Seidel Method, 3.Successive Over
% Relaxation Method (SOR)
for method = 1:3
% jacobi method
if method == 1
jac_error = 9e9; % initial value of error
jacobi_iter = 0;
tic
while jac_error>tol
jacobi_iter = jacobi_iter + 1;
for i = 2:(n-1)
for j = 2:(n-1)
T1(i,j) = (Told1(i-1,j) + Told1(i+1,j) + Told1(i,j-1) + Told1(i,j+1))/4;
end
end
jac_error = max(max(abs(T1-Told1)));
Told1 = T1;
end
analysis_time1 = toc
figure(1)
contourf(x,y,T1,'Showtext','on')
colormap(jet)
colorbar
xlabel('X axis')
ylabel('Y axis')
title(sprintf('2D Steady State Analysis (Jacobi Method) \nTemperature Distribution over a plate \nNumber of iterations = %d \nAnalysis time = %f seconds',jacobi_iter,analysis_time1))
pause(0.02)
% gauss seidel method
elseif method == 2
gauss_error = 9e9; % initial value of error
gauss_iter = 0;
tic
while gauss_error>tol
gauss_iter = gauss_iter + 1;
for i = 2:(n-1)
for j = 2:(n-1)
T2(i,j) = (T2(i-1,j) + Told2(i+1,j) + T2(i,j-1) + Told2(i,j+1))/4;
end
end
gauss_error = max(max(abs(T2-Told2)));
Told2 = T2;
end
analysis_time2 = toc
figure(2)
contourf(x,y,T2,'Showtext','on')
colormap(jet)
colorbar
xlabel('X axis')
ylabel('Y axis')
title(sprintf('2D Steady State Analysis (Gauss Seidel Method) \nTemperature Distribution over a plate \nNumber of iterations = %d \nAnalysis time = %f seconds',gauss_iter,analysis_time2))
pause(0.02)
% SOR method
w = 2/(1+sin(pi*dx)); % relaxation parameter
elseif method == 3
sor_error = 9e9; % initial value of error
sor_iter = 0;
tic
while sor_error>tol
sor_iter = sor_iter + 1;
for i = 2:(n-1)
for j = 2:(n-1)
T3(i,j) = (1-w)*T3(i,j) + w*(T3(i-1,j) + Told3(i+1,j) + T3(i,j-1) + Told3(i,j+1))/4;
end
end
sor_error = max(max(abs(Told3-T3)));
Told3 = T3;
end
analysis_time3 = toc
figure(3)
contourf(x,y,T3,'Showtext','on')
colormap(jet)
colorbar
xlabel('X axis')
ylabel('Y axis')
title(sprintf('2D Steady State Analysis (SOR Method) \nTemperature Distribution over a plate \nNumber of iterations = %d \nAnalysis time = %f seconds',sor_iter,analysis_time3))
pause(0.02)
end
end
Report
2D Steady state heat condution equation:
First we have to descritize the 2D steady state heat conduction equation, so that it can be solved numerically in matlab. After descritizing the equation we can solve the equation by three different iteration techniques. These are Jacobi method, Gauss Seidel method and Successive over relaxation method.
First we define the number of grid points n=50, number of grid points in x direction is equal to the number of grid points in y direction. The metal plate is a thin square plate of size 1 m. After this we define the grit size dx and dy.
Then we deifne the initial temperature condition to 298 kelvin accross entire nxn meshgrid, by using 'ones' command. And then we define the boundary temperature conditions according to the values given to us. Absoulte error criterio is equal to 1e-4, the error in the solution is compared after each iterations with absolute error criterio and if it is more than the absolute error, then the iterations keeps on going until the value of error is less than the absolute error criteria.
1.1 Jacobi iteration method
Formula used:
The above formula is used to calculate the value of T at time interval 'n+1'. In the above formula the RHS is assumed to have some intitial value, so the T at 'n+1' is assumed. Then for the first value fo i,j the value of T at 'n+1' is calculated. After first iteration the difference between the calculated value of T and previously assumed value T is calculated. If it more than the absolute error criteria or tolerance value, then the iteration keeps on going until the difference is less than the tolarance. Told value keep on uptading after each iteration inside the convergence loop.
1.2 Gauss seidel method
Formula used:
The above formula is similar to the jacobi method but when the first iteration is run to calculate the value fo T at a node, then this updated value of T is used to calculate the value at the next node. This technique decreases the convergence time.
1.3 Successive over relaxation method
Formula used:
This technique can be applied to both jacobi and gauss Seidel method, in this method a relaxation parameter 'w' is used which over estimates the value and leads to an increase value in T for each iteration, in this way the convergence is achived much faster than the above two methods.
Formula for relxation parameter:
We used 'tic' and 'toc' command to calculate the time of convergence for each of the three iteration method. After each convergence loop we plot the result on a contour plot between x,y,T. Contour plot visulises the temperature distribution over the plate and also gives us the number of iterations and time needed for convergence.
Graphs
Conclusion:
- The three graphs above shows the temperature distribution over a thin metal plate. And also displays the number of iterations and analysis time.
- From the above graphs we can say that jacobi method has maximum number of iterations than the rest of the methods. And convergence time is also more than other methds.
- Number of iterations and convergence time for Gauss Seidel method is less than jacobi method but more than SOR method.
- The successive over relaxation method has least number of iterations and it took least time for convergence.
- In 2D steady state analysis SOR method is fastest.
2. TRANSIENT STATE ANALYSIS
Matlab code for solving the 2D Transient state heat conduction equation by Explicit method and Implicit method.
clear all
close all
clc
% transient heat conduction
n = 50; % number of grid points
x = linspace(0,1,n);
dx = x(2) - x(1);
dy = dx;
y = linspace(0,1,n);
total_time = 10;
dt = 0.01;
total_time_steps = total_time/dt;
T = 298*ones(n,n); % initial temperature condition
T(:,1) = 400; % K ,temperature left boundary condition
T(end,:) = 600; % K ,temperature top boundary condition
T(:,end) = 800; % K ,temperature right boundary condition
T(1,:) = 900; % K ,temperature bottom boundary condition
% defining the corner temperature for smoother graph
T(1,1) = 500; % k top left boundary condition
T(n,1) = 650; % k bottom left boundary condition
T(1,n) = 700; % k top right boundary condition
T(n,n) = 850; % k bottom right boundary condition
tol = 1e-4;
Alpha = 0.002; % thermal diffusivity
k = Alpha*dt/(dx^2);
CFL = Alpha*dt/(dx^2);
Told1 = T; Told2 = T; Told3 = T; Told4 = T;
T1 = T; T2 = T; T3 = T; T4 = T;
T_time1 = T; T_time2 = T; T_time3 = T; T_time4 = T; % for time loop
% 1.transient explicit, 2.transient implicit jacobi, 3.transient implicit gauss seidel, 4.transient
% implicit SOR
for method = 1:4
% 1.transient explicit
if method ==1
cycle = 0;
tic
for t = 1:total_time_steps
for i = 2:(n-1)
for j = 2:(n-1)
T1(i,j) = Told1(i,j) + k*(Told1(i-1,j) + Told1(i+1,j) + Told1(i,j-1) + Told1(i,j+1)-4*Told1(i,j));
end
end
Told1 = T1;
cycle = cycle + 1;
end
convergence_time1 = toc;
figure(1)
contourf(x,y,T1,'Showtext','on')
colormap(jet)
colorbar
xlabel('X axis')
ylabel('Y axis')
title(sprintf('2D Transient State Analysis (Explicit Method) \nTemperature Distribution over a plate \nNumber of cycles = %d \nTime for convergence = %f seconds',cycle,convergence_time1))
pause(0.02)
% 2.transient implicit jacobi
elseif method == 2
tic
for t = 1:total_time_steps
jacobi_error = 9e9;
jacobi_iter = 0;
while (jacobi_error>tol)
jacobi_iter = jacobi_iter + 1;
for i = 2:(n-1)
for j = 2:(n-1)
T2(i,j) = T_time2(i,j)/(1+4*k) + (k/(1+4*k))*(Told2(i-1,j) + Told2(i+1,j) + Told2(i,j-1) + Told2(i,j+1));
end
end
jacobi_error = max(max(abs(T2-Told2)));
Told2 = T2;
end
T_time2 = T2;
end
convergence_time2 = toc;
figure(2)
contourf(x,y,T2,'Showtext','on')
colormap(jet)
colorbar
xlabel('X axis')
ylabel('Y axis')
title(sprintf('2D Transient State Analysis (Implicit Jacobi Method) \nTemperature Distribution over a plate \nNumber of iterations = %d \nTime for convergence = %f seconds',jacobi_iter,convergence_time2))
pause(0.02)
% 3.transient implicit gauss seidel
elseif method == 3
tic
for t = 1:total_time_steps
GS_error = 9e9;
GS_iter = 0;
while(GS_error>tol)
GS_iter = GS_iter + 1;
for i = 2:(n-1)
for j = 2:(n-1)
T3(i,j) = T_time3(i,j)/(1+4*k) + (k/(1+4*k))*(T3(i-1,j) + Told3(i+1,j) + T3(i,j-1) + Told3(i,j+1));
end
end
GS_error = max(max(abs(Told3-T3)));
Told3 = T3;
end
T_time3 = T3;
end
convergence_time3 = toc;
figure(3)
contourf(x,y,T3,'Showtext','on')
colormap(jet)
colorbar
xlabel('X axis')
ylabel('Y axis')
title(sprintf('2D Transient State Analysis (Implicit Gauss Seidel Method) \nTemperature Distribution over a plate \nNumber of iterations = %d \nTime for convergence = %f seconds',GS_iter,convergence_time3))
pause(0.02)
% 4.transient implicit SOR
elseif method == 4
w = 2/(1+sin(pi*dx)); % relaxation parameter
tic
for t = 1:total_time_steps
SOR_error = 9e9;
SOR_iter = 0;
while (SOR_error>tol)
SOR_iter = SOR_iter + 1;
for i = 2:(n-1)
for j = 2:(n-1)
T4(i,j) =(1-w)*Told4(i,j) + w*(T_time4(i,j))/(1+4*k) + w*(k/(1+4*k))*(T4(i-1,j) + Told4(i+1,j) + T4(i,j-1) + Told4(i,j+1));
end
end
SOR_error = max(max(abs(Told4-T4)));
Told4 = T4;
end
T_time4 = T4;
end
convergence_time4 = toc;
figure(4)
contourf(x,y,T4,'Showtext','on')
colormap(jet)
colorbar
xlabel('X axis')
ylabel('Y axis')
title(sprintf('2D Transient State Analysis (Implicit SOR Method) \nTemperature Distribution over a plate \nNumber of iterations = %d \nTime for convergence = %f seconds',SOR_iter,convergence_time4))
pause(0.02)
end
end
Report
2D Transient State heat conduction equation:
We have to descritize the 2D Transient state heat conduction equation, so that it can be solved numerically in matlab. After descritizing the equation we can solve the equation by explicit and implicit method. In implicit method we can use three different iteration techniques, these are Jacobi method, Gauss Seidel method and Successive over relaxation method.
First we define the number of grid points n=50, number of grid points in x direction is equal to the number of grid points in y direction. The metal plate is a thin square plate of size 1 m. After this we define the grit size dx and dy.
Then we deifne the initial temperature condition to 298 kelvin accross entire nxn meshgrid, by using 'ones' command. And then we define the boundary temperature conditions according to the values given to us. Absoulte error criterio is equal to 1e-4, the error in the solution is compared after each iterations with absolute error criterio and if it is more than the absolute error, then the iterations keeps on going until the value of error is less than the absolute error criteria.
As this is a Transient state analysis, there is a time derivative in the equation, for this we have to define a time loop which will run for the defined number of time steps. We have to define the value of total time for the simulation, time step (dt) and total number of time steps. We also have to take a value for thermal diffusivity alpha = 0.002.
2.1 Explicit method
Formula for explicit method:
The above formula is defined in a time loop. The loop runs for 1000 time steps. Told is updated in each time loop. The initial value of RHS is assumed which gives us the assumed value of T at 'n+1' time interval. For next value of i,j ,T is calculated at the next node. The equations converges in 1000 time steps to give the solution.
In explicit scheme there is simple updating of old value take place and the solution does not depend upon other value at the current level. While implicit scheme contains the information at the current level which requires solution of simaultaneous equations.
2.2 Implicit Jacobi method
Formula for implicit Jacobi method:
The above formula is defined in a convergence loop, which itself is defined in a time loop. Told is updated in convergence loop and T_time values are updated in time loop. The initial values are assumed and the first iteration is done to calculate the T at 'n+1' time interval at the next node. After each convergence loop the difference between previous value of T and calculated value of T is calculated. The error is compared to the tolerance value, the iterations keeps on going until the error value is less than tolarance.
2.3 Implicit Gauss Seidel method
Fromula for implicit Gauss Seidel method:
The above formula is defined in a convergence loop, which itself is defined in a time loop. Told is updated in convergence loop and T_time values are updated in time loop. The above formula is similar to implicit jacobi method, except, when the first iteration is done to calculate the value of T, this updated value of T is used to calculate value of T at next node. The error is compared after each convergence loop.
2.4 Implicit Successive over relaxation (SOR) method
Formula for implicit SOR method:
Formula for relxation parameter:
We used 'tic' and 'toc' command to calculate the time of convergence for explicit and implicit method. After each time loop we plot the result on a contour plot between x,y,T. Contour plot visulises the temperature distribution over the plate and also gives us the number of iterations and time needed for convergence.
Graphs
Conclusion:
- From the above graphs it is clear that explicit method is faster than implicit method to solve 2D transient state heat conduction equation.
- Among implicit method the Guass seidel method took least time for convergence.
- Jacobi method took less time for convergence than SOR but took more time than Gauss seidel method.
- Maxmum time for convergence was taken by SOR method.
- In 2D transient state analysis, while solving implicitly the number of Jacobi, Gauss Seidel and SOR iterations are less than the steady state analysis. But the time of convergence is more, this happens because here thermal diffusivity (Alpha) and time derivative is taken into consideration.