Help! I'm getting NaN + NaNi values...

clear all;
format long
AIRV=2.0;  %DRYING AIR VELOCITY, M/S
 
TAIR=45+273.15;  %TEMPERATURE OF AIR, C
 
RH=20/100;  %DRYING AIR RELATIVE HUMEDITY, %
 
PVSAT=100*exp(27.0214-(6887/(TAIR))-5.31*log((TAIR)/273.16));  %SATURATED VAPOR PRESSURE AT TAIR, PA
OMEGA=0.622*RH*PVSAT/(101000-RH*PVSAT);   %ABSOLUTE HUMIDITY, KG WATER/ KG DRY AIR
 
%DT=0.2; %STEP SIZE FOR TIME, SEC
RUNTIME=1*60*15; %1min * 60seconds * 15mins
STORETIMESTEP=30; %Time step for data storing, measure parameters every 30secs
DT=0.1; %STEP SIZE FOR TIME, SEC
TIME_LOOP=RUNTIME/DT; 
 
TH=(3)/1000;   %THICKNESS OF RODUCT, M
PL=60/1000;   %LENGTH OF RODUCT, M
PW=PL;    %WIDTHE OF RODUCT,M
NN=11;     %NUMBER OF INTERNAL GRID
DX=TH/(NN-1);  %GRID SIZE
SAREA=PL*PW;  %SURFACE AREA OF RODUCT, M^2
PVOLUME=TH*PL*PW;  %VOLUME OF RODUCT, M^3
 
INIMOIST=4.6; %INITIAL MOISTURE CONTENST, KG MOISTURE/KG DP
PMOIST=INIMOIST; %PRESENT MOISTURE CONTENST, KG MOISTURE/KG DP
RO=(11.354*PMOIST^4-92.25*PMOIST^3+277.04*PMOIST^2-402.65*PMOIST)+1110.60; %density OF Okara, kg/m^3
 
CP=(-44.93*PMOIST^4+382.13*PMOIST^3-1246.30*PMOIST^2+2123.60*PMOIST)+1678.10;     %Specific Heat OF Okara
IPWEIGHT=RO*PVOLUME;  %Initial Okara weight, kg
WTMOIST(1)=INIMOIST*IPWEIGHT/(1+INIMOIST);  %WEIGHT OF MOISTURE, KG
WTDP=IPWEIGHT-WTMOIST(1);  %WEIGHT OF DRY RODUCT, KG
 
ROS=1253;   %DENSITY OF B-DRY SOLID RODUCT, KG OF SOLID RODUCT/M^3 OF BD SOLID
ROSP=WTDP/PVOLUME;   %DENSITY OF SOLID IN RODUCT, KG OF SOLID/M^3 OF RODUCT
ROW=WTMOIST(1)/(PVOLUME-WTDP/ROS);   %DENSITY OF WATER, KG OF WATER/M^3 OF WATER
 
% CALCULATION OF PRESENT DENSITY OF SOLID IN RODUCT, KG OF SOLID/M^3 OF RODUCT
for I=1:NN
OROSP(I)=ROSP;  %INITIAL OR OLD DENSITY OF SOLID IN RODUCT, KG OF SOLID/M^3 OF RODUCT
end
% END OF CALCULATION OF PRESENT DENSITY OF SOLID IN RODUCT, KG OF SOLID/M^3 OF RODUCT
 
% CALCULATION OF MOISTURE CONTENT AT GRID POINTS
for I=1:NN
MOISTCONT(I)=WTMOIST(1)/WTDP;  %INITIAL MOISTURE CONTENT AT INNER GRID, KG MOISTURE/KG DP
end
% END OF CALCULATION OF MOISTURE CONTENT AT GRID POINTS
 
MOISTCONT1(1,1:length(MOISTCONT)) = MOISTCONT;
 
ROMP=WTMOIST(1)/PVOLUME;   %MASS CONCENTRATION OF MOISTURE IN RODUCT, KG OF MOISTURE/M^3 OF RODUCT 
 
T=(30+273.15)*ones(NN,1);   %INITIAL TEMPERATURE OF INNER GRIDS, C
PTT=T(NN);   %PREVIOUS TEMPERATURE OF TOP GRID, C
 
KN=(0.0126*PMOIST^3-0.0862*PMOIST^2+0.2018*PMOIST+0.1575);  %Thermal Conductivity of Okara, W/m.K
 
