ode45 returns NaN

1 次查看（过去 30 天）
Talha
Talha 2023-2-13
评论： Talha 2023-2-13
Ran in:
%% 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
  2 个评论
Torsten
Torsten 2023-2-13
编辑：Torsten 2023-2-13
Check which variable makes dTdt become NaN.
Les Beckham
Les Beckham 2023-2-13
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

请先登录，再进行评论。

请先登录，再回答此问题。

采纳的回答

Luca Ferro
Luca Ferro 2023-2-13
编辑：Luca Ferro 2023-2-13
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

更多回答（0 个）

类别

在 Help Center 和 File Exchange 中查找有关 Ordinary Differential Equations 的更多信息

标签

产品


版本

R2022b

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by