ode45 returns NaN

%% Clean the Workspace
clc
clear all

%% Variables
global mS Ts_in
Vceo = 102.687;                 % m3
Ts_in = 104.8;                  % C, inlet steam temp
mS = 85.3;                      % kg/sec
T0 = 57.2;                      % Initial temp
t_interval = [0.1 20];          % Time interval
IC = T0;                        % Initial Condition

%% Solve the DE
[t,Tsol] = ode45(@(t,T) cfunc(t,T) , t_interval , IC);

%% From Tsat to Psat
for i=1:size(Tsol)
Psol(i) = XSteam('psat_T',Tsol(i));
i = i + 1;
end

%% Plotting
figure;
subplot(1,2,1)
plot(t,Tsol)
xlabel('Time (s)')
ylabel('Temperature (C)')
subplot(1,2,2)
plot(t,Psol)
xlabel('Time (s)')
ylabel('Pressure (bar)')

function dTdt = cfunc(t,T)

global mS Ts_in

A = 3300;                       % m2, heat transfer area
U = 9.9616;                     % heat transfer coefficient
T_hi = 104.8;
T_ho = T;
T_ci = 52.2;
T_co = 97.8;
LMTD = ((T_hi-T_co)-(T_ho-T_ci))/(log((T_hi-T_co)/(T_ho-T_ci)));
Q = LMTD*U*A;
deltaHcond = XSteam('hL_T',T) - XSteam('hV_T',T);
cvS = XSteam('CvV_T',T);
cpW = XSteam('CpL_T',T);
mCond = -Q/deltaHcond;

    dTdt =  ((Ts_in*XSteam('CpV_T',Ts_in)*mS-mCond*deltaHcond-mCond*T*cpW) / (cvS*(mS-mCond)) - T)/t ;
    
end

[Torsten]

Check which variable makes dTdt become NaN.

[Les Beckham]

You might want to try entering this command in the command window before you run the code. The debugger will stop when a NaN is detected and then you can examine the various variables to see what is happening.

dbstop if naninf

Answer 1

[Luca Ferro]

The values of t and T from iteration 1 to 5 are:

t= 0.1 T=57.2

t= 0.1 T=59.4

t= 0.1 T=57.49

t= 0.1 T=377.82

t= 0.1 T=NaN

The NaN is cause by the line

deltaHcond = XSteam('hL_T',T) - XSteam('hV_T',T);

and then propagates to everything that follows it.

with that call you go to the XSteam file lines:

case 'hl_t'

T = toSIunit_T(In1); %%%%%%%%%%%%%%%%%%%%% THIS (1)

if T > 273.15 && T < 647.096

p = p4_T(T);

Out = fromSIunit_h(h4L_p(p));

else

Out = NaN;

end

and

case 'hv_t'

T = toSIunit_T(In1); %%%%%%%%%%%%%%%%%%%%% THIS (2)

if T > 273.15 && T < 647.096

p = p4_T(T);

Out = fromSIunit_h(h4V_p(p));

else

Out = NaN;

end

The marked line (1) converts 377.82 to 650.97 determining the validation of the else condition NaN. Same happens on line (2).

the toSIunit function makes sense since it just converts from celsius to kelvin by adding 273.15.

The result is coherent with the conditions you have written

[Talha]

Thank you.