Week 5.1 - Mid term project - Solving the steady and unsteady 2D heat conduction problem
Aim: To write the code which solves the 2D heat conduction equation by using the point iterative techniques by implementing the following methods.
1. jacobi
2. Gauss Siedel
3. Successive over-Relaxation
Introduction:
Considering General Three dimensional Heat conduction equation which has both time derivative and the spatial derivative terms
`(delT)/(delt)+alpha((del^2T)/(delx^2)+(del^2T)/(dely^2)+(del^2T)/(delz^2))=0`
The above equation has 2 terms where term1 represents time derivative and term2 is spatial derivative where alpha is thermal diffusivity.
Now by considering the assumption that the above equation is solved for a thin plate, dimensions in the Z direction are neglected, thereby the above equation is converted into the following form which is as follows:
`delT/del(t)+alpha((del^2T)/(delx^2)+(del^2T)/(dely^2))=0`
The above equation can be calculated by considering two approaches:
1. Steady State Equation:
Under these conditons, the rate of change of temperature remains 0. Therefore in this condition time, the derivative term is equated to zero.
Steady state 2D heat conduction equation:
`((del^2T)/(delx^2)+(del^2T)/(dely^2))=0`
2. Transient State equation:
Under these conditions, both terms are considered therefore the equation is as follows:
`delT/del(t)+alpha((del^2T)/(delx^2)+(del^2T)/(dely^2))=0`
Above two equations can be solved in 2 ways by using 3 iterative methods such as Jacobi, Gauss-Seidel and SOR.
1. Explicit Method
2. Implicit Method
Here in this Project Transient State equation is solved by using both ways.
Assumptions:
1. Assume that the domain is a unit square.
2. Assume nx=ny=10 (For Steady) and 41 (For Unsteady)
3. Boundary condtions for Steady and Transient Case
1. Left=400k
2.Right=800k
3.Top=600k
4.Bottom=900k
4. Tolerance=1e-4
5. Starting error=9e9
6. Thermal Diffusivity (Alpha)=1.88e-5m^2/sec
7.Domain Length = 1m
8. CFL(For Stability)=alpha8(dt/(dx^2))
Solving the Steady State Equation
Considering Two dimensional steady state eqn
`((del^2T)/(delx^2)+(del^2T)/(dely^2))=0`
here discretising the equation by applying central differencing scheme
`T(i,j)=1/k{(Told(i-1,j)+Told(i+1,j))/dx^2)+((Told(i,j-1)+Told(i,j+1))/dy^2)}`
%Solving the steady state equation by employing
%different solvers like Jacobi, Gauss Seidel, SOR
close all
clc
%defining the inputs
nx=10;
ny=nx;
x=linspace(0,1,nx);
y=linspace(0,1,nx);
dx=x(2)-x(1);
dy=y(2)-y(1);
[xx,yy]=meshgrid(x,y);
%defining the boundary conditions
T=ones(nx,ny);
T(1,:)=900;
T(end,:)=600;
T(:,1)=400;
T(:,end)=800;
Told=T;
%error_tol,start error
error_tol=1e-4;
start_error=9e9;
k=2*((dx^2+dy^2))/(dx^2*dy^2);
%defining the solvers
iterative_solver=input('Enter the desired number with respect to solver:');
if iterative_solver==1
tic;
T_jacobi=jacobi_iter(start_error,error_tol,xx,yy,nx,ny,dx,dy,T,k,Told);
time_jacobi_in_sec=toc;
elseif iterative_solver==2
tic;
T_gauss=gauss_iter(start_error,error_tol,xx,yy,nx,ny,dx,dy,T,k,Told);
time_gauss_in_sec=toc;
elseif iterative_solver==3
tic;
T_sor=sor_iter(start_error,error_tol,xx,yy,nx,ny,dx,dy,T,k,Told);
time_sor_in_sec=toc;
end
function T=jacobi_iter(start_error,error_tol,xx,yy,nx,ny,dx,dy,T,k,Told)
jacobi_iter=1;
while(start_error>error_tol)
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
start_error=max(max(abs(Told-T)));
Told=T;
jacobi_iter=jacobi_iter+1;
figure(1)
[c,d]=contourf(xx,yy,T);
colormap(jet)
colorbar
title_text=sprintf('Jacobi iteration number=%d @steady state',jacobi_iter);
title(title_text)
xlabel('X-axis')
ylabel('y-axis')
clabel(c,d)
pause(0.03)
end
fprintf('No of iteration in jacobi=%d',jacobi_iter);
end
function gauss_seidal=gauss_iter(start_error,error_tol,x,y,nx,ny,dx,dy,k,T,Told)
GS_iter=1;
while(start_error>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)+((T(i,j-1)+Told(i,j+1))/dy^2))/k;
end
end
start_error=max(max(abs(Told-T)));
Told=T;
GS_iter=GS_iter+1;
gauss_seidal=figure(2);
[c,d]=contourf(x,y,T);
colormap(jet)
colorbar
title_text=sprintf('GS Iteration number=%d',GS_iter);
title(title_text)
xlabel('X-axis')
ylabel('Y-axis')
clabel(c,d)
pause(0.03)
end
fprintf('No of iteration in GS=%d',GS_iter);
end
function T=sor_iter(start_error,error_tol,x,y,nx,ny,dx,dy,k,T,Told)
sor_iter=1;
omega=1.7;
while(start_error>error_tol)
for i=2:nx-1
for j=2:ny-1
T_gs=(((T(i-1,j)+Told(i+1,j))/dx^2)+((T(i,j-1)+Told(i,j+1))/dy^2))/k;
T(i,j)=omega*(T_gs)+(1-omega)*(Told(i,j));
end
end
start_error=max(max(abs(Told-T)));
Told=T;
sor_iter=sor_iter+1;
figure(3)
[c,d]=contourf(x,y,T);
colormap(jet)
colorbar
title_text=sprintf('SOR iteration number=%d @steady state',sor_iter);
title(title_text)
xlabel('X axis')
ylabel('Y-axis')
clabel(c,d)
pause(0.03)
end
fprintf('No of iteration in SOR=%d',sor_iter);
endFrom all the above plots it can be infereed that the SOR method takes less time and less number of iterations when compared to the other two methods. Therefore SOR stands out to be good and fastest among all.
Solving Transient State Equation:
Transiet state 2D steady state heat conduction equation is solved in both ways i.e. explicity and implicity. Following is the 2D Heat conduction equation.
`delT/del(t)+alpha((del^2T)/(delx^2)+(del^2T)/(dely^2))=0`
Solving the above equation Explicity
Explicit approach is used in muerical analysis for obtaining numerical approximations to the solution of time-dependent oridnary or partial differential equations as in required in computer simulations of a physical process. Explicit methods caluculate 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 explicity for this one unknown using this method.
Here in this equation time derivative terms are calculated using the first order differential equations and whereas the term 2 has second-order derivative terms and elliptical in nature so information needs to be taken from all points equally, So therefore a central differencing scheme is applied for the discretising the term 2.
`T^(n+1)(i,j)=T^n+k1*(T^nold(i+1,j)-2*T^nold(i,j)+T^nold(i-1,j))+k2*[T^nold(i,j-1)-2*T^nold(i,j)+T^nold(i,j+1)]`
`k1=2*((alpha*dt)/dx^2) , k2=2*((alpha*dt)/dy^2)`
close all
clc
nx=41;
ny=nx;
x=linspace(0,1,nx);
y=linspace(0,1,nx);
dx=x(2)-x(1);
dy=y(2)-y(1);
dt=0.01;
alpha=1.88e-5;
[xx,yy]=meshgrid(x,y);
T=ones(nx,ny);
T(1,:)=900;
T(end,:)=600;
T(:,1)=400;
T(:,end)=800;
Told=T;
T(1,1)=(400+900)/2;
T(1,ny)=(900+800)/2;
T(nx,1)=(400+600)/2;
T(nx,ny)=(600+800)/2;
error_tol=1e-4;
start_error=9e9;
k1=alpha*(dt/dx^2);
k2=alpha*(dt/dy^2);
iterative_solver=1;
tic;
if iterative_solver==1
T_jacobi=jacobi_iter_transient(xx,yy,nx,ny,T,k1,k2,Told);
end
time_jacobi_in_sec=toc;
function T=jacobi_iter_transient(xx,yy,nx,ny,T,k1,k2,Told)
jacobi_iter=1;
for k=1:12000
for i=2:nx-1
for j=2:ny-1
term1=Told(i,j);
term2=k1*(Told(i+1,j)-2*Told(i,j)+Told(i-1,j));
term3=k2*(Told(i,j+1)-2*Told(i,j)+Told(i,j-1));
T(i,j)=term1+term2+term3;
end
end
start_error=max(max(abs(Told-T)));
Told=T;
jacobi_iter=jacobi_iter+1;
end
figure(1)
[c,d]=contourf(xx,yy,T);
colormap(jet)
colorbar
title_text=sprintf('Jacobi iteration number=%d @ transient state explicitly',jacobi_iter);
title(title_text)
xlabel('X-Axis')
ylabel('Y-Axis')
clabel(c,d)
pause(0.03)
fprintf('No of iteration in jacobi = %d',jacobi_iter);
end
No of iteration in jacobi=120001
time_jacobi_in_sec=21.6943
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 Y9t) 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))
Here in this equation time derivative terms are calculated using the first order differential equations and whereas the term 2 has second-oder derivative terms and elliptical in nature so information needs to be taken from all the points equally, So therefore a central differencing scheme is applied for the discretising the term2.
`T^(n+1)(i,j)=T^n+(k1*(T^nold(i+1,j)-2*T^nold(i,j)+T^nold(i-1,j)))^(n+1)+(k2*[T^nold(i,j-1)-2*T^nold(i,j)+T^nold(i,j+1)])^(n+1)`
`k1=2*((alpha*dt)/dx^2) , k2=2*((alpha*dt)/dy^2)`
In 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 difference between T at 'n+1' which is obtained and T at 'n+1' is assumed as the value of T which was obtained in the previous iteration.
Once the value of T is found out iteratively whose error is less than the tolerance value, the value of T is calculated at the next time step 'n+2' using the same iterative method.
Equation for Jacobi method are:
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(i,j)*term1)+term2*H+term3*v
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.
term1=(1+(2*k1)+(2*k2))^-1;
term2=k1*term1;
term3=k2*term1;
H=(T(i-1,j)+Told(i+1,j));
v=(T(i,j-1)+Told(i,j+1));
T(i,j)=(Tprev(i,j)*term1)+term2*H+term3*v
SOR: This method can be applied either by using Jacobi or Gauss seidal but the only difference is consideration of factor called omeda(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 fomula
omega=2/(1+sin(pi*dx))
Equations for SOR:
term1=(1+(2*k1)+(2*k2))^-1;
term2=k1*term1;
term3=k2*term1;
H=(T(i-1,j)+Told(i+1,j));
V=(T(i,j-1)+Told(i,j+1));
T(i,j)=omega*(Tprev(i,j)*term1+term2*H+term3*V)+(1-omega)*Tprev(i,j)
close all
clc
nx=41;
ny=nx;
x=linspace(0,1,nx);
y=linspace(0,1,nx);
dx=x(2)-x(1);
dy=y(2)-y(1);
dt=0.01;
alpha=1.88e-5;
[xx,yy]=meshgrid(x,y);
T=298*ones(nx,ny);
T(1,:)=900;
T(end,:)=600;
T(:,1)=400;
T(:,end)=800;
T(1,1)=(400+900)/2;
T(1,ny)=(900+800)/2;
T(nx,1)=(400+600)/2;
T(nx,ny)=(600+800)/2;
Told=T;
Tprev=Told;
error_tol=1e-4;
start_error=9e9;
k1=alpha*(dt/dx^2);
k2=alpha*(dt/dy^2);
omega=2/(1+sin(pi*dx));
CFL=alpha*(dt/(dx^2));
iterative_solver=input('Input the appropriate number to select required the method');
if iterative_solver==1
tic;
T=jacobi_iter_trans_impli(error_tol,xx,yy,nx,ny,T,k1,k2,Told,Tprev);
elseif iterative_solver==2
tic;
T_gs=gs_iter_trans_impli(error_tol,xx,yy,nx,ny,T,k1,k2,Told,Tprev);
elseif iterative_solver==3
tic;
T_sor=sor_iter_trans_impli(error_tol,xx,yy,nx,ny,T,k1,k2,Told,Tprev,omega);
end
function T=jacobi_iter_trans_impli(error_tol,xx,yy,nx,ny,T,k1,k2,Told,Tprev)
jacobi_iter=1;
for k=1:120000
start_error=9e9;
while(start_error>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(i,j)*term1)+term2*H+term3*v;
end
end
start_error=max(max(abs(Told-T)));
Told=T;
jacobi_iter=jacobi_iter+1;
end
Tprev=T;
end
time_jacobi=toc;
contourf(xx,yy,T,'ShowText','on')
colormap(jet)
colorbar
xlabel('X-Axis')
ylabel('Y-axis')
colorbar
title_text=sprintf('Solution for 2D dimensional heat conduction equation solved implicity');
title(title_text)
fprintf('No of iteration = %d and time in sec = %f',jacobi_iter,time_jacobi);
end
function T=gs_iter_trans_impli(error_tol,xx,yy,nx,ny,T,k1,k2,Told,Tprev)
gs_iter=1;
for k=1:120000
start_error=9e9;
while(start_error>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=(T(i,j-1)+Told(i,j+1));
T(i,j)=(Tprev(i,j)*term1)+term2*H+term3*v;
end
end
start_error=max(max(abs(Told-T)));
Told=T;
gs_iter=gs_iter+1;
end
Tprev=T;
end
time_gs=toc;
contourf(xx,yy,T,'ShowText','on')
colormap(jet)
colorbar
xlabel('X-Axis')
ylabel('Y-axis')
colorbar
title_text=sprintf('Solution for 2D dimensional heat conduction equation solved implicity');
title(title_text)
fprintf('No of iteration = %d and time in sec = %f',gs_iter,time_gs);
end
function T=sor_iter_trans_impli(error_tol,xx,yy,nx,ny,T,k1,k2,Told,Tprev,omega)
sor_iter=1;
for k=1:120000
start_error=9e9;
while(start_error>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=(T(i,j-1)+Told(i,j+1));
T(i,j)=omega*(Tprev(i,j)*term1+term2*H+term3*v)+(1-omega)*Tprev(i,j);
end
end
start_error=max(max(abs(Told-T)));
Told=T;
sor_iter=sor_iter+1;
end
Tprev=T;
end
time_sor=toc;
contourf(xx,yy,T,'ShowText','on')
colormap(jet)
colorbar
xlabel('X-Axis')
ylabel('Y-axis')
colorbar
title_text=sprintf('Solution for 2D dimensional heat conduction equation solved implicity');
title(title_text)
fprintf('No of iteration = %d and time in sec = %f',sor_iter,time_sor);
end 
Jacobi
No of iteration = 240826 and time in sec =3.840477
Gauss Seidal
No of iteration = 240001 and time in sec =4.189748
SOR
No of iteration = 241645 and time in sec = 4.617075
Conclusion:
1. In steady state equation SOR methods yileds the results the results very quickly when compared to other two methods whereas jacobi more number of iterations and more time to solve the unknowns.
So the SOR method is most efficient when compared to the other two methods to solve the Steady-State equations.
2. For transient equations, the time required to solve is more when compared to steady-state equations. So, therefore it is recommended to solve the 2nd heat conduction equation with steady state conditions.
In transiet state equatio, explicit methods yields the results quickly when compared to the implict method, but it is recommended under the situation when CFL is less than 0.35.
When transient state equation is solved implicitly Jacobi takes less them when compared to the other two methods whereas SOR is more stable compared to gauss seidal.