Not a Number (NaN) values returned from a loop

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

[madhan ravi]

dbstop if naninf

Answer 1

[Sriram Tadavarty]

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

[Thomas Faulkner]

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