FVM(MATLAB) code for the circular fin(conduction) problem
 

Dear All

I have some problem in FVM code please help me for the particular case i want to do k(thermal conductivity) and A(cross-sectional area) should be constant in circular fin.
close all
clc
%specify the number of control volumes
n=5; %number of control volumes
maxiter=200;%maximum iteration number
np1=n+1;
np2=n+2;
np3=n+3;
%%define calculation domain
tl=1; %total length of the cylinder
delx=tl/n;
dx=ones(1,np2);
dx=delx*dx;
%replace fictitious boundary volume size
dx(1)=1.0e-10;
dx(np2)=1.0e-10;
%assign x-coordinate
x(1)=0;
for m=1:np2
x(m+1)=x(m)+dx(m);
end
%define cross-sectional area
dia=0.25;
for i=1:np3
ac(i)=pi/4*dia^2; %constant cross section
end
%prescribe intitial temperatures for all control volumes
for i=1:np2
te(i)=100;
tep(i)=te(i);
end
%iteration for convergence
iter=0;
iflag=1;
%iteration loop for the convergence
while iflag==1 % end is at the end of program *********
%prescribe thermal conductivity
for i=1:np2
tk(i)=1;
end
%prescribe boundary temperature
te(1)=100;%at the left boundary given temperature
hf=0.2;% right boundary given flux
tf=20;
te(np2)=te(np1)-hf*(te(np2)-tf)/(tk(np1)/...
(.5*dx(np1)));
% te(np2)=100;
%initialize the coefficients for tdma matrix
%evaluate the diffusion conductance and source terms
for i=2:np1
% diffusion conductance
ke=tk(i)*tk(i+1)*(dx(i)+dx(i+1))/...
(dx(i)*tk(i+1)+dx(i+1)*tk(i));
de=2.0*ke*ac(i+1)/(dx(i)+dx(i+1));
ae=de;
kw=tk(i)*tk(i-1)*(dx(i)+dx(i-1))/...
(dx(i-1)*tk(i)+dx(i)*tk(i-1));
dw=2.0*kw*ac(i)/(dx(i-1)+dx(i));
aw=dw;
%source term evaluation
% sp=-1600;
% sc=4.768e5;
per=pi*dia;
sp=-hf*per/ac(i);
sc=hf*per/ac(i)*tf;
ap=ae+aw-sp*dx(i)*0.5*(ac(i+1)+ac(i));
b=sc*dx(i)*0.5*(ac(i)+ac(i+1));
%setting coefficients for tdma matrix
ta(i)=ap;
tb(i)=ae;
tc(i)=aw;
td(i)=b;
%modify using boundary conditions
if i==2
td(i)=td(i)+aw*te(1);
elseif i==np1
td(i)=td(i)+ae*te(np2);
end
end
%solve the simultaneous equations by using tdma
nq=n;
nqp1=nq+1;
nqm1=nq-1;
%forward substitution
beta(2)=tb(2)/ta(2);
alpa(2)=td(2)/ta(2);
for i=3:nqp1
beta(i)=tb(i)/(ta(i)-tc(i)*beta(i-1));
alpa(i)=(td(i)+tc(i)*alpa(i-1))/...
(ta(i)-tc(i)*beta(i-1));
end
%backward substitution
dum(nqp1)=alpa(nqp1);
for j=1:nqm1
i=nqp1-j;
dum(i)=beta(i)*dum(i+1)+alpa(i);
end
%end of tdma
%update the temperature
for i=2:np1
te(i)=dum(i);
end
%check the convergence
for i=1:np2
errote(i)=abs((te(i)-tep(i))/te(i));
end
error=1.0e-6;
if (max(errote)>error)
iter=iter+1;
tep=te;
iflag=1;
else
iflag=0;
end
if iter>maxiter
break
end
end % this end goes with the while iflag==1 at the top********
%compare with exact solution to make sure program works correctly
%locate the midpoint of control volumes
for i=1:np2
xc(i)=0.5*(x(i)+x(i+1));
end
%
m=sqrt(hf*per/(tk(1)*ac(1)));
thetab=te(1)-tf;
for i=1:np2
theta(i)=(cosh(m*(tl-xc(i)))+hf/(m*tk(1))*sinh(m*(tl-xc(i))))/...
(cosh(m*tl)+hf/(m*tk(1))*sinh(m*tl));
end
theta=thetab*theta;
teexact=theta+tf;
teexact=teexact'
%print the results
fprintf('iteration number is %i \n',iter)
disp('circular fin temperatures are')
fprintf('%9.3f \n',te')
%plot the result xc.vs.te
plot(xc,te,'-',xc,teexact,'--o')
grid on
title('circular fin temp'),xlabel('x(m)'),ylabel('T(C)')
legend('numerical solution','analytical solution')

Temperature Output
100.000
86.537
68.128
55.880
48.224
44.180
43.706
the exact solution temperature are
T1=64.22
T2=36.91
T3=26.50
T4=22.60
T5=21.30
in plot delta x is upto 1.4 this should be 1m

thank you