Week 5.1 - Mid term project - Solving the steady and unsteady 2D heat conduction problem
AIM:- Mid term project - Solving the steady and unsteady 2D heat conduction problem
OBJECTIVE:-
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.
1. Jacobi
2. Gauss-seidel
3. Successive over-relaxation
Your absolute error criteria are 1e-4
Make sure to use the "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.
THEORY:-
Numerical Solution schemes are different in being explicit or implicit.
When a direct computation of the dependent variable can be made in terms of known quantities the computation is said to be explicit.
When the dependent variable is defined by sets of equations and either a matrix or iterative.
The technique is n added to obtain the solution to the numerical method is said to be implicit.
EXPLICIT METHOD:-
When the dependent variables are defined by coupled sets o equations, and either a matrix or iterative technique is needed to obtain the solution the numerical method is said to be explicit.
In an explicit numerical method, s would be evaluated in terms of known quantities at the previous time step n.
As a general rule, it can be shown that the condition cd<1 is very nearly equivalent to the stability condition for an explicit approximation.
Another general rule is that the time step sizes for explicit stability and accuracy are usually equivalent
it is this increase damping with the increase in time step size which produces inaccuracies in transient behavior
IMPLICIT METHOD:-
when the dependent variables are defined by coupled sets of equation sets of equations and either a matrix or iterative technique is needed to obtain the solution the numerical method is said to be implicit.
in Computational fluid dynamics (CFD) the governing equations are nonlinear and the number of unknown variables is typically very large.
iteration is used to advance a solution through a sequence of steps from a starting state to a final converged state.
The principal reason for using implicit solution methods which are more complex to program and require more computational either in each solution step is to allow for large time step sizes.
A simple qualitative model will help to illustrate how this works.
In order to solve for dependent variables iterative solvers used are
1. jacobian iterative solver
2. Gauss seidel iterative solver
3. Successive over-relaxation iterative solver.
* jacobian iterative solver
in numerical linear algebra, the Jacobi method is an iterative algorithm for determining the solutions of a strictly diagonally dominant system of numerical linear algebra.
Each diagonal element is solved for and an approximate value is plugged in
This process is then iterative until it converges
The algorithm is a stripped-down version of the Jacobi transformation method of matrix diagonalization.
*Gauss seidel iterative method:- is a technique used to solve a linear equations
The method is named other the German mathematics curl Friedrid gauss and Philip laud was for seidel
The method is similar to the Jacobi method and in the same, a strict or irreducible diagonal dominance of the system is sufficient to ensure coverage meaning the method will work.
* Successive over-relaxation iterative solver.
In this type of iterative solver, the relaxation factor is used in increasing the coverage rate. it goes by the formula
where`ω`>(`∂t/`∂x^2)+(∂t/`∂y^2) is the relaxation factor.
Solving 2D heat conduction equation for steady-state:-
Governing equation for the 2D heat conduction equation for steady-state is (∂t/`∂x^2)+(∂t/`∂y^2) =0
Now in order to solve the above governing equation, we first discretize the equation by using a central differencing scheme are as follows:-
Equation 1 will be solved by using the point iterative method.
The point iterative method which will be used to solve the above equation 1 is the above iteration method which are Jacobi, gauss Siedel, and the SOR method
MATLAB CODE:-
. In the following program, we have to solve equation 1 by using different iterative methods.
Matlab code for steady-state 2D heat conduction profile :
clear all
close all
clc
%specifying length of domain and no of grid points
l=1;
nx=10;
ny=nx;
% Defining space domain
x = linspace(0,1,nx);
dx = x(2)-x(1);
y = linspace(0,1,ny);
dy = y(2)-y(1);
%creating tempeatures across given region
t = 300*ones(nx,ny);
t(:,1) = 400;
t(1,:) = 600;
t(:,end) = 800;
t(end,:) = 900;
% Tempeature at edges
t(1,1) = (600+400)/2;
t(1,end)=(600+800)/2;
t(end,1)= (400+900)/2;
t(end,end) = (800+900)/2;
% updating old values
told = t;
% solving 2D steady state system
[xx,yy] = meshgrid(x,y);
% Assigning tolerence and error values
tolerance=1e-4;
error = 1e8;
% Defining k factor
k = 2*((1/(dx^2))+(1/(dy^2)));
%iteration method input
iter_method = input('Enter the iteration method number= ')
% jacobi method
%Defining iteration varible to store and display no of iterations computated
if iter_method == 1
tic
jacobi_iter =1;
while(error>tolerance)
for i = 2:nx-1
for j = 2:ny-1
term1 = (told(i-1,j)+ told(i+1,j))/(k*(dx^2));
term2 = (told(i,j-1)+told(i,j+1))/(k*(dy^2));
t(i,j) = term1+term2;
end
end
error = max(max(abs(told-t)));
told = t;
jacobi_iter = jacobi_iter+1;
end
time1=toc;
end
%Gauss siedal method
if iter_method==2
tic
gs_iter=1;
while(error>tolerance)
for i=2:nx-1
for j=2:ny-1
% The brackets we have inserted was the problem. The
% equation would be incorrect if we put the extra bracket
% at the starting. I am sending the difference in equations because of the brackets.
% term1 = ((t(i-1,j)+ told(i+1,j))/k*(dx^2));
%term2 = ((t(i,j-1)+ told(i,j+1))/k*(dy^2));
term1 = (t(i-1,j)+ told(i+1,j))/(k*(dx^2));
term2 = (t(i,j-1)+ told(i,j+1))/(k*(dy^2));
t(i,j) = term1+term2;
end
end
error = max(max(abs(told - t)));
told = t;
gs_iter=gs_iter+1;
end
time2 =toc;
end
% SOR Method
if iter_method ==3
omega =1.8;
tic
sor_iter=1;
while(error>tolerance)
for i = 2:nx-1
for j = 2:ny-1
term1 = (t(i-1,j)+told(i+1,j))/(k*(dx^2));
term2 =(t(i,j-1)+told(i,j+1))/(k*(dy^2));
t(i,j) = (1-omega)*told(i,j)+omega*(term1+term2);
end
end
error = max(max(abs(told-t)));
told=t;
sor_iter=sor_iter+1;
end
time3=toc;
end
%plotting the results tempeature
[c,h] = contourf(xx,yy,t);
clabel(c,h);
colorbar;
set(gca,'YDir','reveres')
colormap(jet);
if iter_method ==1
title({'study state tempeature profile';['iteration Method = ',num2str(iter_method)];['No iterations=',num2str(jacobi_iter)];['Simulation time= ',num2str(time1),'S']})
elseif iter_method ==2
title({'Steady state tempeature profile';['iteration Method=',num2str(iter_method)];['No of iterations = ',num2str(gs_iter)];['Simulation time= ',num2str(time2),'S']})
elseif iter_method ==3
title({'Steady state tempeature profile';['iterations Method= ',num2str(iter_method)];['No of iterations = ',num2str(sor_iter)];['simulation time = ',num2str(time3),'S']})
end
xlabel('x_length Domain')
ylabel('y_length Domain')
RESULT:-
Plot for 2D heat conduction for steady-state using Jacobi iterative method:-
Steady-state temperature profile iteration Method = 1
No of iterations = 15940
simulation time = 3.4096s
plot for 2D heat conduction for steady-state using gauss Siedel iterative method:
Steady-state temperature profile
iteration Method =2
No of iterations = 2
Simulation time = 1.9161s
plot for 2D heat conduction for steady-state using SOR iterative method:-
Steady-state temperature profile
iteration Method = 3
Simulation time = 0.34395s
Explicit schemes:-
After applying the first-order forward differencing scheme on L.H.S. and central differencing scheme on R.H.S. on equation 2 we get
Matlab code for unsteady-state 2D heat conduction profile in explicit method:
clear all
close all
clc
l = 1
nx = 60;
ny = 60;
%Defining space domain
x = linspace(0,1,nx);
dx = x(2)-x(1);
y = linspace(0,1,ny);
dy = y(2)-y(1);
%calculation of time step
%dt = (0.5)/((alpha/dx^2)+(alpha/dy^2));
dt = 0.00001;
alpha = 1.4;
% creating tempeatures across given region
t = 300*ones(nx,ny);
t(:,1)= 400;
t(1,:) = 600;
t(:,end) = 800;
t(end,:) = 900;
% Tempeatures at edges
%t(1,1)= (600+400)/2;
%t(1,end) = (600+800)/2;
%t(end,1) = (400+900)/2;
%t(end,end)= (800+900)/2;
%updating at t values
told = t;
CFL_num =(alpha*dt)/((1/dx^2)+(1/dy^2));
% solving 2D unsteady state system
%[xx,yy] = meshgrid(x,y);
k1 = (alpha*dt)/(dx^2);
k2 = (alpha*dt)/(dy^2);
iter_method=1;
tic;
for k = 1:7000
for i = 2:nx-1
for j = 2:ny-1
term1=(1-2*k1-2*k2);
term2 = k1;
term3 = k2;
H = (told(i-1,j)+told(i+1,j));
V = (told(i,j-1)+told(i,j+1));
t(i,j) = (told(i,j)*term1)+(term2*H)+(term3*V);
end
end
told = t;
end
time = toc;
figure(1)
colormap(jet)
contourf(x,y,t,'showtext','on');
colorbar
%clabel(c,h);
title({'unsteady tempeature profile';['simulation time = ',num2str(time),'s']})
xlabel('xlength Domain');
ylabel('ylength Domain');
RESULTS:-
implicit scheme:-
governing equation for 2D heat conduction for the unsteady state is 
∂T/∂t=(∂t/`∂x^2)+(∂t/`∂y^2) =0..........(2)
After applying the first-order forward differencing scheme on L.H.S. and central differencing scheme on R.H.S. on equation 2D, we get.
Matlab code for unsteady-state 2D heat conduction profile for implicit method:-
clear all
close all
clc
nx = 10;
ny = 10;
x = linspace(0,1,nx);
dx = 1/(nx-1);
y = linspace(0,1,ny);
dy = 1/(ny -1);
alpha = 1.4;
%dt =(0.5)/((alpha/dx^2)+(alpha/dx^2));
dt = 0.001;
tolerance = 1e-4;
%error = 1e9;
t = 300*ones(nx,ny);
t(:,1) = 400;
t(1,:) = 600;
t(:,end) = 800;
t(end,:) = 900;
%Tempeature at edges
t(1,1) = (600+400)/2;
t(1,end) = (400+800)/2;
t(end,1) = (400+900)/2;
t(end,end) = (800+900)/2;
%upadting old t values
told = t;
tprev =told;
%solving 2D steady state system
k1 = (alpha*dt)/(dx^2);
k2 = (alpha*dt)/(dy^2);
cFL_num = k1+k2;
nt = 300; %number of time step
iter_method =input('Enter the iteration method number=');
if iter_method == 1
%computing Time loop
tic;
jacobi_iter = 1;
for k =1:nt
%jacobian iterative solve
error =1e3;
while(error>tolerance)
for i = 2:nx-1
for j=2:ny-1
term1 = 1/(1+2*k1+2*k2);
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(i,j)*term1)+(term2*H)+(term3*V);
end
end
% error = mx(max(abs(told - t)));
error = max(max(abs(t - told)));
told = t;
jacobi_iter = jacobi_iter+1;
end
tprev = t;
end
time1 = toc;
% figure(1)
% [c,h] = contourf(x,y,t);
% set(gca, 'YDir','reverse')
% clabel(c,h);
% colormap(jet);
end
%Gauss siedal iteration method
if iter_method == 2
gs_iter =1;
tic
for nt = 1:600
term1 = 1/(1+2*k1+ 2*k2);
term2 = k1*term1;
term3 = k2*term1;
error = 1e9;
while(error>tolerance)
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) = (told(i,j)*term1)+(term2*H)+(term3*V);
end
end
error = max(max(abs(told-t)));
told = t;
gs_iter=gs_iter+1;
end
time2 =toc;
end
end
%SOR iteration method
omega = 1.2;
if iter_method ==3
sor_iter= 1;
tic
for nt = 1:600
term1 = 1/(1+2*k1+2*k2);
term2 = k1*term1;
term3 = k2*term1;
error = 1e9;
while(error>tolerance)
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*told(i,j)+term2*H+term3*V);
end
end
error = max(max(abs(told - t)));
told = t;
sor_iter = sor_iter+1;
end
time3 = toc;
end
end
figure(1)
[c,h] = contourf(x,y,t);
clabel(c,h);
colorbar;
colormap(jet);
if iter_method ==1
title({'unsteady state tempeature profile';['iteration Method =jacobi_iter'];['No of iterations= ',num2str('jacobi_iter')]})
elseif iter_method ==2
title({'unsteady state tempeature profile';['iterations method = gauss siedal']; ['No of iterations = ',num2str(gs_iter)]})
elseif iter_method ==3
title({'unsteady state tempeature profile';['iterations Method = SOR'];['No of iterations = ',num2str(sor_iter)]})
end
Result-
plot for unsteady-state 2D heat conduction in the implicit method for jacobian iteration:-
plot for unsyeady state 2D heat conduction in implicit method for gauss siedel iteratons:
plot for unsteady state 2D heat conduction in implicit method for SOR iterations:-
Conclusion and Explanation:-
Though we are getting the same result for different iterations and computation time taken by the iterative solver to get convergence result.
By observing the above plots it is clear that Jacobi iterations such as gauss Siedel and SOR.
The computation time is also less in jacobian iterations than in the other two iterative methods.
The implicit method requires guess values iteration for convergence as it also based on the the future state values
Which are not known at the time f computations
it checks convergence for each time step.
This caused a delay in the simulation.
The implicit method can be used for problems when converged results are needed in less computation time as they are unconditionally stable.
with an increase in time size, the convergence rate is increased which leads to a steady state with less no of iterations.
The explicit method results in a time-accurate solution as it is governed by stability in criteria in the solution and when is unstable it requires more no iterations