MATLAB: NaN problem with heat transfer program

Hi I am working on a program to calculate the heat lost by a pipe while water is traveling through it.

Here is the code

clear
WaterTemperatureCelcius = [0.01 10:10:100];
WaterTemperature = WaterTemperatureCelcius+273.15;
WaterDensityMatrix= [916.8 999.8 998.3 995.7 992.3 988 983 978 972 965 958];
WaterDynamicViscosityMatrix = [1.78 1.31 1.00 0.798 0.653 0.547 0.467 0.404 0.355 0.314 0.281]; 
WaterKinematicViscosityMatrix = [1.792 1.304 1.004 0.801 0.658 0.553 0.474 0.413 0.365 0.326 0.295]; 
KWaterMatrix = [-0.07 0.088 0.207 0.303 0.385 0.457 0.523 0.585 0.643 0.665 0.752];
PrWaterMatrix = [13.67 9.47 7.01 5.43 4.34 3.56 2.99 2.56 2.23 1.96 1.75];
CWaterMatrix = [4.21 4.193 4.183 4.179 4.179 4.182 4.185 4.191 4.198 4.208 4.219]; 

AirTemperature = 100:50:500;
AirDynamicViscosityMatrix = [0.6924 1.0283 1.3289 1.488 1.893 2.075 2.286 2.484 2.671]; 
AirKinematicViscosityMatrix = [1.923 4.343 4.49 9.49 15.68 20.76 25.9 28.86 37.9];
AirDensityMatrix = [3.601 2.3675 1.7684 1.4128 1.774 0.998 0.8826 0.7833 0.7048];
KAirMatrix = [0.009246 0.013735 0.01809 0.02227 0.02624 0.03003 0.03365 0.03707 0.04038]; 
PrAirMatrix = [0.770 0.753 0.739 0.722 0.708 0.697 0.689 0.683 0.680];

PipeTemperatureMatrixCelcius = [-100 0 100 200 300 400 600 800 1000 1200];
PipeTemperatureMatrix = PipeTemperatureMatrixCelcius+273.15;
KPipeMatrix = [87 73 67 62 55 48 40 36 35 36];

Two = input ('Guess the outside temperature of the pipe wall (K): ');
Twi = input ('Guess the inside temperature of the pipe wall (K): ');
Tm = input ('Guess the temperature in the middle of the flow (K): ');
Do = input ('Please enter the outside diameter of the pipe (m): ');
Di = input ('Please enter the inside diameter of the pupe (m): ');
Flow = input ('Please enter 1 for laminar flow and 2 for turbulent flow: ');

Tinf = 266.48; 
Vinf = 6.7056; 
Tin = 372.04; 
L = 30.48; 
dx = 0.1; 
ivalue = L/dx;
Pi = 3.141593;
TTest1 = 1;
TTest2 = 1;
TTest3 = 1;

Tx = Tin;
Ai = Pi*Di*dx;
Ao = Pi*Do*dx;

