This code performs heat transfer.
Show older comments
Tf=100; % Final Temp
Tint=1000; % Initial Temp
Tinf=20; % Water Temp
h=3000; % W/m^2*K
rho=5550; %kg/m^3
dx=0.005; %m
L=.25/2;%m
num=(L/dx);
T=zeros(num,1);
T(:,1)=Tint;
Cp=600;
k=27.5;
Alpha=(k/(rho*Cp));
Bi=((h*dx)/k);
dt=(dx^2)/(2*Alpha*(1+Bi));
Fo=((Alpha*dt)/(dx^2));
i=1;
while T(num,i)>Tf
for j=1:num
if j==1
T(j,i+1)=(2*Fo(j)*Bi(j))*Tinf+(2*Fo(j))*T(j+1,i)+(1-(2*Fo(j)*Bi(j))-(2*Fo(j)))*T(j,i);
elseif j==num
T(j,i+1)=2*Fo(j)*T((j-1),i)+(1-(2*Fo(j)))*T(j,i);
else
T(j,i+1)=(Fo(j))*(T((j-1),i)+T((j+1),i))+(1-(2*Fo(j)))*T(j,i);
end
if (20<=(T(j,i)))&&((T(j,i))<600)
Cp(j)=425+(.773)*(T(j,i))-.00169*((T(j,i))^2)+.00000222*((T(j,i))^3);
elseif (600<=(T(j,i)))&&(T(j,i))<(735)
Cp(j)=666+(13002/(738-(T(j,i))));
elseif (735<=(T(j,i)))&&(T(j,i))<(900)
Cp(j)=545+(17820/((T(j,i))-731));
elseif (900<=(T(j,i)))&&(T(j,i))<=(1000)
Cp(j)=600;
end
if(20<=(T(j,i)))&&((T(j,i))<800)
k(j)=54-.0333*(T(j,i));
elseif (800<=(T(j,i)))&&(T(j,i))<=(1000)
k(j)=27.3;
end
Alpha=(k(j)/(rho*Cp(j)));
Bi(j)=((h*dx)/k(j));
Fo(j)=((Alpha*dt)/(dx^2));
end
i=i+1;
end
Answers (1)
Your code is written expecting Fo and Bi to be a vectors, but they are scalars, so only ever has 1 element. You get this error any time you try to index using a value greater than 1.
The easiest solution based on the information we have is to not index these 2 variables.
Tf=200; % Final Temp
Tint=1200; % Initial Temp
Tinf=40; % Water Temp
h=4000; % W/m^2*K
rho=7850; %kg/m^3
dx=0.005; %m
L=.5/2;%m
num=(L/dx)+1;
T=zeros(num,1);
T(:,1)=Tint;
Cp=650; % initial Cp
k=27.3; % Initial k
% used to solve for dt
Alpha=(k/(rho*Cp));
Bi=((h*dx)/k);
dt=(dx^2)/(2*Alpha*(1+Bi));
Fo=((Alpha*dt)/(dx^2));
i=1;
while T(num,i)>Tf
for j=1:num
if j==1
T(j,i+1)=(2*Fo*Bi)*Tinf+(2*Fo)*T(j+1,i)+(1-(2*Fo*Bi)-(2*Fo))*T(j,i);
elseif j==num
T(j,i+1)=2*Fo*T((j-1),i)+(1-(2*Fo))*T(j,i);
else
T(j,i+1)=(Fo)*(T((j-1),i)+T((j+1),i))+(1-(2*Fo))*T(j,i);
end
% Cp and k need to be defined each time it finds a temperature since they are both temperature dependent
if (20<=(T(j,i)))&&((T(j,i))<600)
Cp(j)=425+(.773)*(T(j,i))-.00169*((T(j,i))^2)+.00000222*((T(j,i))^3);
elseif (600<=(T(j,i)))&&(T(j,i))<(735)
Cp(j)=666+(13002/(738-(T(j,i))));
elseif (735<=(T(j,i)))&&(T(j,i))<(900)
Cp(j)=545+(17820/((T(j,i))-731));
elseif (900<=(T(j,i)))&&(T(j,i))<=(1200)
Cp(j)=650;
end
if(20<=(T(j,i)))&&((T(j,i))<800)
k(j)=54-.0333*(T(j,i));
elseif (800<=(T(j,i)))&&(T(j,i))<=(1200)
k(j)=27.3;
end
%These values are recalculated each time as well based on the Cp and k
Alpha=(k(j)/(rho*Cp(j)));
Bi=((h*dx)/k(j));
Fo=((Alpha*dt)/(dx^2));
end
i=i+1;
end
ttot=(i-1)*dt
Categories
Find more on Thermal Analysis in Help Center and File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!