MATLAB - Mid term project - Solving the steady and unsteady 2D heat conduction
Problem Statement: -
For A 2D surface maintained at a certain temperature at its edges. Solve using Iteration Methods(Jacobi, Gauss-Seidel, and S.O.R. using Gauss-Seidel) for Temperature distribution due to Heat Conduction that can be observed for the following States: -
1. Transient State until Steady State is Achieved
2. Steady State
The Given Temperature Maintained at the edges are as follows
`T_("Top") = 600 K, T_("Bottom") = 900 K, T_("Left") = 400K, and T_("Right") = 800K`
Solution: -
1. AIM: -
A. To Solve 2D heat Conduction Equation for Steady and Transient State using Iterative Solver( Jacobi, Gauss-Seidel and Successive Over-Relaxation)
2. Given: -
A. `T_("top") = 600 K`
B. `T_("Bottom") = 900 K`
C. `T_("Left") = 400 K`
D. `T_("Right") = 800K`
3. Assumed: -
A. `"2D Transient Heat Conduction" rarr (delT)/(delt) - alpha((del^2T)/(delx^2) + (del^2T)/(dely^2)) = 0`
B. `"2D Steady Heat Conduction" rarr (del^2T)/(delx^2) + (del^2T)/(dely^2) = 0`
C. Dimensions of Surface `=> "length" = 1 m and "Width" = 1 m`
D. Number of Nodes in each Direction = 51
E. Using the Following notation for `T_(i,j)^n`
where `i,j rarr x and y " positions of Node"`
`n rarr "time Step iteration"`
4. Deriving Numerical Methods: -
A. 2D Steady Heat Conduction: -
For 2D Steady Heat Conduction, `(del^2T)/(delx^2) + (del^2T)/(dely^2) = 0`
First Order Numerical Method,
`=> ((T_(i+1,j)^n - 2T_(i,j)^n + T_(i-1,j)^n)/dx^2) + ((T_(i,j+1)^n - 2T_(i,j)^n + T_(i,j-1)^n)/dy^2) = 0`
`=>2T_(i,j)^n(1/dx^2 + 1/dy^2) = ((T_(i+1,j)^n + T_(i-1,j)^n)/dx^2) + ((T_(i,j+1)^n + T_(i,j-1)^n)/dy^2)`
`=>2T_(i,j)^n((dx^2 + dy^2)/(dx^2dy^2)) = ((T_(i+1,j)^n + T_(i-1,j)^n)/dx^2) + ((T_(i,j+1)^n + T_(i,j-1)^n)/dy^2)`
`=>T_(i,j)^n = 1/2[((T_(i+1,j)^n + T_(i-1,j)^n)/(dx^2 + dy^2))dy^2 + ((T_(i,j+1)^n + T_(i,j-1)^n)/(dx^2 + dy^2))dx^2]`
`"Let "(dy^2)/(dx^2 + dy^2) " be " k_1 and "Let "(dx^2)/(dx^2 + dy^2) " be " k_2`
`=>T_(i,j)^n = 1/2[(T_(i+1,j)^n + T_(i-1,j)^n)k_1 + (T_(i,j+1)^n + T_(i,j-1)^n)k_2]`
Now,
`"For " dx = 1/(51-1) and dy = 1/(51-1)`
`k_1 = (dy = 1/50)^2/((dx = 1/50)^2 + (dy = 1/50)^2) = 1/2`
`k_2 = (dx = 1/50)^2/((dx = 1/50)^2 + (dy = 1/50)^2) = 1/2`
`=>T_(i,j)^n = 1/2[(T_(i+1,j)^n + T_(i-1,j)^n)(k_1 = 1/2) + (T_(i,j+1)^n + T_(i,j-1)^n)(k_2 = 1/2)]`
Hence,
`T_(i,j)^n = (T_(i+1,j)^n + T_(i-1,j)^n)/4 + (T_(i,j+1)^n + T_(i,j-1)^n)/4`
| 2D Steady State Heat Conduction Equation | `T_(i,j)^n = (T_(i+1,j)^n + T_(i-1,j)^n)/4 + (T_(i,j+1)^n + T_(i,j-1)^n)/4` |
B. 2D Transient Heat Conduction: -
For 2D Steady Heat Conduction, `(delT)/(delt)-alpha((del^2T)/(delx^2) + (del^2T)/(dely^2)) = 0`
First Order Numerical Method,
`=> ((T_(i,j)^(n+1) - T_(i,j)^n)/dt)-alpha[((T_(i+1,j)^(n+1) - 2T_(i,j)^(n+1) + T_(i-1,j)^(n+1))/dx^2) + ((T_(i,j+1)^(n+1) - 2T_(i,j)^(n+1) + T_(i,j-1)^(n+1))/dy^2)] = 0`
`=> T_(i,j)^(n+1) - T_(i,j)^n = (dt)alpha[((T_(i+1,j)^(n+1) - 2T_(i,j)^(n+1) + T_(i-1,j)^(n+1))/dx^2) + ((T_(i,j+1)^(n+1) - 2T_(i,j)^(n+1) + T_(i,j-1)^(n+1))/dy^2)]`
`=> T_(i,j)^(n+1) = T_(i,j)^n + (dt)alpha[((T_(i+1,j)^(n+1) - 2T_(i,j)^(n+1) + T_(i-1,j)^(n+1))/dx^2) + ((T_(i,j+1)^(n+1) - 2T_(i,j)^(n+1) + T_(i,j-1)^(n+1))/dy^2)]`
`=> T_(i,j)^(n+1)(1 + 2(dt)alpha(1/dx^2+1/dy^2)) = T_(i,j)^n + (dt)alpha[((T_(i+1,j)^(n+1) + T_(i-1,j)^(n+1))/dx^2) + ((T_(i,j+1)^(n+1) + T_(i,j-1)^(n+1))/dy^2)]`
As, `CFL_x = (alpha(dt))/(dx^2)` and `CFL_y = (alpha(dt))/(dy^2)`
For Stable Numerical Method `CFL < 1`
Now,
Let `(alpha(dt))/dx^2 = k_1`
Let `(alpha(dt))/dy^2 = k_2`
Let `1 + 2(dt)alpha(1/dx^2+1/dy^2) = k_3`
`=> k_3 = 1 + 2(k_1 + k_2)`
Hence,
`=> T_(i,j)^(n+1)(k_3) = T_(i,j)^n + k_1(T_(i+1,j)^(n+1) + T_(i-1,j)^(n+1)) + k_2(T_(i,j+1)^(n+1) + T_(i,j-1)^(n+1))`
`=> T_(i,j)^(n+1) = 1/k_3[T_(i,j)^n + k_1(T_(i+1,j)^(n+1) + T_(i-1,j)^(n+1)) + k_2(T_(i,j+1)^(n+1) + T_(i,j-1)^(n+1))]`
| 2D Unsteady State Heat Conduction Equation | `T_(i,j)^(n+1) = 1/k_3[T_(i,j)^n + k_1(T_(i+1,j)^(n+1) + T_(i-1,j)^(n+1)) + k_2(T_(i,j+1)^(n+1) + T_(i,j-1)^(n+1))]` |
5. Explanation Code: -
A. 2D Steady-State Heat Conduction: -
i. Section 1(Defining Variables)
a. Defining Geometry of Surface
`rarr`Defined Length and Width
`rarr`Defined nodes along x and y
`rarr`Defined dx and dy based on nodes and Dimension
`rarr`Created X and Y Vectors with the geometric distance of each node
b. Defining Boundary Condition and Initial Surface Temperature
`rarr`Defined each node of Surface to be at 500 K
`rarr`Defined Top, Bottom, Left, and Right Nodes at a Certain Temperature and the Corners as the average of adjacent nodes.
c. Variables for Iteration
`rarr`Defined Set the tolerance for the Temperature change to be `10^(-5)`
`rarr`Defined initial error = 1
`rarr`Defined iteration counter and set value = 0
d. Method Selection
`rarr`Defined a Variable "Method" whose value is used to select the method for iteration
1 - Jacobi Method, 2 - Gauss-Siedel, and 3 - S.O.R. with Gauss-Siedel(`omega = 1.95`)
ii. Section 2(Computation)
a. Using Jacobi
`rarr` Checked if Method Selected is Jacobi using the "Method" Variable
`rarr` displayed that the iteration in use is Jacobi
`rarr` Created a While loop with condition that the loop terminates after (error for Temperature < Tolerance of Temperature)
`rarr` Applied Jacobi Iteration on each Node for Value of Temperature based on the derived formula.
`rarr` Updated Value of error of Temperature that is maximum error Node calculated from `T_("new") - T_("old")`
b. Using Gauss-Siedel
`rarr` Checked if Method Selected is Gauss-Siedel using the "Method" Variable
`rarr` displayed that the iteration in use is Gauss-Siedel
`rarr` Created a While loop with condition that the loop terminates after (error for Temperature < Tolerance of Temperature)
`rarr` Applied Gauss-Siedel Iteration on each Node for Value of Temperature based on the derived formula.
`rarr` Updated Value of error of Temperature that is maximum error Node calculated from `T_("new") - T_("old")`
c. Using S.O.R. using Gauss-Siedel
`rarr` Checked if Method Selected is S.O.R. using Gauss-Siedel using the "Method" Variable
`rarr` displayed that the iteration in use is S.O.R. using Gauss-Siedel
`rarr` Created a While loop with condition that the loop terminates after (error for Temperature < Tolerance of Temperature)
`rarr` Applied S.O.R. using Gauss-Siedel Iteration on each Node for Value of Temperature based on the derived formula.
`rarr` Updated Value of error of Temperature that is maximum error Node calculated from `T_("new") - T_("old")`
iii. Section 3(OUTPUT and Plotting)
`rarr` displaying the Total number of iterations
`rarr` Displaying the latest updated error
`rarr` Plotted the Graph of Temperature Distribution
B. 2D UnSteady-State Heat Conduction: -
i. Section 1(Defining Variables)
a. Defining Geometry of Surface
`rarr`Defined Length and Width
`rarr`Defined nodes along x and y
`rarr`Defined dx and dy based on nodes and Dimension
`rarr`Created X and Y Vectors with the geometric distance of each node
b. Defining Boundary Condition and Initial Surface Temperature
`rarr`Defined each node of Surface to be at 500 K
`rarr`Defined Top, Bottom, Left, and Right Nodes at a Certain Temperature and the Corners as the average of adjacent nodes.
c. Variables for Iteration
`rarr`Defined Set the tolerance for the Temperature change to be `10^(-4)`
`rarr`Defined initial error = 1
`rarr`Defined iteration counter and set value = 0
d. Time Step Variable
`rarr`Defined value of CFL = 0.9
`rarr`Defined `alpha = 1.2`
`rarr`Calculated dt
`rarr`Defined k1, k2 and k3
`rarr`Defined total time Step and a Vector to define each time step
e. Method Selection
`rarr`Defined a Variable "Method" whose value is used to select the method for iteration
1 - Jacobi Method, 2 - Gauss-Siedel, and 3 - S.O.R. with Gauss-Siedel(`omega = 1.44`)
ii. Section 2(Computation)
a. Time Marching Loop
`rarr` Defined a for loop from 1 to total time steps
`rarr` Recorded current value of T in `T_t` Matrix
`rarr` Created a While Loop based on error value and Tolerance Value
`rarr` Selected Each node and based on the value of "Method" Selected each iteration and applied formula based on derivation.
`rarr` Within each time Step Plotted the Change in temperature
`rarr` Placed a Break Command which is based on Tolerance of Temperature and breaks the loop when Stability is Reached.
iii. Section 3(OUTPUT)
`rarr` Displayed The Type of Iteration under use
`rarr` Displayed the Total number of Iterations and Total Time taken to reach Stability.
6. Code: -
A. 2D Steady-State Heat Conduction: -
% 2D Steady State Heat Conduction using Iterative Methods
clear all
close all
clc
%% Defining Variables
% Geometry of Surface
L = 1;
W = 1;
Node_x = 51;
Node_y = 51;
dx = L/(Node_x-1);
dy = W/(Node_y-1);
X = [0:dx:L];
Y = [0:dy:W];
% Boundary Conditions and Initial Surface Temperature
T = 500*ones(Node_x,Node_y);
T_top = 600;
T_bottom = 900;
T_left = 400;
T_right = 800;
T(2:(Node_x -1),Node_y) = T_top;
T(2:(Node_x - 1),1) = T_bottom;
T(1,2:(Node_y - 1)) = T_left;
T(Node_x,2:(Node_y -1)) = T_right;
T(1,1) = 0.5*(T(1,2)+T(2,1));
T(Node_x,Node_y) = 0.5*(T(Node_x,(Node_y-1))+T((Node_x-1),Node_y));
T(Node_x,1) = 0.5*(T(Node_x,2)+T((Node_x-1),1));
T(1,Node_y) = 0.5*(T(1,(Node_y-1))+T(2,Node_y));
% Variables for Iterations
Tol_T = 1e-5;
err_T = 1;
iter = 0;
% Method Selection
% 1 - Jacobi,
% 2 - Gauss-Siedel
% 3 - S.O.R. with Gauss-Siedel
Method = 2;
%% Computation
% Using Jacobi
if Method == 1
disp('Iteration Method in Use is Jacobi')
while (err_T > Tol_T)
iter = iter+1;
T_old = T;
for i = 2:(Node_x - 1)
for j = 2:(Node_y - 1)
T(i,j) = 0.25*(T_old(i+1,j) + T_old(i-1,j) + T_old(i,j+1) + T_old(i,j-1));
end
end
err_T = max(abs(T - T_old),[],'all');
end
end
% Using Gauss-Siedel
if Method == 2
disp('Iteration Method in Use is Gauss-Siedel')
while (err_T > Tol_T)
iter = iter+1;
T_old = T;
for i = 2:(Node_x - 1)
for j = 2:(Node_y - 1)
T(i,j) = 0.25*(T(i+1,j) + T(i-1,j) + T(i,j+1) + T(i,j-1));
end
end
err_T = max(abs(T - T_old),[],'all');
end
end
% Using S.O.R. and Gauss-Siedel
if Method == 3
disp('Iteration Method in Use is S.O.R.')
omega_val = 1.95;
disp(['Value of Omega = ',num2str(omega_val)])
while (err_T > Tol_T)
iter = iter+1;
T_old = T;
for i = 2:(Node_x - 1)
for j = 2:(Node_y - 1)
T(i,j) = T(i,j)*(1-omega_val) + (omega_val)*0.25*(T(i+1,j) + T(i-1,j) + T(i,j+1) + T(i,j-1));
end
end
err_T = max(abs(T - T_old),[],'all');
end
end
%% OUTPUT and Plotting
disp(['Total Iteration = ',num2str(iter)])
disp(['Error = ',num2str(err_T)])
colormap(jet);
contourf(X,Y,T');
colorbar
title('Temperature Distribution','FontSize',16)
xlabel('Length of Surface(m)','Color',[0.15 0.5 0.5],'FontSize',14)
ylabel('Width of Surface(m)','Color',[0.15 0.5 0.5],'FontSize',14)B. 2D UnSteady-State Heat Conduction: -
% 2D Un-Steady State Heat Conduction using Iterative Methods
clear all
close all
clc
%% Defining Variables
% Geometry of Surface
L = 1;
W = 1;
Node_x = 51;
Node_y = 51;
dx = L/(Node_x-1);
dy = W/(Node_y-1);
X = [0:dx:L];
Y = [0:dy:W];
% Boundary Conditions and Initial Guess
T = 500*ones(Node_x,Node_y);
T_top = 600;
T_bottom = 900;
T_left = 400;
T_right = 800;
T(2:(Node_x -1),Node_y) = T_top;
T(2:(Node_x - 1),1) = T_bottom;
T(1,2:(Node_y - 1)) = T_left;
T(Node_x,2:(Node_y -1)) = T_right;
T(1,1) = 0.5*(T(1,2)+T(2,1));
T(Node_x,Node_y) = 0.5*(T(Node_x,(Node_y-1))+T((Node_x-1),Node_y));
T(Node_x,1) = 0.5*(T(Node_x,2)+T((Node_x-1),1));
T(1,Node_y) = 0.5*(T(1,(Node_y-1))+T(2,Node_y));
% Variables for Iterations
Tol_T = 1e-4;
err_T = 1;
iter = 0;
% Time Step Variables
CFL = 0.9;
alp = 1.2;
dt = CFL*(dx^2)/alp;
k1 = CFL;
k2 = CFL;
k3 = 1 + 2*k1 + 2*k2;
t_steps = 1500;
total_t = [0:dt:(t_steps*dt)];
% Method Selection
% 1 - Jacobi,
% 2 - Gauss-Siedel
% 3 - S.O.R. with Gauss-Siedel
Method = 3;
omega_val = 1.44;
%% Computation
% Time Marching Loop
for t_count = 1:t_steps
T_t = T;
% Iteration using Gauss-Siedel
err_T = 1;
while (err_T > Tol_T)
iter = iter+1;
T_old = T;
for i = 2:(Node_x - 1)
for j = 2:(Node_y - 1)
if Method == 1
T(i,j) = (T_t(i,j) + k1*(T_old(i+1,j) + T_old(i-1,j) + T_old(i,j+1) + T_old(i,j-1)))/k3;
end
if Method == 2
T(i,j) = (T_t(i,j) + k1*(T(i+1,j) + T(i-1,j) + T(i,j+1) + T(i,j-1)))/k3;
end
if Method == 3
T(i,j) = T(i,j)*(1 - omega_val) + (omega_val)*(T_t(i,j) + k1*(T(i+1,j) + T(i-1,j) + T(i,j+1) + T(i,j-1)))/k3;
end
end
end
err_T = max(abs(T - T_old),[],'all');
end
% Plotting after each Iteration
colormap(jet);
contourf(X,Y,T');
colorbar
title(['Temperature Distribution at time = ',num2str(total_t(t_count)),'sec'],'FontSize',12)
xlabel('Length of Surface(m)','Color',[0.15 0.5 0.5],'FontSize',14)
ylabel('Width of Surface(m)','Color',[0.15 0.5 0.5],'FontSize',14)
pause(0.05)
err_T_t = max(abs(T - T_t),[],'all');
% Breaking Time Marching if Change become less than Tolerance
if err_T_t < Tol_T
break
end
end
%% OUTPUTS
if Method == 1
disp('The Iteration method in use was Jacobi')
end
if Method == 2
disp('The Iteration method in use was Gauss-Siedel')
end
if Method == 3
disp('The Iteration method in use was S.O.R. with Gauss-Siedel')
disp(['Value of Omega for S.O.R. = ',num2str(omega_val)])
end
disp(['Total Iteration = ',num2str(iter)])
disp(['Total Time for reaching Stability = ',num2str(total_t(t_count)),'sec'])7. OUTPUTS and Plots Generated: -
A. 2D Steady-State Heat Conduction: -
i. Plot Generated: -
B. 2D UnSteady-State Heat Conduction: -
i. Initial Temperature Distribution
ii. Final Temperature Distribution
