Project

Steady and Unsteady heat conduction using MATLab

NNikhil Darekar5 Nov 2018

#Code For Steady State Conduction:

% Steady conduction
% Square Grid
clc
clear all
close all
L=1;
nx=10;
ny=10;
x=linspace(0,1,nx);
y=linspace(0,1,ny);
dx=x(2)-x(1);
dy=y(2)-y(1);
T_initial=ones(nx,ny);
T_initial(1,1:ny)=900;
T_initial(1:nx,ny)=800;
T_initial(nx,1:ny)=600;
T_initial(1:nx,1)=400;
T_initial(nx,1)=500;
T_initial(1,1)=650;
T_initial(1,ny)=850;
T_initial(nx,ny)=700;
k1= (dy^2)/(2*((dy^2)+(dx^2)));
k2= (dx^2)/(2*((dy^2)+(dx^2)));
error=9e9;
tol=1e-4;
Told=T_initial;
iter=0;
tic
z=input(\'enter your choice 1=Jacobian 2=gauss seidel 3=sor =\');
if z==1
while(error >tol)
for i=2:1:(nx-1)

for j=2:1:(ny-1)

T_initial(i,j)=k1*(Told((i-1),j)+Told(i+1,j))+k2*(Told(i,j-1)+Told(i,j+1));

end

end
error=max(max(abs(Told - T_initial)));
Told=T_initial;
iter=iter+1;
contourf(x,y,T_initial,\'ShowText\',\'on\');
set(gca,\'YDIR\',\'reverse\');
colormap(jet);
colorbar;
title(sprintf(\'No of Iteration=%d\',iter));
xlabel(\'length of slab\')
ylabel(\'height of slab\')
hold on
pause(0.1)
end
time=toc
fprintf(\'\\n no of iterations taken=%d\',iter);
fprintf(\'\\n Time taken=%g\',time);
fprintf(\'\\n the error accounted=%g\',error);
end

tic
if z==2
itergauss=0;
while(error >tol)
for i=2:1:(nx-1)

for j=2:1:(ny-1)

T_initial(i,j)=k1*(T_initial((i-1),j)+T_initial(i+1,j))+k2*(T_initial(i,j-1)+T_initial(i,j+1));

end

end
error=max(max(abs(Told - T_initial)));
Told=T_initial;
itergauss=itergauss+1;
contourf(x,y,T_initial,\'edgecolor\',\'none\');
set(gca,\'YDIR\',\'reverse\');
colormap(jet);
colorbar
xlabel(\'length of slab\')
ylabel(\'height of slab\')
title(sprintf(\'No of Iteration=%d\',itergauss));
hold on
pause(0.1)
end
end

% Steady state SOR;

tic

if z==3
k1= (dy^2)/(2*((dy^2)+(dx^2)));
k2= (dx^2)/(2*((dy^2)+(dx^2)));
T_iter=T_initial;
omeg=2/(1+(sin(pi*dx)));
itersor=0;
while(error >tol)
for i=2:1:(nx-1)

for j=2:1:(ny-1)

T_initial(i,j)=((omeg*((k1*(T_initial((i-1),j)+T_initial(i+1,j)))+(k2*(T_initial(i,j-1)+T_initial(i,j+1))))))+((1-omeg)*T_iter(i,j));

end

end
error=max(max(abs(Told - T_initial)));
Told=T_initial;
T_iter=T_initial;
itersor=itersor+1;
contourf(x,y,T_initial,\'edgecolor\',\'none\');
set(gca,\'YDIR\',\'reverse\');
colormap(jet);
colorbar
xlabel(\'length of slab\')
ylabel(\'height of slab\')
title(sprintf(\'No of Iteration=%d\',itersor));
hold on
pause(0.1)
end
time=toc
fprintf(\'\\n no of iterations taken=%d\',itersor);
fprintf(\'\\n Time taken=%g\',time);
fprintf(\'\\n the error accounted=%g\',error);
end

% Steady state SOR;

tic

Results for Steady State Conduction:

1...............................Jacobi.......................................................

$\"\"$

2...........................Gauss Seidel..........................................................

$\"\"$

3.................................SOR...............................................................

$\"\"$

Explanations:

From above contour plots we can observe that simulation iteration taken for jacobi method was highest.

Gauss seidel take least iteration as compared to jacobi but higher than the SOR.

From above we can conclude that for reaching steady state the SOR method was faster.

# For Unsteady Conduction Program

% unsteady conduction
% Square Grid
clc
clear all
close all
L=1;
nx=10;
ny=10;
x=linspace(0,1,nx);
y=linspace(0,1,ny);
dx=x(2)-x(1);
dy=y(2)-y(1);
T_initial=ones(nx,ny);
T_initial(1,1:ny)=900;
T_initial(1:nx,ny)=800;
T_initial(nx,1:ny)=600;
T_initial(1:nx,1)=400;
T_initial(nx,1)=500;
T_initial(1,1)=650;
T_initial(1,ny)=850;
T_initial(nx,ny)=700;
k1= (dy^2)/(2*((dy^2)+(dx^2)));
k2= (dx^2)/(2*((dy^2)+(dx^2)));
error=9e9;
tol=1e-4;
Told=T_initial;
iter=0;
tic
dt=32;
alpha=1e-4;
k3=(alpha*dt)/(dx^2);
k4=(alpha*dt)/(dy^2);
z=input(\'enter your choice 1=Jacobian 2+gauss seidel 3=sor=\');
if z==1
for i=1:200

for i=2:1:(nx-1)

for j=2:1:(ny-1)

T_initial(i,j)=k3*(Told((i-1),j)+Told(i+1,j)-2*Told(i,j))+k2*(Told(i,j-1)+Told(i,j+1)-2*Told(i,j))+Told(i,j);

end

end
error1=max(max(Told-T_initial));
Told=T_initial;
iter=iter+1;
contourf(x,y,T_initial,\'edgecolor\',\'none\');
set(gca,\'YDIR\',\'reverse\');
colormap(jet);
colorbar
xlabel(\'length of slab\')
ylabel(\'height of slab\')
title(sprintf(\'No of Iteration=%d\',iter));
hold on
pause(0.1)
end
time=toc
fprintf(\'\\n no of iterations taken=%d\',iter);
fprintf(\'\\n Time taken=%g\',time);
fprintf(\'\\n the error accounted=%g\',error1);
end

if z==2
%Unstedy Exlicit Method.
k1= alpha*dt/((dx^2));
k2= alpha*dt/((dy^2));
term1=((2*k1)+(2*k2)+1)^(-1);
term2=term1*k1;
term3=term1*k2;
Told=T_initial;
T_olddt=Told;
iterun=0;
for i=1:200
while (error> tol)

for i=2:nx-1
for j=2:ny-1
H=(Told((i-1),j)+Told(i+1,j))*term2;
V=(Told(i,j-1)+Told(i,j+1))*term3;
T_initial(i,j)=H+V+T_olddt(i,j)*term1;
end
end
error1=max(max(abs(Told- T_initial)));
error=max(max(abs(Told- T_initial)));
Told=T_initial;
iterun=iterun+1;
contourf(x,y,T_initial,\'edgecolor\',\'none\');
colorbar
xlabel(\'length of slab\')
ylabel(\'height of slab\')
set(gca,\'YDIR\',\'reverse\');
colormap(jet);
colorbar;
title(sprintf(\'No of Iteration=%d\',iterun));
hold on
pause(0.1);
end
T_olddt=T_initial;
error=9e9;
end
time=toc
fprintf(\'\\n no of iterations taken=%d\',iterun);
fprintf(\'\\n Time taken=%g\',time);
fprintf(\'\\n the error accounted=%g\',error1);
end

if z==3
%Unstedy Implicit Method.
k1= alpha*dt/((dx^2));
k2= alpha*dt/((dy^2));
term1=((2*k1)+(2*k2)+1)^(-1);
term2=term1*k1;
term3=term1*k2;
Told=T_initial;
T_olddt=Told;
iterun=0;
for v=1:200
while (error> tol)

for i=2:nx-1
for j=2:ny-1
H=(T_initial((i-1),j)+Told(i+1,j))*term2;
V=(T_initial(i,j-1)+Told(i,j+1))*term3;
T_initial(i,j)=H+V+T_olddt(i,j)*term1;
end
end
error=max(max(abs(Told- T_initial)));
Told=T_initial;
iterun=iterun+1;
contourf(x,y,T_initial,\'edgecolor\',\'none\');
colorbar
xlabel(\'length of slab\')
ylabel(\'height of slab\')
set(gca,\'YDIR\',\'reverse\');
colormap(jet);
title(sprintf(\'No of Iteration=%d\',iterun));
hold on
pause(0.1)
end
T_olddt=T_initial;
error=9e9;
end
time=toc
fprintf(\'\\n no of iterations taken=%d\',iterun);
fprintf(\'\\n Time taken=%g\',time);
fprintf(\'\\n the error accounted=%g\',error);
end

if z==4
%Unstedy Implicit Method SOR.
omeg=2/(1+(sin(pi*dx)));
k1= (alpha*dt)/((dx^2));
k2= (alpha*dt)/((dy^2));
term1=1/((2*k1)+(2*k2)+1);
term2=term1*k1;
term3=term1*k2;
Told=T_initial;
T_olddt=T_initial;
iterun=0;
for v=1:200
while (error> tol)

for i=2:nx-1
for j=2:ny-1
H=(T_initial((i-1),j)+Told(i+1,j))*term2*omeg;
V=(T_initial(i,j-1)+Told(i,j+1))*term3*omeg;
T_initial(i,j)=H+V+(T_olddt(i,j)*(1-omeg))+(T_olddt(i,j)*term1*(omeg));
end
end
error=max(max(abs(Told- T_initial)));
Told=T_initial;

iterun=iterun+1;
contourf(x,y,T_initial,\'edgecolor\',\'none\');
colorbar
xlabel(\'length of slab\')
ylabel(\'height of slab\')
set(gca,\'YDIR\',\'reverse\');
colormap(jet);
title(sprintf(\'No of Iteration=%d\',iterun));
hold on
pause(0.1)
end
T_olddt=T_initial;
error=9e9;
end
time=toc
fprintf(\'\\n no of iterations taken=%d\',iterun);
fprintf(\'\\n Time taken=%g\',time);
fprintf(\'\\n the error accounted=%g\',error);
end

Output:

1.............Jacobi...................................................

$\"\"$

2.....................Gauss Seidel...........................................

$\"\"$

3..........................SOR............................ $\"\"$

From above results it can be seen that while time marching of solution is taking quite more time to reach steady state.

No of iterations are successively reduced from 2322 to 1570 to 1140.

Courant number C= alpha*dt/dx^2

=1e-4*32/0.1111^2

=0.2592

As the C<=1 provides the stability to the solution and convergence is achieved.

Steady and Unsteady heat conduction using MATLab

Ship work like this
in your own course.

Talk to a career counsellor

Ship work like thisin your own course.

Talk to a career counsellor

Ship work like this
in your own course.