Week 5.1 - Mid term project - Solving the steady and unsteady 2D heat conduction problem
%% ========== IMPLICIT SCHEME USING JACOBI METHOD ========== %%
% Solving 2D Heat Conduction Problem using Jacobi Method (Steady)
% Equation = d^2T/dx^2 + d^2T/dy^2
% Top Boundary = 600 K
% Bottom Boundary = 900 K
% Left Boundary = 400 K
% Right Boundary = 800 K
function T = jacobi(nx,ny)
% Input values
% value for nx and ny will be defined in the main function
% ASSIGNING BOUNDARY CONDITIONS
T = 300*ones(nx,ny);
% Assigning Top BC
T(1,:) = 600;
% Assigning Left BC
T(:,1) = 400;
% Assigning Bottom BC
T(nx,:)= 900;
% Assigning Right BC
T(:,nx)= 800;
% Defining BC for corner points
% For Left Top corner
T(1,1)= (T(2,2)+ T(1,2))/2;
% For Left Bottom corner
T(end,1)= (T(end-1,1) + T(end,2))/2;
% For Right Top corner
T(1,end) = (T(1,end-1) + T(2,end))/2;
% For Right Bottom corner
T(nx,end)= (T(nx,end-1) + T(end-1,nx))/2;
% Creating a copy of T
x= linspace(0,1,nx);
dx= x(2)-x(1);
y= linspace(0,1,ny);
dy= y(2)-y(1);
k= 2*(((dx.^2)+(dy.^2))/((dx.^2)*(dy.^2)));
% JACOBI METHOD
error = 9e9;
tol = 1e-4;
Told = T;
jacobi_iteration = 1;
while (error>tol)
for i=2:(nx-1)
for j= 2:(ny-1)
T(i,j)=((Told(i-1,j)+Told(i+1,j))./(dx^2*k) ) + ((Told(i,j-1)+Told(i,j+1))./(dy^2*k));
end
end
error = max(max(abs(Told-T)));
Told= T;
jacobi_iteration = jacobi_iteration + 1;
end
%% ========== IMPLICIT SCHEME USING GAUSS SEIDEL METHOD ========== %%
% Solving 2D Heat Conduction Problem using Gauss-Seidel Method (Steady)
% Equation = d^2T/dx^2 + d^2T/dy^2
% Top Boundary = 600 K
% Bottom Boundary = 900 K
% Left Boundary = 400 K
% Right Boundary = 800 K
function T = GS(nx,ny)
% Input values
% value for nx and ny will be defined in the main function
% ASSIGNING BOUNDARY CONDITIONS
T = 300*ones(nx,ny);
% Assigning Top BC
T(1,:) = 600;
% Assigning Left BC
T(:,1) = 400;
% Assigning Bottom BC
T(nx,:)= 900;
% Assigning Right BC
T(:,nx)= 800;
% Defining BC for corner points
% For Left Top corner
T(1,1)= (T(2,2)+ T(1,2))/2;
% For Left Bottom corner
T(end,1)= (T(end-1,1) + T(end,2))/2;
% For Right Top corner
T(1,end) = (T(1,end-1) + T(2,end))/2;
% For Right Bottom corner
T(nx,end)= (T(nx,end-1) + T(end-1,nx))/2;
% Creating a copy of T
x= linspace(0,1,nx);
dx= x(2)-x(1);
y= linspace(0,1,ny);
dy= y(2)-y(1);
k= 2*(((dx.^2)+(dy.^2))/((dx.^2)*(dy.^2)));
% GAUSS SEIDEL METHOD
error = 9e9;
tol = 1e-4;
Told = T;
GS_iteration = 1;
while (error>tol)
for i=2:(nx-1)
for j= 2:(ny-1)
T(i,j)=((T(i-1,j)+Told(i+1,j))./(dx^2*k) ) + ((T(i,j-1)+Told(i,j+1))./(dy^2*k));
end
error = max(max(abs(Told-T)));
Told= T;
GS_iteration = GS_iteration + 1;
end
end
%% ========== IMPLICIT SCHEME USING SUCCESSIVE OVER RELAXATION METHOD ========== %%
% Solving 2D Heat Conduction Problem using Successive Over-relaxation
% Method (SOR) Method (Steady)
% Equation = d^2T/dx^2 + d^2T/dy^2
% Top Boundary = 600 K
% Bottom Boundary = 900 K
% Left Boundary = 400 K
% Right Boundary = 800 K
function T = SOR(nx,ny)
% Input values
% value for nx and ny will be defined in the main function
% ASSIGNING BOUNDARY CONDITIONS
T = 300*ones(nx,ny);
% Assigning Top BC
T(1,:) = 600;
% Assigning Left BC
T(:,1) = 400;
% Assigning Bottom BC
T(nx,:)= 900;
% Assigning Right BC
T(:,nx)= 800;
% Defining BC for corner points
% For Left Top corner
T(1,1)= (T(2,2)+ T(1,2))/2;
% For Left Bottom corner
T(end,1)= (T(end-1,1) + T(end,2))/2;
% For Right Top corner
T(1,end) = (T(1,end-1) + T(2,end))/2;
% For Right Bottom corner
T(nx,end)= (T(nx,end-1) + T(end-1,nx))/2;
% Creating a copy of T
x= linspace(0,1,nx);
dx= x(2)-x(1);
y= linspace(0,1,ny);
dy= y(2)-y(1);
k= 2*(((dx.^2)+(dy.^2))/((dx.^2)*(dy.^2)));
% SUCCESSIVE OVER RELAXATION METHOD
error = 9e9;
tol = 1e-4;
Told = T;
alfa = 1.1;
sor_iteration = 1;
while (error>tol)
for i=2:(nx-1)
for j= 2:(ny-1)
T(i,j)= (1-alfa)*Told(i,j) + (alfa*(T(i-1,j)+Told(i+1,j))./(dx^2*k)) + (alfa*(T(i,j-1)+Told(i,j+1))./(dy^2*k));
end
error = max(max(abs(Told-T)));
Told= T;
sor_iteration = sor_iteration + 1;
end
end
%% ========== MAIN CODE FOR JACOBI, GS, SOR ========== %%
a = input('Input value for nx: ');
nx=a;
ny=a;
Temp1 = jacobi(nx,ny);
Temp2 = GS(nx,ny);
Temp3 = SOR(nx,ny);
x= linspace(0,1,nx);
y= linspace(0,1,ny);
%% code for making the contour plot
[XX, YY]= meshgrid(x,y);
figure(1)
contourf(XX, -YY, Temp1)
colormap(jet)
colorbar
[XX, YY]= meshgrid(x,y);
figure(2)
contourf(XX, -YY, Temp2)
colormap(jet)
colorbar
[XX, YY]= meshgrid(x,y);
figure(3)
contourf(XX, -YY, Temp3)
colormap(jet)
colorbar
%% ========== THEORY ========= %%
- When using IMPLICIT SCHEME, one cannot solve equation like traditional method.
- While solving governing equation via implicit scheme, we only have one unkown unlike explicit method.
- In order to solve such equation, we use iterative method. Also know as fixed point method.
- Iterative method discussed here are: 1) Jacobi 2) Gauss Seidel 3) Successive over relaxation.
- The given problem is solved for 2 dimensional head conduction eqution for both steady and transient state.
%% ========== CODE EXPLANATION FOR STEADY STATE IMPLICIT SCHEME ========== %%
- An array is created with initial value of all the nodes as 300K.
- In the section 'ASSIGNING THE BOUNDARY CONDITION', the boundary conditions are applied to all the nodes as specified in the problem statement.
- A copy of this is stored in variable Told.
- An error term and tolerance value is assigned. The loop will run until the error drops down or is less than tolerance value.
- It will ask user for 'nx', which is used to create the initial array.
- In the main code the functions are called and contour is plotted using contourf and meshgrid command.
%% ===== DIFFERENCE BETWEEN JACOBI, GAUSS-SEIDEL AND SUCCESSIVE OVER RELAXATION SCHEME ===== %%
- JACOBI: In this method all the value of temperature are taken from the the previous iteration.
- GAUSS-SEIDEL: In this method the value of temperature is taken from the same iteration for the node which is before 'ith' node. Meaning, if we are computing for 4th node, the value for the 3rd node can be taken from the same iteration, whereas for 5th node the values cannot be taken from the same iteration since it has not been computed yet.
- SUCCESSIVE OVER RELAXATION: In this method we take value from previous iteration and combine it with either Jacobi or Gauss-Seidel method.
%% ========== DERIVATION ========== %%
Consider 1-D Heat Conduction equation:
`del`^2T/`del`x^2
Discitize this equaiton:
(T(i+1) - 2T(i) + T(i-1))/2
If we consider 5 node points :
-----x-----x-----x-----x-----x----- here x = nodes
1 2 3 4 5
Therefore,
For i = 2, (T3 - 2T2 + T1)/`triangle`x^2 = 0
For i = 3, (T4 - 2T3 - T2)/`triangle`x^2 = 0
for i = 4, (T5 - 2T4 + T3)/`triangle`x^2 = 0
Now simplify,
-2T2+T1+T3 = 0
-2T3+T4+T5 = 0
-2T4+T3+T5 = 0
On further simplification,
T2(k+1) = ((T1+T3)/2)k
T3(k+1) = ((T2+T4)/2)k
T4(k+1) = ((T3+T5)/2)k
k+1 represents '(k+1)th' iteration and k represents 'kth' iteration.
Note: It is not a part of equation, it for understanding purposes.
Considering RHS, if we use the values of temperature at 'kth' iteration, then we are using Jacobi Method.
If we are using '(k+1)th' iteration for T1, T2, T3 of RHS and 'kth' iteration for T3, T4, T5, for respective T2, T3, T4 on LHS, then we are using Gauss-Seidel Method. As shown below:
T2(k+1) = ((T1(k+1)+T3(k))/2)
T3(k+1) = ((T2(k+1)+T4(k))/2)
T4(k+1) = ((T3(k+1)+T5(k))/2)
For Successive Over Relaxation Method, we use previous iteration value with combination of either Jacobi or Gauss-Seidel method. As shown below:
Considering Node 3,
T3(k+1) = T3(k) + alfa*(Tgs - T3(k))*(number of iteration) %% If Gauss_Seidel is used.
T3(k+1) = T3(k) + alfa*(Tjack - T3(k))*(number of iteration) %% If Jacobi is used.
Here Tgs = ((T2(k+1) + T4(k))/2); and Tjack = ((T2(k)+T4(k))/2).
==============================================================================
%% ========== EXPLICIT TRANSIENT METHOD ========= %%
% Solving 2D Heat Conduction Problem (Transient Explicit)
% Equation = dT/dt + alfa*(d^2T/dx^2 + d^2T/dy^2)
% Top Boundary = 600 K
% Bottom Boundary = 900 K
% Left Boundary = 400 K
% Right Boundary = 800 K
function T = two_D_HC_transient_explicit(nx,ny,nt,k1,k2)
% Input values
% value for nx and ny will be defined in the main function
% ASSIGNING BOUNDARY CONDITIONS
T = 300*ones(nx,ny);
% Assigning Top BC
T(1,:) = 600;
% Assigning Left BC
T(:,1) = 400;
% Assigning Bottom BC
T(nx,:)= 900;
% Assigning Right BC
T(:,nx)= 800;
% Defining BC for corner points
% For Left Top corner
T(1,1)= (T(2,2)+ T(1,2))/2;
% For Left Bottom corner
T(end,1)= (T(end-1,1) + T(end,2))/2;
% For Right Top corner
T(1,end) = (T(1,end-1) + T(2,end))/2;
% For Right Bottom corner
T(nx,end)= (T(nx,end-1) + T(end-1,nx))/2;
% Creating a copy of T
x= linspace(0,1,nx);
dx= x(2)-x(1);
y= linspace(0,1,ny);
dy= y(2)-y(1);
Told = T;
% Time Loop
for k= 1:nt
% Space Loop
for i=2:(nx-1)
for j= 2:(ny-1)
T(i,j)= Told(i,j) + (k1*(Told(i-1,j) - (2*Told(i,j)) + Told(i+1,j))) + (k2*(Told(i,j-1) - (2*Told(i,j)) + Told(i,j+1)));
end
%error = max(max(abs(Told-T)))
Told= T;
end
% Plotting
[XX, YY]= meshgrid(x,y);
figure(1)
contourf(XX, -YY, T)
colormap(jet)
colorbar
xlabel('X')
ylabel('Y')
title('Temperature Contour')
end
%% ========= MAIN CODE FOR EXPLICIT TRANSIENT ========= %%
% Main code (Function calling for explicit)
nx = input('Input value for nx: ');
ny = nx;
% Calculating number of time step
% Considering value of CFL as 0.4 %%==== CFL = alfa*dt/dx^2 or alfa*dt/dy^2
CFL = 0.3;
alfa = 1.5;
x = linspace(0,1,nx);
dx = x(2)-x(1);
y = linspace(0,1,ny);
dy = y(2)-y(1);
dt = (CFL*(dx^2))/(alfa);
run_time = input('Input Value for run time: ');
nt = run_time/dt;
nt = round(nt);
k1 = (alfa*dt)./(dx^2);
k2 = (alfa*dt)./(dy^2);
Temperature_explicit = two_D_HC_transient_explicit(nx,ny,nt,k1,k2);
%% ========== CODE EXPLANATION ========== %%
- An array is created based on users input. The value at each location is 300K.
- In section 'ASSIGNING THE BOUNDARY CONDITION', the boundary conditions are applied as specified in the problem.
- Since this is 2D problem, two indices have been assigned 'i' and 'j' in order to account spatial behaviour.
- In order to make it transient, time loop is introduced inside which the spatial loop runs.
- In order to run this code, user must specift the value of 'nx'; which decides the size of array and value of run time must also be specified.
- Values of CFL and alfa are predefined (although one can change and see the effect).
- With the help of CFL and alfa value one can calculate step size 'dt'.
- Once 'dt' is calculated, and having run time specified by the user, one can calculate the value of number of time step.
- k1 and k2 are constants. (also refered as CFL).
- value 'nx', 'ny' (nx = ny), 'nt', 'k1' ,'k2' are supplied to the function to calculate or capture the transient behaviour.
%% ========== EQUATION FOR TRANSIENT EXPLICIT METHOD ========== %%
Tp(n+1) = Tp(n) - alfa*dt*[((Tl - 2Tp + Tr )/`triangle`x^2) + ((Tt - 2Tp + Tb)/`triangle`y^2)]
Tp = current node
Tl = left node
Tr = right node
Tt = top node
Tb = bottom node
NOTE: 1) Since the derivative of temperature is elliptical in nature, we have to take information from all nodes.
2) n+1 is the current time interval and n is previous time interval.
==============================================================================
%% ========== IMPLICIT TRANSIENT USING JACOBI METHOD ========== %%
% Solving 2D Heat Conduction Problem using Transient Iterative Method (Jacobi)
% Equation = dT/dt + alfa*(d^2T/dx^2 + d^2T/dy^2)
% Top Boundary = 600 K
% Bottom Boundary = 900 K
% Left Boundary = 400 K
% Right Boundary = 800 K
function T = two_D_HC_transient_implicit_jacobi(nx,ny,nt,k1,k2)
% Input values
% value for nx and ny will be defined in the main function
% ASSIGNING BOUNDARY CONDITIONS
T = 300*ones(nx,ny);
% Assigning Top BC
T(1,:) = 600;
% Assigning Left BC
T(:,1) = 400;
% Assigning Bottom BC
T(nx,:)= 900;
% Assigning Right BC
T(:,nx)= 800;
% Defining BC for corner points
% For Left Top corner
T(1,1)= (T(2,2)+ T(1,2))/2;
% For Left Bottom corner
T(end,1)= (T(end-1,1) + T(end,2))/2;
% For Right Top corner
T(1,end) = (T(1,end-1) + T(2,end))/2;
% For Right Bottom corner
T(nx,end)= (T(nx,end-1) + T(end-1,nx))/2;
% Creating a copy of T
x= linspace(0,1,nx);
dx= x(2)-x(1);
y= linspace(0,1,ny);
dy= y(2)-y(1);
% JACOBI METHOD
tol = 1e-4;
Told = T;
T_previous = T;
term1 = 1/(1+2*k1+2*k2);
term2 = k1*(term1);
term3 = k2*(term1);
jacobi_iteration = 1;
% Time Loop
for a= 1:nt
error = 9e9;
while (error>tol)
% Space Loop
for i= 2:(nx-1)
for j= 2:(ny-1)
H = (Told(i-1,j) + Told(i+1,j));
V = (Told(i,j-1) + Told(i,j+1));
T(i,j)= (term1*(T_previous(i,j)))+(H*term2)+(V*term3);
end
end
error = max(max(abs(Told-T)))
Told= T;
jacobi_iteration = jacobi_iteration + 1;
end
T_previous = T;
% Plotting
[XX, YY]= meshgrid(x,y);
figure(1)
contourf(XX, -YY, T)
colormap(jet)
colorbar
xlabel('X')
ylabel('Y')
title('Temperature Contour')
end
%% ========== IMPLICIT TRANSIENT USING GAUSS-SEIDEL METHOD ========== %%
% Solving 2D Heat Conduction Problem using Transient Iterative Method (GAUSS SEIDEL)
% Equation = dT/dt + alfa*(d^2T/dx^2 + d^2T/dy^2)
% Top Boundary = 600 K
% Bottom Boundary = 900 K
% Left Boundary = 400 K
% Right Boundary = 800 K
function T = two_D_HC_transient_implicit_GS(nx,ny,nt,k1,k2)
% Input values
% value for nx and ny will be defined in the main function
% ASSIGNING BOUNDARY CONDITIONS
T = 300*ones(nx,ny);
% Assigning Top BC
T(1,:) = 600;
% Assigning Left BC
T(:,1) = 400;
% Assigning Bottom BC
T(nx,:)= 900;
% Assigning Right BC
T(:,nx)= 800;
% Defining BC for corner points
% For Left Top corner
T(1,1)= (T(2,2)+ T(1,2))/2;
% For Left Bottom corner
T(end,1)= (T(end-1,1) + T(end,2))/2;
% For Right Top corner
T(1,end) = (T(1,end-1) + T(2,end))/2;
% For Right Bottom corner
T(nx,end)= (T(nx,end-1) + T(end-1,nx))/2;
% Creating a copy of T
x= linspace(0,1,nx);
dx= x(2)-x(1);
y= linspace(0,1,ny);
dy= y(2)-y(1);
% GAUSS SEIDEL METHOD
tol = 1e-4;
Told = T;
T_previous = T;
term1 = 1/(1+2*k1+2*k2);
term2 = k1*(term1);
term3 = k2*(term1);
GS_iteration = 1;
% Time Loop
for a= 1:nt
error = 9e9;
while (error>tol)
% Space Loop
for i= 2:(nx-1)
for j= 2:(ny-1)
H = (T(i-1,j) + Told(i+1,j));
V = (T(i,j-1) + Told(i,j+1));
T(i,j)= (term1*(T_previous(i,j)))+(H*term2)+(V*term3);
end
end
error = max(max(abs(Told-T)))
Told= T;
GS_iteration = GS_iteration + 1;
end
T_previous = T;
% Plotting
[XX, YY]= meshgrid(x,y);
figure(2)
contourf(XX, -YY, T)
colormap(jet)
colorbar
xlabel('X')
ylabel('Y')
title('Temperature Contour')
end
%% ========== IMPLICIT TRANSIENT USING SUCCESSIVE OVER RELAXATION METHOD ========== %%
% Solving 2D Heat Conduction Problem using Transient Iterative Method (SUCCESSIVE OVER RELAXATION)
% Equation = dT/dt + alfa*(d^2T/dx^2 + d^2T/dy^2)
% Top Boundary = 600 K
% Bottom Boundary = 900 K
% Left Boundary = 400 K
% Right Boundary = 800 K
function T = two_D_HC_transient_implicit_SOR(nx,ny,nt,k1,k2)
% Input values
% value for nx and ny will be defined in the main function
% ASSIGNING BOUNDARY CONDITIONS
T = 300*ones(nx,ny);
% Assigning Top BC
T(1,:) = 600;
% Assigning Left BC
T(:,1) = 400;
% Assigning Bottom BC
T(nx,:)= 900;
% Assigning Right BC
T(:,nx)= 800;
% Defining BC for corner points
% For Left Top corner
T(1,1)= (T(2,2)+ T(1,2))/2;
% For Left Bottom corner
T(end,1)= (T(end-1,1) + T(end,2))/2;
% For Right Top corner
T(1,end) = (T(1,end-1) + T(2,end))/2;
% For Right Bottom corner
T(nx,end)= (T(nx,end-1) + T(end-1,nx))/2;
% Creating a copy of T
x= linspace(0,1,nx);
dx= x(2)-x(1);
y= linspace(0,1,ny);
dy= y(2)-y(1);
% SUCCESSIVE OVER RELAXATION METHOD
tol = 1e-4;
Told = T;
T_previous = T;
term1 = 1/(1+2*k1+2*k2);
term2 = k1*(term1);
term3 = k2*(term1);
omega = 0.99;
SOR_iteration = 1;
% Time Loop
for a= 1:nt
error = 9e9;
while (error>tol)
% Space Loop
for i= 2:(nx-1)
for j= 2:(ny-1)
H = (T(i-1,j) + Told(i+1,j));
V = (T(i,j-1) + Told(i,j+1));
T(i,j)= ((1-omega)*Told(i,j)) + (omega*(term1*(T_previous(i,j)))+(H*term2)+(V*term3));
end
end
error = max(max(abs(Told-T)))
Told= T;
SOR_iteration = SOR_iteration + 1;
end
T_previous = T;
% Plotting
[XX, YY]= meshgrid(x,y);
figure(3)
contourf(XX, -YY, T)
colormap(jet)
colorbar
xlabel('X')
ylabel('Y')
title('Temperature Contour')
end
%% ========== MAIN CODE FOR CALLING ALL THE IMPLICIT FUNCTIONS ==========%%
% Main code (Function calling for implicit)
nx = input('Input value for nx: ');
ny = nx;
% Calculating number of time step
% Considering value of CFL as 0.4 %%==== CFL = alfa*dt/dx^2 or
% alfa*dt/dx^2
CFL = 0.4;
alfa = 1.5;
x = linspace(0,1,nx);
dx = x(2)-x(1);
y = linspace(0,1,ny);
dy = y(2)-y(1);
dt = (CFL*(dx^2))/(alfa);
run_time = input('Input Value for run time: ');
nt = run_time/dt;
nt = round(nt);
k1 = (alfa*dt)./(dx^2);
k2 = (alfa*dt)./(dy^2);
Temperature_implicit_jacobi = two_D_HC_transient_implicit_jacobi(nx,ny,nt,k1,k2);
Temperature_implicit_GS = two_D_HC_transient_implicit_GS(nx,ny,nt,k1,k2);
Temperature_implicit_SOR = two_D_HC_transient_implicit_SOR(nx,ny,nt,k1,k2);
%%========== CODE EXPLANATION ========== %%
- An array is created based on users input. The value at each location is 300K.
- In section 'ASSIGNING THE BOUNDARY CONDITION', the boundary conditions are applied as specified in the problem.
- Since this is 2D problem, two indices have been assigned 'i' and 'j' in order to account spatial behaviour.
- In order to make it transient, time loop is introduced inside which the spatial loop runs.
- For each iteration it makes sure that the error value is less than tolerance value.
- In order to run this code, user must specift the value of 'nx'; which decides the size of array and value of run time must also be specified.
- Values of CFL and alfa are predefined (although one can change and see the effect).
- With the help of CFL and alfa value one can calculate step size 'dt'.
- Once 'dt' is calculated, and having run time specified by the user, one can calculate the value of number of time step.
- k1 and k2 are constants. (also refered as CFL).
- value 'nx', 'ny' (nx = ny), 'nt', 'k1' ,'k2' are supplied to the function to calculate or capture the transient behaviour.
%% ========== EQUATION FOR IMPLICIT TRANSIENT METHOD ========= %%
Tp(n+1) =((Tprevious(n))/(1+2*k1+2*k2)) + (1/1+2*k1+2*k2)*[(k1*(Tl(n+1)+Tr(n+1)) + (k2*(Tt(n+1) + Tb(n+1)))))]
Tp = current node
Tl = left node
Tr = right node
Tt = top node
Tb = bottom node
Tprevious = Teperaaturer from previous time interval
NOTE: 1) Since the derivative of temperature is elliptical in nature, we have to take information from all nodes.
2) n+1 is the current time interval and n is previous time interval.