Week 5 - Mid term project - Solving the steady and unsteady 2D heat conduction problem
The objective was to solve the steady state and transient 2D heat conduction equation.
Steady State Equation:
The steady state equation is given as:
`\(partial^2T)\/(\partialx^2) + \(partial^2T)\/(\partialy^2)=0`
So by numerically discretising the equation we have,
`Tp = (1/K)*[(T_L+T_R)/(\Deltax^2) + (T_T+T_B)/(\Deltay^2)]`
where K is `2((\Deltax^2+\Deltay^2)/(\Deltax^2*\Deltay^2))`
The domain is assumed to be a unit square. The temperature boundary conditions are:
1. T_Left:400K
2. T_Top:600K
3. T_Right:800K
4. T_Bottom:900K
Initially all the points except boundary are assumed to be at 298K. The absolute error criteria is 1e-4.
Following are the results obtained from the 3 solvers:
Jacobi Solver:
The jacobi solver uses the old values to update the current values.
As can be seen the temperature reached steady state in 826 iterations with simulation time of 0.017s.
The above graph shows the drop in error as iterations progress. The rate is slow at first but increases as iterations progress.
Gauss-Seidel Solver:
The Gauss-Seidel Solver uses the new values updated upto last iteration to update current values.
The above contour shows the steady state temperature contour plot as solved by Gauss Seidel solver. The solver converged in 440 iterations with a simulation time of 0.012s, which is less than Jacobi Solver.
As can be seen in the error vs iterations plot,the error drops faster in case of Gauss-Seidel solver as iterations progress.
SOR Solver:
The SOR uses relaxation parameter omega which is set to 1.7313 to achieve faster convergence.
The SOR solver reaches the solution in 59 iterations with simulation time= 0.005s, which is the least of the 3 solvers. This is due to the fact that over-relaxation was used with parameter omega=1.7313. The optimal parameter was obtained by trying omega from 0.5 to 2.
As is evident from error vs iteration graph, the error drops very fast due to over-relaxation.
Transient Equation:
The transient 2D heat conduction equation is given as:
`\(partialT)/(\partialt) = alpha((\partial^2T)/(\partialx^2)+(\partial^2T)/(\partialy^2))`
So by numerically discretising the equation we have,
Explicit:`T_p^(n+1)= T_p^(n) + alpha\Deltat((T_R-2T_p+T_L)/(\Deltax^2)+(T_T-2T_p+T_B)/(\Deltay^2))`
Implicit:`T_p^(n+1) = (T_p^(n)+k_1*(T_L^(n+1)+T_R^(n+1))+k_2*(T_T^(n+1)+T_B^(n+1)))/(1+2*k_1+2*k_2)`
where: k1=`alpha(\Deltat)/(\Deltax^2)`
k2=`alpha(\Deltat)/(\Deltay^2)`
Explicit: In case of explicit formulation, the terms on right side are known, so there is no need of iterative solver.
As can be seen the Explicit solver converges after 820 timesteps or iterations with simulation time of 0.024s.
Implicit: In case of implicit formulation there are unknown terms present on both sides of equation, hence it gives a coupled system of equations, which requires an iterative solver to reach desired solution.
Jacobi Solver:
The Jacobi implicit solver uses the older values to iterate in each timestep.
As the implicit solver is unconditionally stable due to it's formulation, the timestep value is taken very large as 1. As can be seen it reaches solution in 1770 iterations with simulation time of 0.028s.
GS Solver:
The GS solver uses the updated values to solve for next iteration within each timestep.
As can be seen the GS implicit solver converges in 986 iterations with simulation time=0.013s which is still lesser than Jacobi solver.
SOR Solver:
The SOR solver uses the relaxation parameter omega which is set at 1.7313.
As can be seen the solution converges in 256 iterations with a simulation time of 0.012s which is least of the solvers for implicit formulation
%Gaurav Gokhale
close all
clear all
clc
%Steady State
% omega=linspace(0.1,2,100); % loop for value of omega
% iter_final=1000;
% for k=1:length(omega)
% omega1=omega(k);
%Inputs
L=1;
nx=20;
ny=20;
method=3; % method=1 for jacobi, 2 for Gauss-Seidel & 3 for SOR
tol=1e-4;
error=1;
%Calculations
dx=L/nx;
dy=L/ny;
x=linspace(0,1,20);
y=linspace(0,1,20);
[X,Y]=meshgrid(x,y);
k=2*(dx^2+dy^2)/(dx*dy)^2;
T=298*ones(nx);
%BCs
T(1,:)=900;
T(ny,:)=600;
T(:,1)=400;
T(:,nx)=800;
T(1,1)=650;
T(1,nx)=850;
T(ny,1)=500;
T(ny,nx)=700;
Told=T;
%Jacobi
if (method==1)
tic
iter=1;
while(error>tol)
for i=2:nx-1
for j=2:nx-1
term1=(Told(i,j-1)+Told(i,j+1))/dx^2;
term2=(Told(i-1,j)+Told(i+1,j))/dy^2;
T(i,j)=(1/k)*(term1+term2);
end
end
% figure(1)
% [C,h]=contourf(X,Y,T);
% colormap jet
% clabel(C,h,'FontSize',10,'Color','black');
% title_text=sprintf('Iteration Number %d',iter);
% title(title_text)
%pause(0.03)
error(iter)=max(max(abs(T-Told)));
Told=T;
% figure(2)
% iter1=1:iter;
% loglog(iter1,error)
% xlabel('Iterations')
% ylabel('Error')
% title('Error Vs Iterations Jacobi Solver')
iter=iter+1;
end
sim_time=toc;
figure(1)
[C,h]=contourf(X,Y,T);
colormap jet
clabel(C,h,'FontSize',10,'Color','black');
contourcbar
xlabel('X-axis');
ylabel('Y-axis')
title_text=sprintf('Steady State Jacobi Sovler, Iterations=%d, simulation time=%8.3f s',iter-1,sim_time);
title(title_text);
end
%GS
if method==2
tic
iter=1;
while(error>tol)
for i=2:nx-1
for j=2:nx-1
term1=(T(i,j-1)+T(i,j+1))/dx^2;
term2=(T(i-1,j)+T(i+1,j))/dy^2;
T(i,j)=(1/k)*(term1+term2);
end
end
% figure(1)
% [C,h]=contourf(X,Y,T);
% colormap jet
% clabel(C,h,'FontSize',10,'Color','black');
% title_text=sprintf('Iteration Number %d',iter);
% title(title_text)
% pause(0.03)
error(iter)=max(max(abs(T-Told)));
Told=T;
% figure(2)
% iter1=1:iter;
% loglog(iter1,error)
% xlabel('Iterations')
% ylabel('Error')
% title('Error Vs Iterations GS Solver')
iter=iter+1;
end
sim_time=toc;
figure(1)
[C,h]=contourf(X,Y,T);
colormap jet
clabel(C,h,'FontSize',10,'Color','black');
contourcbar
xlabel('X-axis');
ylabel('Y-axis')
title_text=sprintf('Steady State Gauss-Seidel Sovler, Iterations=%d, simulation time=%8.3f s',iter-1,sim_time);
title(title_text);
end
%SOR
if method==3
tic
iter=1;
omega1=1.7313;
while(error>tol)
for i=2:nx-1
for j=2:nx-1
term1=(T(i,j-1)+T(i,j+1))/dx^2;
term2=(T(i-1,j)+T(i+1,j))/dy^2;
T(i,j)=Told(i,j)*(1-omega1)+omega1*(1/k)*(term1+term2);
end
end
% figure(1)
% [C,h]=contourf(X,Y,T);
% colormap jet
% clabel(C,h,'FontSize',10,'Color','black');
% title_text=sprintf('Iteration Number %d',iter);
% title(title_text)
% pause(0.03)
error(iter)=max(max(abs(T-Told)));
Told=T;
% figure(2)
% iter1=1:iter;
% loglog(iter1,error)
% xlabel('Iterations')
% ylabel('Error')
% title('Error Vs Iterations SOR Solver')
%
iter=iter+1;
% if iter>iter_final
% break
% end
end
% if iter<iter_final
% omega_final=omega1
% iter_final=iter
% end
sim_time=toc;
figure(1)
[C,h]=contourf(X,Y,T);
colormap jet
clabel(C,h,'FontSize',10,'Color','black');
colorbar;
xlabel('X-axis')
ylabel('Y-axis')
title_text=sprintf('Steady State SOR Sovler, Iterations=%d, simulation time=%8.3f s',iter-1,sim_time);
title(title_text);
end
% end
%Transient
close all
clear all
clc
% omega=linspace(0.1,2,100); % loop for value of omega
% iter_final=1000;
% for l=1:length(omega)
% omega1=omega(l);
%Inputs
L=1;
nx=20; %grid points
ny=20;
alpha=1.4;
m=2; %m=1 for explicit and m=2 for implicit
method=3; %method=1 for jacobi 2 for Gauss-Seidel & 3 for SOR
tol=1e-4; %absolute error value
error=1;
%Calculations
dx=L/nx;
dy=L/ny;
dt=4.5e-4;
t=1;
k1=alpha*dt/(dx)^2;
k2=alpha*dt/(dy)^2;
x=linspace(0,1,20);
y=linspace(0,1,20);
[X,Y]=meshgrid(x,y);
T=298*ones(nx);
%BCs
T(1,:)=900;
T(ny,:)=600;
T(:,1)=400;
T(:,nx)=800;
T(1,1)=650;
T(1,nx)=850;
T(ny,1)=500;
T(ny,nx)=700;
Told=T;
T_prev=T;
%explicit
if m==1
tic
for k=1:1400 %timesteps
for i=2:ny-1
for j=2:nx-1
term1=T_prev(i,j-1)-2*T_prev(i,j)+T_prev(i,j+1);
term2=T_prev(i+1,j)-2*T_prev(i,j)+T_prev(i-1,j);
T(i,j)=T_prev(i,j)+k1*term1+k2*term2;
end
end
% figure(2)
% residual(k)=max(max(abs(T_prev-T)));
% loglog(residual);
% xlabel('timesteps')
% ylabel('residual')
% title_text1=sprintf('residual Vs timesteps for Explicit for dt=%d',dt);
% title(title_text1)
%
if max(max(abs(T_prev-T)))<1e-4
break
end
T_prev=T;
% figure(1)
% [C,h]=contourf(X,Y,T);
% colormap jet
% clabel(C,h,'FontSize',10,'Color','black');
% title_text=sprintf('Timestep %d',k);
% title(title_text)
% pause(0.03)
end
sim_time=toc;
figure(1)
[C,h]=contourf(X,Y,T);
colormap jet
clabel(C,h,'FontSize',10,'Color','black');
colorbar;
xlabel('X-axis')
ylabel('Y-axis')
title_text=sprintf('Unsteady State Explicit Sovler, timesteps=%d, simulation time=%8.3f s',k,sim_time);
title(title_text);
end
%implicit
if m==2
dt=1;
k1=alpha*dt/(dx)^2;
k2=alpha*dt/(dy)^2;
%Jacobi
if (method==1)
tic
iter=1;
for k=1:1400 %timesteps
error=1;
T_prev=T;
while (error>tol)
term1=(1+2*k1+2*k2)^-1;
term2=k1*term1;
term3=k2*term1;
for i=2:nx-1
for j=2:nx-1
H=Told(i,j-1)+Told(i,j+1);
V=Told(i+1,j)+Told(i-1,j);
T(i,j)=T_prev(i,j)*term1+H*term2+V*term3;
end
end
error=max(max(abs(T-Told)));
Told=T;
iter=iter+1;
end
% figure(1)
% [C,h]=contourf(X,Y,T);
% colormap jet
% clabel(C,h,'FontSize',10,'Color','black');
% title_text=sprintf('Timestep %d',k);
% title(title_text)
% pause(0.03)
% figure(2)
% residual(k)=max(max(abs(T_prev-T)));
% loglog(residual);
% xlabel('timesteps')
% ylabel('residual')
% title('residual Vs timesteps for Implicit Jacobi Solver')
if max(max(abs(T_prev-T)))<1e-4
break
end
end
sim_time=toc;
figure(1)
[C,h]=contourf(X,Y,T);
colormap jet
clabel(C,h,'FontSize',10,'Color','black');
colorbar;
xlabel('X-axis')
ylabel('Y-axis')
title_text=sprintf('Unsteady State Implicit Jacobi Sovler, iterations=%d, simulation time=%8.3f s',iter,sim_time);
title(title_text);
end
%GS
if (method==2)
iter=1;
tic
for k=1:1400 %timesteps
error=1;
T_prev=T;
while (error>tol)
term1=(1+2*k1+2*k2)^-1;
term2=k1*term1;
term3=k2*term1;
for i=2:nx-1
for j=2:nx-1
H=T(i,j-1)+T(i,j+1);
V=T(i+1,j)+T(i-1,j);
T(i,j)=T_prev(i,j)*term1+H*term2+V*term3;
end
end
error=max(max(abs(T-Told)));
Told=T;
iter=iter+1;
end
% figure(1)
% [C,h]=contourf(X,Y,T);
% colormap jet
% clabel(C,h,'FontSize',10,'Color','black');
% title_text=sprintf('Timestep %d',k);
% title(title_text)
% pause(0.03)
%
% figure(2)
% residual(k)=max(max(abs(T_prev-T)));
% loglog(residual);
% xlabel('timesteps')
% ylabel('residual')
% title('residual Vs timesteps for Implicit GS Solver')
if max(max(abs(T_prev-T)))<1e-4
break
end
end
sim_time=toc;
figure(1)
[C,h]=contourf(X,Y,T);
colormap jet
clabel(C,h,'FontSize',10,'Color','black');
colorbar;
xlabel('X-axis')
ylabel('Y-axis')
title_text=sprintf('Unsteady State Implicit Gauss-Seidel Sovler, iterations=%d, simulation time=%8.3f s',iter,sim_time);
title(title_text);
end
%SOR
if (method==3)
iter=1;
tic
for k=1:1400 %timesteps
error=1;
T_prev=T;
omega1=1.7313;
while(error>tol)
term1=(1+2*k1+2*k2)^-1;
term2=k1*term1;
term3=k2*term1;
for i=2:nx-1
for j=2:nx-1
H=(1-omega1)*(Told(i,j-1)+Told(i,j+1))+omega1*(T(i,j-1)+T(i,j+1));
V=(1-omega1)*(Told(i+1,j)+Told(i-1,j))+omega1*(T(i+1,j)+T(i-1,j));
T(i,j)=T_prev(i,j)*term1+H*term2+V*term3;
end
end
error=max(max(abs(T-Told)));
Told=T;
iter=iter+1;
% if iter>iter_final
% break
% end
end
% figure(1)
% [C,h]=contourf(X,Y,T);
% colormap jet
% clabel(C,h,'FontSize',10,'Color','black');
% title_text=sprintf('Timestep %d',k);
% title(title_text)
% % pause(0.03)
%
% residual(k)=max(max(abs(T_prev-T)));
% figure(2)
% loglog(residual);
% xlabel('Timesteps')
% ylabel('Residual')
% title('Residual Vs Timestep Implicit SOR')
if max(max(abs(T_prev-T)))<1e-4
break
end
end
% if iter<iter_final
% omega_final=omega1
% iter_final=iter
% end
sim_time=toc;
figure(1)
[C,h]=contourf(X,Y,T);
colormap jet
clabel(C,h,'FontSize',10,'Color','black');
title_text=sprintf('Unsteady State Implicit SOR Sovler, iterations=%d, simulation time=%8.3f s',iter,sim_time);
title(title_text);
end
end
% end