Week 5.1 - Mid term project - Solving the steady and unsteady 2D heat conduction problem
Considering General Three dimensional Heat Conduction Equation which has both time derivative and the spatial derivative terms.
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:
The above equation can be calculated by considering two approaches:
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.
Steady-state 2D heat conduction equation:
∂2T∂x2+∂2T∂y2=0∂2T∂x2+∂2T∂y2=0
2.Transient State equation:
Under these conditions, both terms are considered therefore the equation is as follows:
∂T∂t+α(∂2T∂x2+∂2T∂y2)=0∂T∂t+α(∂2T∂x2+∂2T∂y2)=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.
The assumption for solving the problem:
- Assume that the domain is a unit square.
- Assume nx = ny=10(for steady ) and 41(for unsteady) [Number of points along the x-direction is equal to the number of points along the y direction]
- Boundary conditions for steady and transient case
- Left:400K
- Top:600K
- Right:800K
- Bottom:900K
4.Tolerance=1e-4
5.Starting error=9e9
6.Thermal diffusivity(Alpha)=1.88e-5 m^2/sec
7.Doamin length=1m
8.CFL(for stability)=alpha*(dt/(dx^2))
Solving the Steady-state equation:
Considering Two dimensional Steady State Equation
here discretising the equation by applying the central differencing scheme
Code For Steady State Heat Transfer:-
%different solvers like jacobi,gauss seidal,Sucessive over relaxation
clear all
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,x,y,nx,ny,dx,dy,k,T,Told);
time_gauss_in_sec=toc
elseif iterative_solver==3
tic;
T_sor= sor_iter(start_error,error_tol,x,y,nx,ny,dx,dy,k,T,Told);
time_sor_in_sec=toc
end
%Function code for Jacobi:
%function for jacobi method to solve to steady state eqauation
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 code for Gauss-Seidal:
%function for gauss_seidal method to solve steady state method
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 @steady state',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 code for SOR:
%function for gauss_seidal method to solve steady state method
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);
end
Code for implicit Heat Transfer:-
%solving 2D Transient heat conduction equation im implicitly
clear all
close all
clc
%defining the inputs
nx=41;
ny=nx;
L=1;
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;%thermal diffusivity
[xx,yy] = meshgrid(x,y);
%defining the boundary conditions
T=298*ones(nx,ny);
T(1,:)=900;
T(end,:)=600;
T(:,1)=400;
T(:,end)=800;
%rounding off the corner value by taking the average value
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,start error
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));
%defining the solvers
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
%Code for jacobi function:
%function for jacobi method to solve to transient state eqauation(explicit)
function T=jacobi_iter_trans_impli(error_tol,xx,yy,nx,ny,T,k1,k2,Told,Tprev)
jacobi_iter=1;
for k= 1:120000 %time loop
start_error=9e9;
while(start_error > error_tol)%convergence loop
for i=2:nx-1 %nodal loop
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 implicitly using jacobi method with no of iteration=%d and time to solve=%f sec',jacobi_iter,time_jacobi);
title(title_text)
fprintf('No of Iteration in Jacobi =%d and time in sec=%f',jacobi_iter,time_jacobi);
end
%Code for Gauss Seidal function:
%function for GS method to solve to transient state eqauation(explicit)
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 implicitly using GS method with no of iteration=%d and time to solve=%f sec',gs_iter,time_gs));
title(title_text)
fprintf('No of Iteration in GS=%d and and time in sec=%f',gs_iter,time_gs);
end
%Code for SOR function:
%function for SOR method to solve to transient state eqauation(explicit)
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 implicitly using SOR method with no of iteration=%d and time to solve=%f sec',sor_iter,time_sor));
title(title_text)
fprintf('No of Iteration in SOR=%d and time in sec=%f',sor_iter,time_sor);
end
Code for explicit heat transfer:-
%solving 2D Transient heat conduction equation explicitly
clear all
close all
clc
%defining the inputs
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;%thermal diffusivity
[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;
%rounding off the corner value by taking the average value
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,start error
error_tol=1e-4;
start_error=9e9;
k1= alpha*(dt/dx^2);
k2= alpha*(dt/dy^2);
%defining the solvers
iterative_solver = 1;
tic;
if iterative_solver==1
T_jacobi= jacobi_iter_trasient_expli(xx,yy,nx,ny,T,k1,k2,Told);
end
time_jacobi_in_sec=toc
%Function code for Jacobi:
%function for jacobi method to solve to transient state eqauation(explicit)
function T= jacobi_iter_trasient_expli(xx,yy,nx,ny,T,k1,k2,Told)
jacobi_iter=1;
for k=1:120000
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\n',jacobi_iter);
end
Conclusion:
1.In Steady state equation SOR methods yields the reults very quickly when compare to other two methods.whereas Jacobi takes more number of iterations and more time to solve the unknows.
So the SOR method is most efficient when compared to the other two methods to solve the steady-state equation.
- For transient equations, the time required to solve is more when compared to steady-state equations.So, therefore it is recommended to solve the 2d heat conduction equation with steady state conditions.