TD=-8*10^-10*PMOIST^4+7*10^-9*PMOIST^3-2*10^-8*PMOIST^2+4*10^-8*PMOIST+9*10^-8; %Thermal Diffusivity of Okara, m^2/s 
ROM=ROMP*ones(NN,1);  %LOCAL MASS CONCENTRATION OF MOISTURE, KG OF MOISTURE/M^3 OF RODUCT
 
% CALCULATION OF HEAT AND MASS TRANSFER COEFFICIENT IN AIR
ROAIR=-0.00000003510101*(TAIR-273.15)^3+0.00001583982684*(TAIR-273.15)^2-0.00469952020202*(TAIR-273.15)+1.29213571428571;  %DENSITY OF AIR (10<=TAIR<=80 C), KG/M^3
MEWAIR=(0.00000017676768*(TAIR-273.15)^3-0.00005541125541*(TAIR-273.15)^2+0.04983297258297*(TAIR-273.15)+17.1964285714285)*0.000001;  %ABSOLUTE VISCOSITY OF AIR (10<=TAIR<=80 C), KG/M S
NEWAIR=MEWAIR/ROAIR;  %KINEMATIC VISCOSITY OF AIR (10<=TAIR<=80 C), M^2/S
PRAIR=-0.00000002272727*(TAIR-273.15)^3+0.0000041991342*(TAIR-273.15)^2-0.00035335497835*(TAIR-273.15)+0.719;  %AIR PRANDTL NUMBER (10<=TAIR<=80 C) 
TCONAIR=(0.00000068181818*(TAIR-273.15)^3-0.0001474025974*(TAIR-273.15)^2+0.08029112554113*(TAIR-273.15)+24.0835714285714)*0.001;  %THERMAL CONDUCTIVITY OF AIR (10<=TAIR<=80 C), W/M C
 
REYMOND=AIRV*PL/NEWAIR; %REYNOLDS NUMBER OF DRYING AIR
HA=(TCONAIR/PL)*0.664*(REYMOND^0.5)*((ROAIR)^(1/3));   %CONVECTIVE HEAT TRANSFER COEFF. OF AIR, J/S/M^2/K
 
DIFAIR=0.000026; %DIFFUSION COEFFICIENT OF WATER VAPOR IN AIR, M^2/S
SCHMIDT=NEWAIR/DIFAIR;  %SCHMIDT NUMBER OF AIR
HM=(DIFAIR/PL)*0.664*(REYMOND^0.5)*((SCHMIDT)^(1/3));   %MASS TRANSFER COEFF. BET. SURFACE AND AIR, M/S
% END OF CALCULATION OF HEAT AND MASS TRANSFER COEFFICIENT IN AIR
 
% CALCULATION OF HEAT OF VAPORIZATION
PMOISTL=MOISTCONT(NN);
INITIAL_MOISTURE_CONTENT=INIMOIST;
%break
PRESENTTIME=0;
 
%%%%%%%%%%%%%%%%%%%%%%%%
%Calculation of average moisture content
NEWPRESENTTIME=0;
JJJ=1;
TOTALMOISTURE=0;
for I=1:NN
TOTALMOISTURE=TOTALMOISTURE+MOISTCONT(I);
end
AVMOISTCONT(JJJ)=TOTALMOISTURE/NN;
TIMESTEP(JJJ)=JJJ-1;
%%%%%%%%%%%%%%%%%%%
NEWPRESENTTIME=0;
JJJ=1;
TOTALTEMP=0;
for I=1:NN
TOTALTEMP=TOTALTEMP+T(I);
end
AVTEMP(JJJ)=TOTALTEMP/NN;
TIMESTEP(JJJ)=JJJ-1;
%%%%%%%%%%%%%%%%%%%
 
for LOOP1=1:TIME_LOOP
    
if PMOISTL < 0.2
    HFG=1000*(-2.394*(T(NN)-273.15)+2502.1+(-8.206767191142*PMOISTL^4+4.00032779720271*PMOISTL^3-0.61613964160838*PMOISTL^2+0.02368187645688*PMOISTL+0.00116308333333)*1000000); %LATENT HEAT OF VAPORIZATION, J/KJ
else
    HFG=1000*(-2.394*(T(NN)-273.15)+2502.1);
end
% END OF CALCULATION OF HEAT OF VAPORIZATION
 
