Week 5.1 - Mid term project - Solving the steady and unsteady 2D heat conduction problem
STEADY STATE HEAT CONDUCTION:
In heat conduction under steady conditions, the temperature of a body at any point does not change with time. Thus, it can be concluded that this type of heat conduction is a function of position only.
The formula for steady state heat conduction without heat generation is given by the formula:
We make use of numerical methods making use of boundary elements.The solutions are obtained at discrete points in our square domain.
Matlab Program:
clear all
close all
clc
%initial conditions
L=1; %length of the domain
nx=input('For the uniformly space grid please enter the number of grid points:');
ny=nx;
x=linspace(0,L,nx);
y=linspace(L,0,ny);
dx=L/(nx-1); %providing the spacing between the grid points along x axis
dy=L/(ny-1); %providing the spacing between the grid points along y axis
error=9e9;
tol=1e-4;
%Defining Boundary Conditions
T_L=400;
T_R=800;
T_T=600;
T_B=900;
%Temperature Matrix
T=298*ones(nx,ny);
T(1,:)=T_B;
T(end,:)=T_T;
T(:,1)=T_L;
T(:,end)=T_R;
%Edge Refining
T(1,1)=(900+400)/2;
T(1,end)=(900+800)/2;
T(end,1)=(400+600)/2;
T(end,end)=(800+600)/2;
%Calculation of relaxation factor
%omega=2/(1+sin(pi*dx));
omega=1.2;
%Creating a copy of old values
[X,Y]=meshgrid(x,y);
Told=T;
iterative_solver=input('1) Jacobi 2) Gauss Seidel 3) SOR \nEnter an option for the type of solver you would like to use: \n');
k=2/dx^2 + 2/dy^2;
tic
if iterative_solver==1;
jacobi_iter=1;
while(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
error=max(max(abs(Told-T)));
Told=T;
jacobi_iter=jacobi_iter+1;
figure(1)
[A,B]=contourf(X,Y,T);
set(gca,'YDIR','reverse')
colorbar
colormap(jet)
clabel(A,B);
text=sprintf('Jacobi Iterations = %d',jacobi_iter);
ylabel('Y')
title(text)
pause(0.03)
end
end
if iterative_solver==2;
gs_iter=1;
while(error>tol)
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)=term1+term2;
end
end
error=max(max(abs(Told-T)));
Told=T;
gs_iter=gs_iter+1;
figure(1)
[A,B]=contourf(X,Y,T);
set(gca,'YDIR','reverse')
colorbar
colormap(jet)
clabel(A,B);
text=sprintf('Gauss Seidel Iterations = %d ',gs_iter);
title(text)
ylabel('Y')
pause(0.03)
end
end
if iterative_solver==3;
sor_iter=1;
while(error>tol)
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;
figure(1)
[A,B]=contourf(X,Y,T);
set(gca,'YDIR','reverse')
colorbar
colormap(jet)
clabel(A,B);
text=sprintf('SOR Iterations = %d ',sor_iter);
title(text)
ylabel('Y')
pause(0.03)
end
end
Time_Taken=tocOUTPUT:
Jacobi Iteration
Gauss Seidel
SOR
Transient state heat cconduction
In heat conduction under unsteady state conditions, the temperature of a body at any point varies with time as well as position in one- dimensional and multidimensional systems. This type of heat conduction is also known as transient heat conduction. The 2D transient heat conduction equation is given by:
In Method of Charactersitics we learned that when having a time derivative in the governing equation we make use of the time marching solution. This solution is implemented by making use of numerical schemes at each and every discrete points on the square domain.
Explicit Method Program
clear all
close all
clc
L=1;
nx=input('For the uniformly space grid please enter the number of grid points:');
ny=nx;
x=linspace(0,L,nx);
y=linspace(L,0,ny);
dx=L/(nx-1);
dy=L/(ny-1);
dt=0.001;
%t=0.1;
%nt=t/dt;
error=9e9;
tol=1e-4;
T_L=400;
T_R=800;
T_T=600;
T_B=900;
T=298*ones(nx,ny);
T(1,:)=T_B;
T(end,:)=T_T;
T(:,1)=T_L;
T(:,end)=T_R;
T(1,1)=(900+400)/2;
T(1,end)=(900+800)/2;
T(end,1)=(400+600)/2;
T(end,end)=(800+600)/2;
alpha=1.4; %Diffusivity
k1=alpha*dt/dx^2;
k2=alpha*dt/dy^2;
[X,Y]=meshgrid(x,y);
Told=T;
CFL=k1+k2 %To verify if the solution is stable within the limit CFL(x)+CFL(y)=0.5
unsteady_explicit_iter=1;
tic
for nt=1:1400
while(error>tol)
for i=2:nx-1
for j=2:ny-1
term1=Told(i-1,j)-2*Told(i,j)+Told(i+1,j);
term2=Told(i,j-1)-2*Told(i,j)+Told(i,j+1);
T(i,j)=Told(i,j)+k1*term1+k2*term2;
end
end
error=max(max(abs(Told-T)));
Told=T;
unsteady_explicit_iter=unsteady_explicit_iter+1;
figure(1)
[A,B]=contourf(X,Y,T);
set(gca,'YDIR','reverse')
colorbar
colormap(jet)
clabel(A,B);
text=sprintf('Unsteady Explicit Iterations = %d',unsteady_explicit_iter);
ylabel('Y')
title(text)
pause(0.03)
end
end
Time_Taken=toc
Output:
Implicit Method Program
clear all
close all
clc
L=1;
nx=input('For the uniformly space grid please enter the number of grid points:');
ny=nx;
x=linspace(0,L,nx);
y=linspace(L,0,ny);
dx=1/(nx-1);
dy=1/(ny-1);
alpha=1.4;
dt=0.001;
%error=9e9
tol=1e-4;
T_L=400;
T_R=800;
T_T=600;
T_B=900;
T=298*ones(nx,ny);
T(1,:)=T_B;
T(end,:)=T_T;
T(:,1)=T_L;
T(:,end)=T_R;
T(1,1)=(900+400)/2;
T(1,end)=(900+800)/2;
T(end,1)=(400+600)/2;
T(end,end)=(800+600)/2;
k1=alpha*dt/dx^2;
k2=alpha*dt/dy^2;
CFL=k1+k2
Told=T;
T_prev_dt=T;
iter=1;
omega=1.1;
iterative_solver=input('1) Jacobi 2) Gauss Seidel 3) SOR \nEnter an option for the type of solver you would like to use: \n');
tic
for nt=1:1400
if iterative_solver==1
term1=(1+2*k1+2*k2)^(-1);
term2=k1*term1;
term3=k2*term1;
error=9e9;
while(error>tol)
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) = (T_prev_dt(i,j)*term1)+(H*term2)+(V*term3);
end
end
error=max(max(abs(Told-T)));
Told=T;
iter=iter+1;
end
end
if iterative_solver==2
term1=(1+2*k1+2*k2)^(-1);
term2=k1*term1;
term3=k2*term1;
error=9e9;
while(error>tol)
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) = (T_prev_dt(i,j)*term1)+(H*term2)+(V*term3);
end
end
error=max(max(abs(Told-T)));
Told=T;
iter=iter+1;
end
end
if iterative_solver==3
term1=(1+2*k1+2*k2)^(-1);
term2=k1*term1;
term3=k2*term1;
error=9e9;
while(error>tol)
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*((T_prev_dt(i,j)*term1)+(H*term2)+(V*term3));
end
end
error=max(max(abs(Told-T)));
Told=T;
iter=iter+1;
end
end
T_prev_dt=T;
end
Time_Taken=toc
figure(1)
[P Q]=contourf(x,y,T);
set(gca,'YDIR','reverse')
colorbar;
colormap(jet);
clabel(P,Q);
text=sprintf('Implicit Iterations = %d',iter);
ylabel('Y')
title(text)
pause(0.03)Output:
Jacobi :
Gauss seidel:
SOR :
Result:
In steady state methods, there is no time loop and the results converge within a few iterations. We include time loop for implicit methods and it runs for a large number of iterations.
Transient solvers take a lot of time for computation as it is a function of both time and position, when comparing to steady state solver which is only a function of position.
The SOR provides us the answer with less number of iterations and takes less computation time.