MATLAB: Not a Number (NaN) values returned from a loop

forheat transfernanrunge kutta

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

Best Answer

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
Related Question