for i=1:ivalue,
    if i == 1
        Tx=Tin;
    end
    if i > 1
         Tx = Tout;
    end
    while (abs(TTest1-Two) > 0.001) && (abs(TTest2-Twi) > 0.001) && (abs(TTest3-Tm) > 0.001)
        disp(['Two is ' num2str(Two)])
        disp(['Twi is ' num2str(Twi)])
        disp(['Tm is ' num2str(Tm)])
        Tfo = (Two+Tinf)/2;
        Tfi = (Twi+Tin)/2;
        PipeTemperature = Two-Twi;

          WaterKinematicViscosity = (interp1(WaterTemperature, WaterKinematicViscosityMatrix, Tfi))*(10^(-6)); 
          WaterDynamicViscosity = (interp1(WaterTemperature, WaterDynamicViscosityMatrix, Tfi))*(10^(-2)); 
          WaterDensity = interp1(WaterTemperature, WaterDensityMatrix, Tfi); 
          KWater = interp1(WaterTemperature, KWaterMatrix, Tfi); 
          PrWater = interp1(WaterTemperature, PrWaterMatrix, Tfi);
          CWater = (interp1(WaterTemperature, CWaterMatrix, Tfi)*1000);

          AirKinematicViscosity = (interp1(AirTemperature, AirKinematicViscosityMatrix, Tfo))*(10^(-6)); 
          AirDynamicViscosity = (interp1(AirTemperature, AirDynamicViscosityMatrix, Tfo))*(10^(-5)); 
          AirDensity = interp1(AirTemperature, AirDensityMatrix, Tfo);
          KAir = interp1(AirTemperature, KAirMatrix, Tfo);
          PrAir = interp1(AirTemperature, PrAirMatrix, Tfo);

          KPipe = interp1(PipeTemperatureMatrix, KPipeMatrix, PipeTemperature);

          if Flow==1
              ReDi = 2000;
              hi = (4.363*KWater/Di);
          else 
              ReDi = 5000;
              hi= (KWater*0.023*(ReDi^0.8)*(PrWater^0.3));
          end

          Vm = (ReDi*WaterKinematicViscosity)/Di;
          ReDf = Vinf*Do/AirKinematicViscosity;
          disp(['ReDf is ' num2str(ReDf)])
          disp(['ReDi is ' num2str(ReDi)])
          disp(['hi is ' num2str(hi)])
          if 0.4 < ReDf <= 4
              c=0.989;
              n=0.33;
          elseif 4 < ReDf <= 40
              c=0.911;
              n=0.385;
          elseif 40 < ReDf <= 4000
              c=0.683;
              n=0.466;
          elseif 4000 < ReDf <= 40000
              c=0.193;
              n=0.618;
          elseif 40000 < ReDf <= 400000
              c=0.0266;
              n=0.805;
          end
          disp(['c is ' num2str(c)])
          disp(['n is ' num2str(n)])
          NuDf = c*(ReDf^n)*(PrAir^(1/3));
          disp(['NuDf is ' num2str(NuDf)])
          ho = NuDf*KAir/Do;
          MassFlowRate = WaterDensity*Vm*Pi*Di*Di/4;
          disp(['Mass flow rate is ' num2str(MassFlowRate) ' kg/s'])

          q = (Tx-Tinf)/(((1/(ho*Ao))+((1/(2*Pi*KPipe*dx))*log(Do/Di))+(1/(hi*Ai))+(1/(2*MassFlowRate*CWater))));
          TwoNew = ((q/(ho*Ao))+Tinf);  
          TwiNew = (q*(log(Do/Di))/(2*Pi*KPipe*dx))+TwoNew;
          TmNew = Tx-(q/(2*MassFlowRate*CWater));
          disp(['TwoNew is ' num2str(TwoNew)])
          disp(['TwiNew is ' num2str(TwiNew)])
          disp(['TmNew is ' num2str(TmNew)])
          TTest1 = Two;
          TTest2 = Twi;
          TTest3 = Tm;
          Two = TwoNew;
          Twi = TwiNew;
          Tm = TmNew;

       end
          Txdx = (2*TmNew) - Tx;
          Tout = Tx + Txdx;
  end

q1 = ho*(Twi-Tinf)*Ao;
q2 = (Twi-Two)/((1/(2*Pi*KPipe*dx))*log(Do/Di));
q3 = hi*(Tm-Twi)*Ai;
q4 = MassFlowRate*CWater*(Tx-((2*Tm)-Tx));

disp (['The value of q1 is ' num2str(q1) ' W'])
disp (['The value of q2 is ' num2str(q2) ' W'])
disp (['The value of q3 is ' num2str(q3) ' W'])
disp (['The value of q4 is ' num2str(q4) ' W'])
disp (['The value of q5 is ' num2str(q) ' W'])

I keep getting NaN output for the temperature, which stops the loop

    Guess the outside temperature of the pipe wall (K): 320
    Guess the inside temperature of the pipe wall (K): 330
    Guess the temperature in the middle of the flow (K): 340
    Please enter the outside diameter of the pipe (m): 1
    Please enter the inside diameter of the pupe (m): 0.9
    Please enter 1 for laminar flow and 2 for turbulent flow: 1
    Two is 320
    Twi is 330
    Tm is 340
    ReDf is 451765.1016
    ReDi is 2000
    hi is 3.0572
    c is 0.989
    n is 0.33
    NuDf is 64.8209
    Mass flow rate is 0.51629 kg/s
    TwoNew is NaN
    TwiNew is NaN
    TmNew is NaN
    The value of q1 is NaN W
    The value of q2 is NaN W
    The value of q3 is NaN W
    The value of q4 is NaN W
    The value of q5 is NaN W

Any advice on solving this would help.

Thanks.

Answer 1

[Walter Roberson]

Your PipeTemperature is 320 - 330 = -10, which is outside of the range of your PipeTemperatureMatrix, so the interp1() returns nan for the interpolation yielding KPipe.

[Hemdeep]

Thank you!

Answer 2

[Andrew Newell]

Try running this code after the command

dbstop if naninf

to see where the problem is occurring.

[Hemdeep]

Thanks for the feedback. This technique will be very helpful for my future codes.

Answer 3

[Hung Pham]

good document