for I=1:NN
DIFG(I)=(1.125*10^-6*MOISTCONT(I))+(2.11875*10^-10*T(I))-9.00924*10^-8;
end
%END OF CALCULATION OF MOISTURE DIFFUSIVITY AT DIFFERENT GRID POINT OF POTATO
 
MW=18.015;  %MOLECUALR WT. OF WATER, 
 
PMOISTL=MOISTCONT(NN);
 
AW=0.000007*PMOISTL^3-0.001*PMOISTL^2+0.0479*PMOISTL+0.2112 
 
RU=8315;    %UNIVERSAL GAS CONSTANT
 
PVS=AW*100*exp(27.0214-(6887/(T(NN)))-5.31*log((T(NN))/273.16));    %PARTIAL PRESSURE OF SAT. VAPOR AT SURFACE,PA 
PVAIR=RH*100*exp(27.0214-(6887/(TAIR))-5.31*log((TAIR)/273.16));    %PARTIAL PRESSURE OF SAT. VAPOR AT AIR,
JMS=HM*(MW/RU)*(PVS/(T(NN))-PVAIR/(TAIR));   %MOISTURE FLUX AT SURFACE, KG OF MOISTURE/M^2/SEC 
 
%CALCULATION OF ALPHA AT DIFFERENT GRID POINT OF POTATO
   ALPHAI(1)=-2*KN/(DX*DX*CP*RO);
   ALPHAI(2)=2*KN/(DX*DX*CP*RO);
 
   J=2;
   for I=2:NN-1
      ALPHAI(J+1)=KN/(DX*DX*CP*RO);
      ALPHAI(J+2)=-2*KN/(DX*DX*CP*RO);
      ALPHAI(J+3)=KN/(DX*DX*CP*RO);
      J=J+3;
   end
   ALPHAI(3*(NN-2)+2+1)=2*KN/(DX*DX*CP*RO);
   ALPHAI(3*(NN-2)+2+2)=-2*(HA+KN/DX)/(DX*CP*RO);
%END OF CALCULATION OF ALPHA AT DIFFERENT GRID POINT OF POTATO
   COEFFT=2*(HA*TAIR-HFG*JMS)/(DX*CP*RO);      
   COEFFC=0;
 
%CALCULATION OF BITA AT DIFFERENT GRID POINT OF POTATO
BITAI(1)=-2*DIFG(1)/(DX*DX);
BITAI(2)=2*DIFG(1)/(DX*DX);
J=2;
for I=2:NN-1
    BITAI(J+1)=DIFG(I)/(DX*DX); 
    BITAI(J+2)=-2*DIFG(I)/(DX*DX); 
    BITAI(J+3)=DIFG(I)/(DX*DX);
    J=J+3;
end
 
BITAI(3*(NN-2)+2+1)=2*DIFG(NN)/(DX*DX);       
BITAI(3*(NN-2)+2+2)=-2*DIFG(NN)/(DX*DX); 
%END OF CALCULATION OF BITA AT DIFFERENT GRID POINT OF POTATO
COEFFR=-2*JMS/DX;
 
MT=[ALPHAI(1) ALPHAI(2) 0 0 0 0 0 0 0 0 0;...
      ALPHAI(3) ALPHAI(4) ALPHAI(5) 0 0 0 0 0 0 0 0;...
      0 ALPHAI(6) ALPHAI(7) ALPHAI(8) 0 0 0 0 0 0 0;...
      0 0 ALPHAI(9) ALPHAI(10) ALPHAI(11) 0 0 0 0 0 0;...
      0 0 0 ALPHAI(12) ALPHAI(13) ALPHAI(14) 0 0 0 0 0;...
      0 0 0 0 ALPHAI(15) ALPHAI(16) ALPHAI(17) 0 0 0 0;...
      0 0 0 0 0 ALPHAI(18) ALPHAI(19) ALPHAI(20) 0 0 0;...
      0 0 0 0 0 0 ALPHAI(21) ALPHAI(22) ALPHAI(23) 0 0;...
      0 0 0 0 0 0 0 ALPHAI(24) ALPHAI(25) ALPHAI(26) 0;...
      0 0 0 0 0 0 0 0 ALPHAI(27) ALPHAI(28) ALPHAI(29);...
      0 0 0 0 0 0 0 0 0 ALPHAI(30) ALPHAI(31)];
 
