Must return a column vector error when using "ode45"

Question

Talha 2024-8-20

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/2146499-must-return-a-column-vector-error-when-using-ode45

评论： Talha 2024-8-20

XSteam.m

%% Clean the Workspace
clc
clear all
%% Format
format long
%% Variables
global m_ExtraSteam U0_shell V rho_0 U_steam
T0_shell = 105;                                       % initial shell side temp, degC
P0_shell = 1.2;                                  % initial pressure, bara
U0_shell = XSteam('uV_p', P0_shell);
Plimit  = 3.3;                                        % max pressure, bara
MW      = 0.018015;                             % kg/mol
R       = 0.00008205736;                        % bar.m3/mol.K
V       = 233.21;                                   % m3
rho_0   = XSteam('rhoV_p', P0_shell); 
tb      = 20;                                  % Upper time interval limit
m_ExtraSteam  = 35.08;           % Flowrate of extra steam due to ST trip, kg/sec
t_interval = [0.001 tb];                                     % Time interval
%% Solve the DE
[tsol,Usol] = ode45(@(t, U_steam) cfunc(t, U_steam) , t_interval , U0_shell);      % Temperature with time, degC
Warning: The value of local variables may have been changed to match the globals.  Future versions of MATLAB will require that you declare a variable to be global before you use that variable.
Warning: Trust-region-dogleg algorithm of FSOLVE cannot handle non-square systems; using Levenberg-Marquardt algorithm instead.
Equation solved at initial point.

fsolve completed because the vector of function values at the initial
point is near zero as measured by the value of the function tolerance,
and the problem appears regular as measured by the gradient.
Error using odearguments (line 94)
@(T,U_STEAM)CFUNC(T,U_STEAM) must return a column vector.

Error in ode45 (line 104)
    odearguments(odeIsFuncHandle,odeTreatAsMFile, solver_name, ode, tspan, y0, options, varargin);
%% Dynamic Parameter Calculation
Ptube_in = 7;
Ptube_out = 6.5;
m_tube = [707.9 707.9 707.9 707.9];
T_tube_in = [91.7 91.7 91.7 91.7];
    for m=1:lenght(Usol)
        U_steam = Usol(m);
        Tsol(m) = fsolve(@TfinderwU,100);
    end
    for p=5:length(Tsol)
T_tube_in(p) = 80.1;
m_tube(p) = 1131.2;
    end
T_tube_out = 98;
        for k = 1:length(Tsol)
            deltaH_tube(k) = XSteam('h_pT', Ptube_out, T_tube_out) - XSteam('h_pT', Ptube_in, T_tube_in(k));        %  Specific enthalpy difference between the inlet and outlet of the tube side, kJ/kg
   
            if ((m_tube(k)*deltaH_tube(k))/(XSteam('hV_T', Tsol(k)) - XSteam('hL_T', Tsol(k)))) >= m_ExtraSteam
                m_ExtraSteam_Condensation(k) = m_ExtraSteam;
            else
                m_ExtraSteam_Condensation(k) = (m_tube(k)*deltaH_tube(k))/(XSteam('hV_T', Tsol(k)) - XSteam('hL_T', Tsol(k)));
            end
        end
m_steam_accumulated = (m_ExtraSteam - m_ExtraSteam_Condensation)';                 % Steam accumulated in the shell side, kg/sec
nfor(1) = m_steam_accumulated(1)*(tsol(1))/MW;
        for k1 = 2:length(Tsol)
            nfor(k1) =  m_steam_accumulated(k1)*(tsol(k1)-tsol(k1-1))/MW;
        end
%% Final Pressure Calculation
%% Plotting
figure('Name','Temperature, Pressure vs Time','NumberTitle','off');
yyaxis left % subplot(1,2,1)
plot(tsol,Tsol)
xlabel('Time (s)')
ylabel('Temperature (C)')
yyaxis right % subplot(1,2,2)
plot(tsol,Psol)
xlabel('Time (s)')
ylabel('Pressure (bara)')
%% Display Results
fprintf('Max pressure during the observed time interval = %.3f bara.\n', max(Psol))
%% Functions
    % Main Function
function dUdt = cfunc(t,U_steam)
% Variables
global m_ExtraSteam U0_shell V rho_0 U_steam
Ptube_in = 7;
Ptube_out = 6.5;
T_tube_out = 98;
if t<3.3
m_tube = 707.9;
T_tube_in = 91.7;
else 
m_tube = 1131.2;
T_tube_in = 80.1;
end
% Find the Temp with Internal Energy 
T = fsolve(@TfinderwU,100);
deltaH_tube = XSteam('h_pT', Ptube_out, T_tube_out) - XSteam('h_pT', Ptube_in, T_tube_in);        %  Specific enthalpy difference between the inlet and outlet of the tube side, kJ/kg
if ((m_tube*deltaH_tube)/(XSteam('hV_T', T) - XSteam('hL_T', T))) >= m_ExtraSteam
        m_ExtraSteam_Condensation = m_ExtraSteam;
else
   m_ExtraSteam_Condensation = (m_tube*deltaH_tube)/(XSteam('hV_T', T) - XSteam('hL_T', T)); 
