Not a Number (NaN) values returned from a loop

3 次查看（过去 30 天）
Thomas Faulkner
Thomas Faulkner 2020-3-20
评论： Thomas Faulkner 2020-3-20
Hi everyone,
I am having trouble with the values being returned in my loops.
my k, w and f terms are stored in the workspace as NaN after i run my code.
I appreciate there may be a lot wrong with my loop and the terms within it but any help would be great.
Thanks.
clc
clear
% Parameters & Boundary Conditions
dTdz=30.12; Ri=0.05; Ro=0.075; Rm=0.0621; hcoeff=100; To=300; Ti=0; Zo=0; Zi=1.3*Ri*hcoeff*(Ti-To); uRi=0; uRo=0; tawi=4.3465; tawo=3.792; npoints=1000;
h=(Ro-Ri)/npoints;
r=(Ri+h:h:Ro)';
uR= zeros(npoints,1);
T= zeros(npoints,1);
Z = zeros(npoints,1);
%
% Call Runge Kutta Function Fourth Order
for i=1:(npoints-1)
T(1)=Ti;
Z(1)=Zi;
uR(1)= uRi;
if r(i)<Rm
k1=h*dUdrI(r(i),tawi,mu(T(i)),rho(T(i)),epsilon_m_(r(i),T(i),Rm,Ri,Ro,tawi,tawo),Rm,Ri);
w1=h*dTdr(Z(i),r(i),k(T(i)),rho(T(i)),Cp(T(i)),epsilon_m_(r(i),T(i),Rm,Ri,Ro,tawi,tawo));
f1=h*dZdr(r(i),rho(T(i)),Cp(T(i)),uR(i),dTdz);
k2=h*dUdrI((r(i)+h/2),tawi,mu(T(i)),rho(T(i)),epsilon_m_(r(i),T(i),Rm,Ri,Ro,tawi,tawo),Rm,Ri);
w2=h*dTdr(Z(i),(r(i)+h/2),k(T(i)),rho(T(i)),Cp(T(i)),epsilon_m_(r(i),T(i),Rm,Ri,Ro,tawi,tawo));
f2=h*dZdr((r(i)+h/2),rho(T(i)),Cp(T(i)),uR(i),dTdz);
k3=h*dUdrI((r(i)+h/2),tawi,mu(T(i)),rho(T(i)),epsilon_m_(r(i),T(i),Rm,Ri,Ro,tawi,tawo),Rm,Ri);
w3=h*dTdr(Z(i),(r(i)+h/2),k(T(i)),rho(T(i)),Cp(T(i)),epsilon_m_(r(i),T(i),Rm,Ri,Ro,tawi,tawo));
f3=h*dZdr((r(i)+h/2),rho(T(i)),Cp(T(i)),uR(i),dTdz);
k4=h*dUdrI((r(i)+h),tawi,mu(T(i)),rho(T(i)),epsilon_m_(r(i),T(i),Rm,Ri,Ro,tawi,tawo),Rm,Ri);
w4=h*dTdr(Z(i),(r(i)+h),k(T(i)),rho(T(i)),Cp(T(i)),epsilon_m_(r(i),T(i),Rm,Ri,Ro,tawi,tawo));
f4=h*dZdr((r(i)+h),rho(T(i)),Cp(T(i)),uR(i),dTdz);
uR(i+1) = uR(i) + (k1+2*k2+2*k3+k4)/6;
T(i+1) = T(i) + (w1+2*w2+2*w3+w4)/6;
Z(i+1) = Z(i) + (f1+2*f2+2*f3+f4)/6;
end
end
for i=(npoints-1):h:1
T(1)=To;
Z(1)=Zo;
uR(1)=uRo;
if r(i)>=Rm
k1=h*dUdrO(r(i),tawo,mu(T(i)),rho(T(i)),epsilon_m_(r(i),T(i),Rm,Ri,Ro,tawi,tawo),Rm,Ro);
w1=h*dTdr(Z(i),r(i),k(T(i)),rho(T(i)),Cp(T(i)),epsilon_m_(r(i),T(i),Rm,Ri,Ro,tawi,tawo));
f1=h*dZdr(r(i),rho(T(i)),Cp(T(i)),uR(i),dTdz);
k2=h*dUdrO((r(i)-h/2),tawo,mu(T(i)),rho(T(i)),epsilon_m_(r(i),T(i),Rm,Ri,Ro,tawi,tawo),Rm,Ro);
w2=h*dTdr(Z(i),(r(i)-h/2),k(T(i)),rho(T(i)),Cp(T(i)),epsilon_m_(r(i),T(i),Rm,Ri,Ro,tawi,tawo));
f2=h*dZdr((r(i)-h/2),rho(T(i)),Cp(T(i)),uR(i),dTdz);
k3=h*dUdrO((r(i)-h/2),tawo,mu(T(i)),rho(T(i)),epsilon_m_(r(i),T(i),Rm,Ri,Ro,tawi,tawo),Rm,Ro);
w3=h*dTdr(Z(i),(r(i)-h/2),k(T(i)),rho(T(i)),Cp(T(i)),epsilon_m_(r(i),T(i),Rm,Ri,Ro,tawi,tawo));
f3=h*dZdr((r(i)-h/2),rho(T(i)),Cp(T(i)),uR(i),dTdz);
k4=h*dUdrO((r(i)-h),tawo,mu(T(i)),rho(T(i)),epsilon_m_(r(i),T(i),Rm,Ri,Ro,tawi,tawo),Rm,Ro);
w4=h*dTdr(Z(i),(r(i)-h),k(T(i)),rho(T(i)),Cp(T(i)),epsilon_m_(r(i),T(i),Rm,Ri,Ro,tawi,tawo));
f4=h*dZdr((r(i)-h),rho(T(i)),Cp(T(i)),uR(i),dTdz);
uR(i+1)=uR(i)+(k1+2*k2+2*k3+k4)/6;
T(i+1)=T(i)+(w1+2*w2+2*w3+w4)/6;
Z(i+1)=Z(i)+(f1+2*f2+2*f3+f4)/6;
end
end
figure, plot (r,uR)
xlabel('Radius')
ylabel('Velocity')
title('Velocity vs Radius')
figure, plot (r,T)
xlabel('Radius')
ylabel('Temperature')
title('Temperature vs Radius')
figure, plot (r,Z)
xlabel('Radius')
ylabel('Heat transfer per unit length')
title('Heat transfer per unit length vs Radius')
% Heat Transfer per Unit Length
function [dZdr_] = dZdr (r,rho,Cp,uR,dTdz)
dZdr_ = -r*rho*Cp*uR*dTdz;
end
% Temperature Profile
function [dTdr_] = dTdr (Z,r,k,rho,Cp,epsilon_m_)
dTdr_ = Z/(r*(k+rho*Cp*epsilon_m_));
end
% Velocity Profile
function [dUdrI_] = dUdrI (tawi,mu,rho,epsilon_m_,Rm,r,Ri)
dUdrI_= tawi./(mu+rho*epsilon_m_)*((Rm^2-r^2)/(Rm^2-Ri^2))*(Ri/r);
end
function [dUdrO_] = dUdrO (tawo,mu,rho,epsilon_m_,r,Rm,Ro)
dUdrO_= tawo./(mu+rho*epsilon_m_)*((r^2-Rm^2)/(Ro^2-Rm^2))*(Ro/r);
end
%eddy diffusivities of momentum
function [epsilon_m] = epsilon_m_ (r,T,Rm,Ri,Ro,tawi,tawo)
if r<Rm
radratio=(r-Ri)/(Rm-Ri);
rEpsI=6.516e-4+3.9225e-1*radratio+(-6.0949e-1)*(radratio).^2+(2.3391e-1)*(radratio).^3+(1.410e-1)*(radratio).^4+(-9.6098e-2)*(radratio).^5;
epsilon_m=rEpsI*sqrt(tawi./rho(T))*(Rm-Ri);
else
radratio=(Ro-r)/(Ro-Rm);
rEpsO=6.516e-4+3.9225e-1*radratio+(-6.0949e-1)*(radratio).^2+(2.3391e-1)*(radratio).^3+(1.410e-1)*(radratio).^4+(-9.6098e-2)*(radratio).^5;
epsilon_m=rEpsO*sqrt(tawo./rho(T))*(Ro-Rm);
end
end
function [Cp_]=Cp(T)
Cp_=1010.4755 + 0.1151*(T) + 4.00e-5*(T).^2;
end
function [k_]=k(T)
k_=0.0243+(6.548e-5)*(T) - (1.65e-8)*(T).^2;
end
%Dynamic Viscosity
function [mu_]=mu(T)
mu_=1.747e-5 + 4.404e-8*(T) - 1.645e-11*((T).^2);
end
function [Pr_]=Pr(T)
Pr_=0.7057*10^(2.06e-5*(T));
end
function [rho_]=rho (T)
rho_ =1e5/(287*(T));
end
function [vis_]=vis(T)
vis_=1.380e-5 + (8.955e-8)*(T) - (1.018e-10)*(T)^2;
end

请先登录，再回答此问题。

采纳的回答

Sriram Tadavarty
Sriram Tadavarty 2020-3-20
Hi Thomas,
Update Ti with a value other than 0 and then this would provide some values other than Nan. When Ti is 0, the value of rho is calculated as 1e5/(287*0) and therefore set to Nan.
Hope this helps.
Regards,
Sriram
  1 个评论
Thomas Faulkner
Thomas Faulkner 2020-3-20
Hi Sriram,
Thanks for the help!
This has solved the issue for my NaN values however i am now getting complex numbers for these values.
Also, my second loop is not changing value and is stuck at zero.
What would you suggest I change here?
Kind regards,
Thomas

请先登录，再进行评论。

更多回答（0 个）

类别

在 Help Center 和 File Exchange 中查找有关 Mathematics 的更多信息

标签

Community Treasure Hunt

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

Start Hunting!

Translated by