ST=[COEFFC; 0; 0; 0; 0; 0; 0; 0; 0; 0; COEFFT];
 
MRO=[BITAI(1) BITAI(2) 0 0 0 0 0 0 0 0 0;...
      BITAI(3) BITAI(4) BITAI(5) 0 0 0 0 0 0 0 0;...
      0 BITAI(6) BITAI(7) BITAI(8) 0 0 0 0 0 0 0;...
      0 0 BITAI(9) BITAI(10) BITAI(11) 0 0 0 0 0 0;...
      0 0 0 BITAI(12) BITAI(13) BITAI(14) 0 0 0 0 0;...
      0 0 0 0 BITAI(15) BITAI(16) BITAI(17) 0 0 0 0;...
      0 0 0 0 0 BITAI(18) BITAI(19) BITAI(20) 0 0 0;...
      0 0 0 0 0 0 BITAI(21) BITAI(22) BITAI(23) 0 0;...
      0 0 0 0 0 0 0 BITAI(24) BITAI(25) BITAI(26) 0;...
      0 0 0 0 0 0 0 0 BITAI(27) BITAI(28) BITAI(29);...
      0 0 0 0 0 0 0 0 0 BITAI(30) BITAI(31)];
 
SRO=[0; 0; 0; 0; 0; 0; 0; 0; 0; 0; COEFFR];
      
      K1=DT*(MT*T+ST);
      K2=DT*(MT*(T+K1/2)+ST);
      K3=DT*(MT*(T+K2/2)+ST);
      K4=DT*(MT*(T+K3)+ST);
      
      T=T+(K1+2*K2+2*K3+K4)/6;
            
      L1=DT*(MRO*ROM+SRO);
      L2=DT*(MRO*(ROM+L1/2)+SRO);
      L3=DT*(MRO*(ROM+L2/2)+SRO);
      L4=DT*(MRO*(ROM+L3)+SRO);
      
      ROM=ROM+(L1+2*L2+2*L3+L4)/6;      
      % CALCULATION OF MOISTURE CONTENT AT GRID POINTS 
      
      for I=1:NN
MOISTCONT(I)=ROM(I)*DX*PL*PW/(WTDP/(NN-1));   %MOISTURE CONTENT AT INNER GRID, KG MOISTURE/KG DP
end
% END OF CALCULATION OF MOISTURE CONTENT AT GRID POINTS
PMOISTL=MOISTCONT(NN);
 
PRESENTTIME=PRESENTTIME+DT
 
%Calculation of average moisture content
NEWPRESENTTIME=NEWPRESENTTIME+DT;
 
if NEWPRESENTTIME >= STORETIMESTEP
JJJ=JJJ+1;
TOTALMOISTURE=0;
TOTALTEMP=0;
for I=1:NN
TOTALMOISTURE=TOTALMOISTURE+MOISTCONT(I);
TOTALTEMP=TOTALTEMP+T(I);
end
AVMOISTCONT(JJJ)=TOTALMOISTURE/NN;
AVTEMP(JJJ)=TOTALTEMP/NN;
TIMESTEP(JJJ)=JJJ-1;
NEWPRESENTTIME=0;
end
%%%%%%%%%%%%%%%%%%%
 
end      %End of LOOP1
TIMESTEP=TIMESTEP';
AVMOISTCONT=AVMOISTCONT';
AVTEMP=AVTEMP';
 
INITIAL_MOISTURE_CONTENT;
FINAL_MOISTURE_CONTENT=MOISTCONT
      T
      ROM

Note: I'm getting this NaN + NaNi value when I input the equation highlighted in bold (DIFG(I)=(1.125*10^-6*MOISTCONT(I))+(2.11875*10^-10*T(I))-9.00924*10^-8;).. Can anyone please advise, thanks!

[Edwin Yong]

Still seeking for help to this problem if anyone can help to resolve..

[Rik]

Have a read here and here. It will greatly improve your chances of getting an answer.

I will now edit your question.

clear all should never be used in a script like this. Use clear or clearvars instead.

Answer 1

[KALYAN ACHARJYA]

Till TIME_LOOP=5, its perfectly works. The code is so lengthy, it quite difficult to bebug. Please check carefully. The main issue to debug the code with the following lines

MOISTCONT(I)=ROM(I)*DX*PL*PW/(WTDP/(NN-1));

The variable is changing twice within the same loop.