Week 5.1 - Mid term project - Solving the steady and unsteady 2D heat conduction problem
Aim - Solving the steady and unsteady 2D heat conduction problem.
Objective - Steady state analysis || Iterative solver - Jacobi, Gauss -seidel, Relaxation methods
Theory -
Solving this by
1. Jacobi
2. Gauss-seidel
3. Successive Over Relaxation SOR
Expression -
1. Steady State equation:
Under these conditions, the rate of change of temperature remains 0. Therefore in this condition time, the derivative term is equated to zero.
`(d^2T)/dx^2 + (d^2T)/dy^2 = 0`
2.Transient State equation:
Under these conditions, both terms are considered therefore the equation is as follows:
`dT/dt+alpha((d^2T)/dx^2 + (d^2T)/dy^2) = 0`
Explicit and implicit methods are approaches used in numerical analysis for obtaining numerical approximations to the solutions of time-dependent ordinary and PDEs, as is required in computer simulations of physical process. Explicit methods calculate the state of a system at a later time from the state of the system at the current time, while implicit methods find a solution by solving an equation involving both the current state of the system and the later one.
Jacobi :
The unknowns on the right-hand side are assumed to have an initial value, i.e. T at 'n+1' at all nodes is assumed. Using the equation, then the LHS is calculated. On the next value of i,j similar method is applied to calculate LHS. Thus in the first iteration, value of T for time 'n+1' is calculated at all nodes. The difference between T at 'n+1' which is obtained and T at 'n+1' which was assumed previously is taken. If it is greater than a specified tolerance value, iterations are carried out to calculate T at 'n+1'. But in the second iteration, T at 'n+1' is assumed as the value of T which was obtained in the previous iteration.
Gauss siedal :
The gauss siedal method is similar to the Jacobi method, but the only difference is- For the first iteration when the value of T at a particular node is calculated, this updated value at that node is used to calculate the value at the next node. But, the error term is calculated similar to Jacobi method, i.e. T at 'n+1' calculated minus the T at 'n+1' previously assumed.
SOR :
This method can be applied either by using Jacobi or Gauss seidal but the only difference is consideration of factor called omega(scaling factor) which in turn makes the process faster by which the number of iterations required to solve will be decreased
The optimum value of omega is calculated by using the formula : omega = 2/(1+sin(pi*dx))
Program Codes -
% Solving Steady and unsteady 2D heat conduction problem || Steady state analysis
% Iterative solver - Jacobi, Gauss -seidel, Relaxation methods
clear
close
clc
% initial conditions
nx = 10;
ny = 10;
x = linspace(0,1,nx);
y = linspace(0,1,ny);
dx = x(2) - x(1);
dy = y(2) - y(1);
% given data of boundary condition and initial temp condition
t = 300*ones(nx,ny);
t(:,1) = 400;
t(1,:) = 600;
t(:,ny) = 800;
t(nx,:) = 900;
% defining the temperature at corner points
t(1,1) = 500;
t(nx,ny) = 850;
t(1,ny) = 700;
t(nx,1) = 650;
% time step
dt = 1e-2;
told=t;
error = 100;
tol = 1e-4;
k = 2*(dx^2+dy^2)/(dx^2*dy^2);
% Iterative method selection
method = input('Iterative method = ');
% JACOBI ITERATIVE METHOD
if (method==1)
tic;
jacobi = 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))+((told(i,j+1)+told(i,j-1))/(dy^2)))/k;
end
end
error = max(max(abs(t-told)));
told = t;
jacobi = jacobi+1;
time_taken = toc;
end
% Plotting Code
figure(1)
contourf(x,y,t,'ShowText','on');
colormap(jet);
colorbar;
xlabel('X')
ylabel('Y')
text1 = sprintf('number of iteration = %d', jacobi);
text2 = sprintf('simulation time = %fn', time_taken);
title({'Steady-state Conduction Analysis by jacobi'; text1; text2});
end
% GAUSS-SEIDEL ITERATIVE METHOD
if (method==2)
tic;
gauss = 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))+((told(i,j+1)+t(i,j-1))/(dy^2)))/k;
end
end
error = max(max(abs(t-told)));
told = t;
gauss = gauss+1;
time_taken = toc;
end
% Plotting code
figure(2)
contourf(x,y,t,'ShowText','on');
colormap(jet);
colorbar;
xlabel('X')
ylabel('Y')
text1 = sprintf('number of iteration = %d', gauss);
text2 = sprintf('simulation time = %fn', time_taken);
title({'Steady-state Conduction Analysis by Gauss-Seidel'; text1; text2});
end
% SUCCESSIVE OVER RELAXATION(SOR) ITERATIVE METHOD
if (method==3)
tic;
sor = 1;
omega = 1.4; %over relaxation factor using omega
while(error>tol)
for i = 2:(nx-1)
for j = 2:(ny-1)
term1 =(t(i-1,j)+told(i+1,j))/(dx^2);
term2 =(told(i,j+1)+t(i,j-1))/(dy^2);
t(i,j) = ((told(i,j)*(1-omega)))+(omega*(term1+term2))/k;
end
end
error = max(max(abs(t-told)));
told = t;
sor = sor+1;
time_taken = toc;
end
% Plotting code
figure(3)
contourf(x,y,t,'ShowText','on');
colormap(jet);
colorbar;
xlabel('X')
ylabel('Y')
text1 = sprintf('Number of iterations = %d', sor);
text2 = sprintf('Simulation time = %f', time_taken);
text3 = sprintf('Over relaxation factor = %fn', omega);
title({'Steady-state Conduction Analysis by Successive Over Relaxation'; text1; text2; text3});
end
Output Graphs -
Observation -
From above results we can see that the SOR method gives a better in terms of iterations carried out by them and no. of simulations time.
----------------------------------------------------------------------------------------------------------------------------------------------
Objective - Transient state solution (IMPLICIT) || Iterative solver - Jacobi, Gauss -seidel, Relaxation methods
Theory -
Solving the above equation Implicitly :
Implicit Approach is used in numerical analysis for obtaining numerical approximations to the solution of time-dependent ordinary or partial differential equations as in required in computer simulations of a physical process.
Implicit methods find the solution by solving an equation involving both the current state of the system and the later one. Mathematically if Y(t) is the current system state and Y(T+dt) is the state at the later time dt, then, for an explicit method Y(T+dT)=F(Y(t))
Program Codes -
% Solving Steady and Transient 2D heat conduction problem || Transient state solution (IMPLICIT)
% Iterative solver - Jacobi, Gauss -seidel, Relaxation methods
clear
close
clc
% initial conditions
nx = 10;
ny = 10;
nt = 300; % No.of steps
x = linspace(0,1,nx);
y = linspace(0,1,ny);
dx = x(2) - x(1);
dy = y(2) - y(1);
% given data of boundary condition
% initial temp condition
t = 300*ones(nx,ny);
t(:,1) = 400;
t(1,:) = 600;
t(:,ny) = 800;
t(nx,:) = 900;
% defining the temperature at corner points
t(1,1) = 500;
t(nx,ny) = 850;
t(1,ny) = 700;
t(nx,1) = 650;
dt = 0.01; %time step
told=t;
alpha = 0.3; % Thermal diffusivity
tprev_dt = told;
tol = 1e-4;
k1 = alpha*dt/(dx^2);
k2 = alpha*dt/(dy^2);
% iterative method selection
method = input('Iterative method = ');
% JACOBI ITERATIVE METHOD
if (method==1)
tic;
jacobi = 1;
for p = 1:nt
error =100;
while(error>tol)
for i = 2:(nx-1)
for j = 2:(ny-1)
term1 = (1+(2*k1)+(2*k2))^-1;
term2 = k1*term1;
term3 = k2*term1;
h = (told(i-1,j)+told(i+1,j));
v = (told(i,j+1)+told(i,j-1));
t(i,j) = ((tprev_dt(i,j)*term1)+(h*term2)+(v*term3));
end
end
error = max(max(abs(t-told)));
told = t;
jacobi = jacobi+1;
end
tprev_dt = t ;
end
time_taken = toc;
% Plotting Code
figure(1)
contourf(x,y,t,'ShowText','on');
colormap(jet);
colorbar;
xlabel('X')
ylabel('Y')
text1 = sprintf('number of iteration = %d', jacobi);
text2 = sprintf('simulation time = %fn', time_taken);
title({'Transient-state Conduction Analysis Implicit - Jacobi'; text1; text2});
end
% GAUSS-SEIDEL ITERATIVE METHOD
if (method==2)
tic;
gauss = 1;
for p = 1:nt
error =100;
while(error>tol)
for i = 2:(nx-1)
for j = 2:(ny-1)
term1 = (1+(2*k1)+(2*k2))^-1;
term2 = k1*term1;
term3 = k2*term1;
h = (t(i-1,j)+told(i+1,j));
v = (told(i,j+1)+t(i,j-1));
t(i,j) = ((tprev_dt(i,j)*term1)+(h*term2)+(v*term3));
end
end
error = max(max(abs(t-told)));
told = t;
gauss = gauss+1;
end
tprev_dt = t;
end
time_taken = toc;
% Plotting Code
figure(1)
contourf(x,y,t,'ShowText','on');
colormap(jet);
colorbar;
xlabel('X')
ylabel('Y')
text1 = sprintf('number of iteration = %d', gauss);
text2 = sprintf('simulation time = %fn', time_taken);
title({'Transient-state Conduction Analysis Implicit - Gauss-seidel'; text1; text2});
end
% SUCCESSIVE OVER RELAXATION ITERATIVE METHOD
if (method==3)
tic;
sor = 1;
omega = 1.2; % Over relaxation factor
for p = 1:nt
error =100;
while(error>tol)
for i = 2:(nx-1)
for j = 2:(ny-1)
term1 = (1+(2*k1)+(2*k2))^-1;
term2 = k1*term1;
term3 = k2*term1;
h = (t(i-1,j)+told(i+1,j));
v = (told(i,j+1)+t(i,j-1));
t(i,j) = (told(i,j)*(1-omega))+(omega*((tprev_dt(i,j)*term1)+(h*term2)+(v*term3)));
end
end
error = max(max(abs(t-told)));
told = t;
sor = sor+1;
end
tprev_dt = t ;
end
time_taken = toc;
% Plotting Code
figure(1)
contourf(x,y,t,'ShowText','on');
colormap(jet);
colorbar;
xlabel('X')
ylabel('Y')
text1 = sprintf('number of iteration = %d', sor);
text2 = sprintf('simulation time = %fn', time_taken);
title({'Transient-state Conduction Analysis Implicit - SOR'; text1; text2});
end
Output Graphs -
Observation -
From the above simulation it is observed that the output is as silimar as what was in the steady state, and it was observed that SOR methods gives better results with no. of simulation and time elapsing on the run case.
----------------------------------------------------------------------------------------------------------------------------------------------
Objective - Transient state solution (EXPLICIT)
Theory -
Solving the above equation Explicitly :
Explicit Approach is used in numerical analysis for obtaining numerical approximations to the solution of time-dependent ordinary or partial differential equations as in required in computer simulations of a physical process. Explicit methods calculate the state of a system at a later time from the state of the system at the current time.
In this method, each differential equation approach contains only one unknown and can be solved explicitly for this one unknown using this method.
Program Codes -
% Solving Steady and Transient 2D heat conduction problem
% Transient state solution (EXPLICIT)
clear
close
clc
% initial conditions
nx = 10;
ny = 10;
nt = 300; % No.of steps
x = linspace(0,1,nx);
y = linspace(0,1,ny);
dx = x(2) - x(1);
dy = y(2) - y(1);
% initial temp condition
t = 300*ones(nx,ny);
t(:,1) = 400;
t(1,:) = 600;
t(:,ny) = 800;
t(nx,:) = 900;
% defining the temperature at corner points
t(1,1) = 500;
t(nx,ny) = 850;
t(1,ny) = 700;
t(nx,1) = 650;
dt = 0.01; %time step
told=t;
alpha = 0.3; % Thermal diffusivity
k1 = (alpha*dt)/(dx^2);
k2 = (alpha*dt)/(dy^2);
CFL = alpha*dt/(dx^2);
% EXPLICIT METHOD
% running the loop
tic;
explicit = 1;
for p = 1:nt
for i = 2:(nx-1)
for j = 2:(ny-1)
term1 = (told(i-1,j)-2*told(i,j)+told(i+1,j));
term2 = (told(i,j-1)-2*told(i,j)+told(i,j+1));
t(i,j) = ((told(i,j))+(k1*term1)+(k2*term2));
end
end
told = t;
explicit = explicit+1;
end
time_taken = toc;
% Plotting Code
figure(1)
contourf(x,y,t,'ShowText','on');
colormap(jet);
colorbar;
xlabel('X')
ylabel('Y')
text1 = sprintf('time N = %d', dt*nt);
text2 = sprintf('simulation time = %fn', time_taken);
title({'Transient-state Conduction Analysis Explicit'; text1; text2});
Output Graphs -
Observation -
Explicit method is converging faster than implicit methods.
Conclusion - After simulating all methods of convergence, It is concluded that explicit method is converging faster than implicit methods which is taking less time for simulations.