Must return a column vector error when using “ode45”
%% 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
%% 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%% 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
%% 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 %% 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
%% 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 ode45, function, error MATLAB Answers — New Questions