end
m_steam_accumulated = m_ExtraSteam - m_ExtraSteam_Condensation;                 % Steam accumulated in the shell side, kg/sec
steam_mass = rho_0*V + m_steam_accumulated*t;
U_condensate = XSteam('uL_T', T);
CpL = XSteam('CpL_T', T);
Qcv = m_tube*CpL*(T_tube_out-T_tube_in);
% Differential Equation
dUdt = (-Qcv + m_ExtraSteam*U0_shell - m_ExtraSteam_Condensation*U_condensate - U_steam*m_steam_accumulated)/steam_mass;
end
    % Temperature Finder Function with Internal Energy 
function Obj = TfinderwU(T)
global U_steam
uV_T = XSteam('uV_T', T);
Obj = U_steam- uV_T;
end

2 个评论
显示无隐藏无

Stephen23 2024-8-20

编辑：Stephen23 2024-8-20

在 MATLAB Online 中打开

XSteam.m

dUdt is empty when the error occurs. Possible causes (that should be adressed anyway):

Using GLOBAL (with default = empty) instead of parameterizing the function calls: https://www.mathworks.com/help/matlab/math/parameterizing-functions.html. If a variable is used in multiple functions then use nested functions.
using matrix operations everywhere, a complete lack of array operations: https://www.mathworks.com/help/matlab/matlab_prog/array-vs-matrix-operations.html. This makes it unlikely that the function CFUNC was written to fulfill the requirements given in the ODE45 documentation: "The function dydt = odefun(t,y), for a scalar t and a column vector y, must return a column vector..." source: https://www.mathworks.com/help/matlab/ref/ode45.html#bu00_4l_sep_shared-odefun.

Lets test it right now:

global m_ExtraSteam U0_shell V rho_0 U_steam
T0_shell = 105;                                       % initial shell side temp, degC
P0_shell = 1.2;                                  % initial pressure, bara
U0_shell = XSteam('uV_p', P0_shell);
Plimit  = 3.3;                                        % max pressure, bara
MW      = 0.018015;                             % kg/mol
R       = 0.00008205736;                        % bar.m3/mol.K
V       = 233.21;                                   % m3
rho_0   = XSteam('rhoV_p', P0_shell); 
tb      = 20;                                  % Upper time interval limit
m_ExtraSteam  = 35.08;           % Flowrate of extra steam due to ST trip, kg/sec
cfunc(1,[1,2]) % fails the requirements
Warning: The value of local variables may have been changed to match the globals.  Future versions of MATLAB will require that you declare a variable to be global before you use that variable.
Warning: Trust-region-dogleg algorithm of FSOLVE cannot handle non-square systems; using Levenberg-Marquardt algorithm instead.
Equation solved at initial point.

fsolve completed because the vector of function values at the initial
point is near zero as measured by the value of the function tolerance,
and the problem appears regular as measured by the gradient.


ans =

     []
function dUdt = cfunc(t,U_steam)
% Variables
global m_ExtraSteam U0_shell V rho_0 U_steam
Ptube_in = 7;
Ptube_out = 6.5;
T_tube_out = 98;
if t<3.3
m_tube = 707.9;
T_tube_in = 91.7;
else 
m_tube = 1131.2;
T_tube_in = 80.1;
end
% Find the Temp with Internal Energy 
T = fsolve(@TfinderwU,100);
deltaH_tube = XSteam('h_pT', Ptube_out, T_tube_out) - XSteam('h_pT', Ptube_in, T_tube_in);        %  Specific enthalpy difference between the inlet and outlet of the tube side, kJ/kg
if ((m_tube*deltaH_tube)/(XSteam('hV_T', T) - XSteam('hL_T', T))) >= m_ExtraSteam
        m_ExtraSteam_Condensation = m_ExtraSteam;
else
   m_ExtraSteam_Condensation = (m_tube*deltaH_tube)/(XSteam('hV_T', T) - XSteam('hL_T', T)); 
end
m_steam_accumulated = m_ExtraSteam - m_ExtraSteam_Condensation;                 % Steam accumulated in the shell side, kg/sec
steam_mass = rho_0*V + m_steam_accumulated*t;
U_condensate = XSteam('uL_T', T);
CpL = XSteam('CpL_T', T);
Qcv = m_tube*CpL*(T_tube_out-T_tube_in);
% Differential Equation
dUdt = (-Qcv + m_ExtraSteam*U0_shell - m_ExtraSteam_Condensation*U_condensate - U_steam*m_steam_accumulated)/steam_mass;
end
    % Temperature Finder Function with Internal Energy 
function Obj = TfinderwU(T)
global U_steam
uV_T = XSteam('uV_T', T);
Obj = U_steam- uV_T;
end

Talha 2024-8-20

Thanks. Most likely, this happens because of the global variable "U_steam" which is empty at the start.

请先登录，再进行评论。

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

Must return a column vector error when using "ode45"

2 个评论
显示无隐藏无

回答（0 个）

另请参阅

类别

标签

产品

版本

Community Treasure Hunt

Must return a column vector error when using "ode45"

2 个评论 显示 无隐藏 无

回答（0 个）

另请参阅

类别

标签

产品

版本

Community Treasure Hunt

2 个评论
显示无隐藏无