2D Steady and Tranisent heat conduction problem
Objective:
The aim of this project to solve the 2D heat conduction steady and transient equation using different types of iterative solvers.
A domain length and height of `1xx1 m^2`was used for the probem. The equation for the temperature at each time step for the tranisent equation was estimated using the central difference scheme. The iterative methods used were Jacobian, Gauss-Siedel and SOR on Gauss method.
The general equation for the Heat conduction is given by:
`(delT)/(delt) =alpha[(del^2T)/(del^2x)+(del^2T)/(del^2y)]`
This is the equation for the transient state.
The equation for the steady state is given by:
`(del^2T)/(del^2x)+(del^2T)/(del^2y)=0`
The above equation is a second order derivative or an elliptical equation. So a central differnce scheme can be used to solve this.
`(T_(i+1,j)-2T_(i,j)+T_(i-1,j))/(Deltax^2) + (T_(i,j+1)-2T_(i,j)+T_(i,j-1))/(Deltay^2)=0`
After solving this equation for `T_(i,j)` we get the equation as
`T_(i,j)=1/k[(T_(i+1,j)+T_(i-1,j))/(Deltax^2) + (T_(i,j+1)+T_(i,j-1))/(Deltay^2)]`
where `k=[2/(Deltax^2) + 2/(Deltay^2)]`
The equation for the transient state is given by:
`(delT)/(delt)=alpha[(del^2T)/(del^2x)+(del^2T)/(del^2y)]`
We have the time coordinates as well the spatial coordinates in this case. So we would be using the time marching solution(forward differencing) for the time loops and central difference for the spatial coordinates
on applying the central difference scheme on this equation we get:
`[(T_(i,j)^(n+1)- T_(i,j)^(n))/(Deltat)]=alpha[[T_(i+1,j)-2T_(i,j)+T_(i-1,j))/(Deltax^2) + (T_(i,j+1)-2T_(i,j)+T_(i,j-1))/(Deltay^2)]^(n)`
On solving this equation we get:
`T_(i,j)^(n+1)=T_(i,j) + k_(1)xx(T_(i+1,j) +T_(i-1,j))+ k_(2)xx(T_(i,j+1) + T_(i,j-1))`
where `k_(1)=alpha((Deltat)/(Deltax^2))`
`k_(2)=alpha((Deltat)/(Deltay^2))`
The above equations are used for explicit schemes.
For implicit scheme
`(delT)/(delt)=alpha[(del^2T)/(del^2x)+(del^2T)/(del^2y)]`
We have the time coordinates as well the spatial coordinates in this case. So we would be using the time marching solution(forward differencing) for the time loops and central difference for the spatial coordinates
on applying the central difference scheme on this equation we get:
`[(T_(i,j)^(n+1)- T_(i,j)^(n))/(Deltat)]=alpha[[T_(i+1,j)-2T_(i,j)+T_(i-1,j))/(Deltax^2) + (T_(i,j+1)-2T_(i,j)+T_(i,j-1))/(Deltay^2)]^(n+1)`
On solving this equation we get:
`T_(i,j)^(n+1)=k_(3)xxT_(i,j) + k_(1)xxk_(3)xx(T_(i+1,j) +T_(i-1,j))+ k_(2)xxk_(3)xx(T_(i,j+1) + T_(i,j-1))`
where `k_(1)=alpha((Deltat)/(Deltax^2))`
`k_(2)=alpha((Deltat)/(Deltay^2))`
`k_(3)=1/(1+2k_(1) +2k_(2))`
In the iterative solvers for the Jacobian matrix we use the old values for T while in the Gauss Siedel method we use the updated value of T.
For SOR we use the relaxtion factor `omega` which is given by the equation
`omega=2/(1+sin(pixxdx))`
The Matlab Code for the Steady state equation is given below:
clear all
close all
clc
% inputs for iterative solving method for unsteady state equation
L=1;
B=1;
alpha=0.05;
nx=30;
ny=nx;
x=linspace(0,L,nx);
y=linspace(0,B,ny);
dx=x(2)-x(1);
dy=y(2)-y(1);
theta=(cos(pi/nx))+(cos(pi/ny));
omega=4/(2+sqrt(4-theta^2));
T1=400*ones(nx,ny);
T1(2:end-1,1)=400;
T1(1,2:end-1)=600;
T1(nx,2:end-1)=900;
T1(2:end-1,ny)=800;
T1(1,1)=500;
T1(1,end)=700;
T1(end,1)=650;
T1(end,end)=850;
T=T1;
Told=T;
iterative_solver=3;
error=9e9;
tol=1e-4;
k=2*(dx^2 +dy^2)/(dx^2*dy^2);
omega=2/(1+sin(pi*dx));
if iterative_solver==1
tic
T=Jacobiansteady(T,x,y,nx,ny,dx,dy,error,tol,k);
time_taken=toc;
else if iterative_solver==2
tic
T=Gausssteady(T,x,y,nx,ny,dx,dy,error,tol,k);
time_taken=toc;
else if iterative_solver==3
tic
T=SORsteady(T,x,y,nx,ny,dx,dy,error,tol,k,omega);
time_taken=toc;
end
end
end
The Jacobian Function:
function T=Jacobiansteady(T,x,y,nx,ny,dx,dy,error,tol,k)
[X,Y]=meshgrid(x,y);
Told=T;
jacobiiter=1;
while (error>tol)
for i=2:nx-1
for j=2:ny-1
T1(i,j)= (Told(i+1,j)+Told(i-1,j))/(k*dx^2);
T2(i,j)= (Told(i,j+1)+Told(i,j-1))/(k*dy^2);
T(i,j)=T1(i,j)+T2(i,j);
end
end
error=max(max(abs(Told-T)));
Told=T;
jacobiiter=jacobiiter+1;
figure(1)
C=contourf(X,Y,T);
set(gca,\'YDIR\',\'reverse\');
colormap(jet)
clabel(C,\'Fontsize\',14);
title(sprintf(\'No of iterations for steady state jacobi: %d\',jacobiiter));
end
end The Gauss Siedel function:
function T=Gausssteady(T,x,y,nx,ny,dx,dy,error,tol,k)
[X,Y]=meshgrid(x,y);
Told=T;
gaussiiter=1;
while (error>tol)
for i=2:nx-1
for j=2:ny-1
T1(i,j)= (T(i+1,j)+T(i-1,j))/(k*dx^2);
T2(i,j)= (T(i,j+1)+T(i,j-1))/(k*dy^2);
T(i,j)=T1(i,j)+T2(i,j);
end
end
error=max(max(abs(Told-T)));
Told=T;
gaussiiter=gaussiiter+1;
figure(1)
C=contourf(X,Y,T);
set(gca,\'YDIR\',\'reverse\');
colormap(jet)
clabel(C,\'Fontsize\',14);
title(sprintf(\'No of iterations for steady state GaussSiedel: %d\',gaussiiter));
end
end The SOR on Guass function:
function T=SORsteady(T,x,y,nx,ny,dx,dy,error,tol,k,omega)
[X,Y]=meshgrid(x,y);
Told=T;
sor=1;
while (error>tol)
for i=2:nx-1
for j=2:ny-1
T1(i,j)= (T(i+1,j)+T(i-1,j))/(k*dx^2);
T2(i,j)= (T(i,j+1)+T(i,j-1))/(k*dy^2);
T(i,j)=omega*(T1(i,j)+T2(i,j))-(omega-1)*(T(i,j));
end
end
error=max(max(abs(Told-T)));
Told=T;
sor=sor+1;
figure(1)
C=contourf(X,Y,T);
set(gca,\'YDIR\',\'reverse\');
colormap(jet)
clabel(C,\'Fontsize\',14);
title(sprintf(\'No of iterations for steady state SOR: %d\',sor));
end
end The Results obtained from the Code are presented in the figure below:
Jacobian Steady state Iterations
Gauss Siedel Method:
SOR on Gauss Siedel Method:
Matlab Code for the Transient Heat conduction equation:
clear all
close all
clc
% inputs for iterative solving method for unsteady state equation
L=1;
B=1;
alpha=0.05;
nx=30;
ny=nx;
x=linspace(0,L,nx);
y=linspace(0,B,ny);
dx=x(2)-x(1);
dy=y(2)-y(1);
dt=1e-3;
theta=(cos(pi/nx))+(cos(pi/ny));
omega=4/(2+sqrt(4-theta^2));
T1=400*ones(nx,ny);
T1(2:end-1,1)=400;
T1(1,2:end-1)=600;
T1(nx,2:end-1)=900;
T1(2:end-1,ny)=800;
T1(1,1)=500;
T1(1,end)=700;
T1(end,1)=650;
T1(end,end)=850;
T=T1;
Told=T;
T_prev_dt=Told;
iterative_solver=5;
error=9e9;
ts=0;
t=3;
nt=t/dt;
tol=1e-4;
k1=alpha*(dt/(dx^2));
k2=alpha*(dt/(dy^2));
k3=1/(1+(2*k1)+(2*k2));
omega=2/(1+sin(pi*dx));
if iterative_solver==1
tic
T=Jacobiansolver(T,x,y,nx,ny,nt,error,tol,k1,k2);
time_taken=toc;
else if iterative_solver==2
tic
T=GaussSeidelsolver(T,x,y,nx,ny,nt,error,tol,k1,k2);
time_taken=toc;
else if iterative_solver==3
tic
T=SORgausssolver(T,x,y,nx,ny,nt,error,tol,k1,k2);
time_taken=toc;
else if iterative_solver==4
tic
T=Jacobian_implicit_solver(T,Told,T_prev_dt,x,y,nx,ny,nt,error,tol,k1,k2,k3);
time_taken=toc;
else if iterative_solver==5
tic
T=GaussSiedel_implicit_solver(T,Told,T_prev_dt,x,y,nx,ny,nt,error,tol,k1,k2,k3);
time_taken=toc;
else if iterative_solver==6
tic
T=SORGaussSiedel_implicit_solver(T,Told,T_prev_dt,x,y,nx,ny,nt,error,tol,k1,k2,k3,omega);
time_taken=toc;
end
end
end
end
end
end
The Jacobian Expllicit function for Transient equation:
function u = Jacobiansolver(T,x,y,nx,ny,nt,error,tol,k1,k2);
for k=length(nt);
jacobi_iter=1;
Told=T;
while (error > tol)
for i=2:nx-1;
j=2:ny-1;
H=Told(i+1,j)-(2*Told(i,j))+Told(i-1,j);
V=Told(i,j+1)-(2*Told(i,j))+Told(i,j-1);
T(i,j)= Told(i,j)+(k1*H)+(k2*V);
end
error = max(max(abs(Told-T)));
Told=T;
jacobi_iter=jacobi_iter+1;
[X,Y]=meshgrid(x,y);
figure(1)
C=contourf(X,Y,T);
set(gca,\'YDIR\',\'reverse\');
colormap(jet);
clabel(C,\'fontsize\',14);
xlabel(\'grid points\');
ylabel(\'Temperature\');
title_text=sprintf(\'No of iterationsfor 2D transient Jacobian method=%d\',jacobi_iter);
title(title_text);
end
endGauss Siedel Function explicit scheme for Transient equation:
function T = GaussSeidelsolver(T,x,y,nx,ny,nt,error,tol,k1,k2);
for k=length(nt);
gauss_iter=1;
Told=T;
while (error > tol)
for i=2:nx-1;
j=2:ny-1;
H=T(i+1,j)-(2*T(i,j))+T(i-1,j);
V=T(i,j+1)-(2*T(i,j))+T(i,j-1);
T(i,j)= Told(i,j)+(k1*H)+(k2*V);
end
error = max(max(abs(Told-T)));
Told=T;
gauss_iter=gauss_iter+1;
[X,Y]=meshgrid(x,y);
figure(1)
C=contourf(X,Y,T);
set(gca,\'YDIR\',\'reverse\');
colormap(jet);
clabel(C,\'fontsize\',14);
xlabel(\'grid points\');
ylabel(\'Temperature\');
title_text=sprintf(\'No of iterations for 2D transient=%d\',gauss_iter);
title(title_text);
end
endThe function for SOR scheme on Gauss Siedel for Transient:
function T = SORgausssolver(T,x,y,nx,ny,nt,error,tol,k1,k2);
for k=length(nt);
Told=T;
dx=1/49;
omega=2/(1+sin(pi*dx));
gauss_iter=1;
while (error > tol)
for i=2:nx-1;
j=2:ny-1;
H=T(i+1,j)-(2*T(i,j))+T(i-1,j);
V=T(i,j+1)-(2*T(i,j))+T(i,j-1);
T(i,j)= ((omega)*(T(i,j)+(k1*H)+(k2*V)))-((omega-1)*T(i,j));
end
error = max(max(abs(Told-T)));
Told=T;
gauss_iter=gauss_iter+1;
[X,Y]=meshgrid(x,y);
figure(1)
C=contourf(X,Y,T);
set(gca,\'YDIR\',\'reverse\');
colormap(jet);
clabel(C,\'fontsize\',14);
xlabel(\'grid points\');
ylabel(\'Temperature\');
title_text=sprintf(\'No of iterationsfor 2D transient=%d\',gauss_iter);
title(title_text);
end
endfunction for Implicit Scheme for Jacobian on Transient case:
function T = Jacobian_implicit_solver(T,Told,T_prev_dt,x,y,nx,ny,nt,error,tol,k1,k2,k3);;
for k=1:nt
implicitjacob_iter=1
while (error > tol)
Told=T;
T_prev_dt=Told;
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)= (k3*T_prev_dt(i,j))+(k1*k3*H)+(k3*k2*V);
end
end
error = max(max(abs(Told-T)));
Told=T;
[X,Y]=meshgrid(x,y);
figure(1)
C=contourf(X,Y,T);
set(gca,\'YDIR\',\'reverse\');
colormap(jet);
clabel(C,\'fontsize\',14);
xlabel(\'grid points\');
ylabel(\'Temperature\');
title_text=sprintf(\'No of iterations for 2D transient Jacobian Implicit Scheme=%d\',implicitjacob_iter);
title(title_text);
pause(0.000003)
implicitjacob_iter=implicitjacob_iter+1;
end
T_prev_dt=T;
end
end Gauss Siedel Implicit Solver for function for Transient Case:
function T = GaussSiedel_implicit_solver(T,Told,T_prev_dt,x,y,nx,ny,nt,error,tol,k1,k2,k3);
for k=1:nt
implicitgauss_iter=1;
while (error > tol)
Told=T;
T_prev_dt=Told;
for i=2:nx-1;
for j=2:ny-1;
H(i,j)=T(i+1,j)+T(i-1,j);
V(i,j)=T(i,j+1)+T(i,j-1);
T(i,j)= (k3*T_prev_dt(i,j))+(k1*k3*H(i,j))+(k2*k3*V(i,j));
end
end
error = max(max(abs(Told-T)));
Told=T;
[X,Y]=meshgrid(x,y);
figure(1)
C=contourf(X,Y,T);
set(gca,\'YDIR\',\'reverse\');
colormap(jet);
clabel(C,\'fontsize\',14);
xlabel(\'grid points\');
ylabel(\'Temperature\');
title_text=sprintf(\'No of iterations for 2D transient=%d\',implicitgauss_iter);
title(title_text);
pause(0.003)
implicitgauss_iter=implicitgauss_iter+1;
end
T_prev_dt=T;
end
endSOR on gauss Siedel implicit function for Transient case:
function T = SORGaussSiedel_implicit_solver(T,Told,T_prev_dt,x,y,nx,ny,nt,error,tol,k1,k2,k3,omega);
tic
for k=1:nt
implicitsorjacob_iter=1;
while (error > tol)
Told=T;
T_prev_dt=Told;
for i=2:nx-1;
for j=2:ny-1;
H(i,j)=T(i+1,j)+T(i-1,j);
V(i,j)=T(i,j+1)+T(i,j-1);
T(i,j)= omega*(((k1*k3*H(i,j))+(k2*k3*V(i,j))+k3*T_prev_dt(i,j)))-((omega-1)*(T_prev_dt(i,j)));
end
end
error = max(max(abs(Told-T)));
Told=T;
[X,Y]=meshgrid(x,y);
figure(1)
C=contourf(X,Y,T);
set(gca,\'YDIR\',\'reverse\');
colormap(jet);
clabel(C,\'fontsize\',14);
xlabel(\'grid points\');
ylabel(\'Temperature\');
title_text=sprintf(\'No of iterations for 2D transient for SOR on Gauss Siedel Scheme=%d\',implicitsorjacob_iter);
title(title_text);
pause(0.003)
implicitsorjacob_iter=implicitsorjacob_iter+1;
end
T_prev_dt=Told;
end
time_taken=toc;
endGraphs For the Transient case:
Jacobian Method Explicit:
Gauss Siedel Method Explicit:
SOR on Gauss Siedel Explicit Transient:
Graphs for Implicit Schemes:
Gauss Siedel Implcicit Schemes for transient case:
SOR on Gauss Siedel Implicit Scheme for transient case:
Conclusion:
In the steady state equation SOR has the less number of iterations and converges pretty quickly than Jacobian and gauss Siedel method.
The Explicit Scheme for the transinet scheme takes less number of iteartions when compared to the Implicit scheme.
This might depend on the number of grids and time step selected for the problem.
But it can be concluded from above that SOR is very useful in converging the solution quickly which can prove to be cost effective and save a lot of computational time.