This code performs heat transfer.

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
Index exceeds the number of array elements. Index must not exceed 1.

3 Comments

Torsten
Torsten on 25 Apr 2023
Edited: Torsten on 25 Apr 2023
Sorry, but your code is full of errors. Maybe it makes more sense if you explain what you are trying to do.
This is part of the code to find how long it takes a slab of steele to cool from 1200° C to 200° C by using water jets at 40° C
Seems you try to solve the transient heat conduction equation in 1d.
MATLAB's "pdepe" is the correct tool to do this.

Sign in to comment.

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
ttot = 826.5942

Categories

Asked:

on 25 Apr 2023

Edited:

on 28 Apr 2023

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!