Solving the equations of Simscape's Ejector with vsolve

Question

Paris 2024-6-13

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/2128131-solving-the-equations-of-simscape-s-ejector-with-vsolve

评论： Paris 2024-6-17

Dear Everyone,

I have written a Matlab script based on the equations of the Ejector in Simscape (https://www.mathworks.com/help/hydro/ref/ejectorg.html). I am not good at using Simscape and I am having difficulties running the Ejector in Simscape.

My script calculates Pm (pressure after mixing ) using fsolve. I have modified the equations a bit as I think γ should be different for the two streams, therefore I have kp (γ of the primary stream) and ks (γ of the secondary stream). I also included an iteration to calculate the stagnation temperatures of the two streams since the documentation of the Ejector states that the block assumes that Tp and Ts (temperatures of primary and secondary streams) are equal to their primary flow stagnation temperature. These temperatures are needed for the calculation of Tm (temperature of mixing). The stagnation temperature depends on Mm (Mach number of mixing) which depends on Tm which in turn depends on the stagnation temperatures.

I am considering a case where the primary stream is CO2 at 1MPa and 20°C and the secondary stream is air at 0.1 MPa and 20°C, and I calculate Pm = 0.0012 MPa and Mm = 0.033. I think these results are wrong, could someone please confirm this? It seems that my solver calculates the same Pm regardless of Pp (primary pressure) and a Pm always close to Ps (secondary pressure).

I am attaching the m files needed to run the script below. "ejector_solver.m". "Ejector_equations_star" uses the same equations as in "Ejector_equations" but replaces Pm with Pm_star to calculate the critical diffuser outlet pressure.

Also, does anyone know how the cross-sectional areas of the three ports are used? They need to be specified for the Ejector in Simscape but I am not sure where they appear in the equations.

I am grateful for any help you can provide.

Best regards,

clearvars; close all; clc;
% Primary stream
Primary_fluid = 'CO2';
Pp = 1; % [MPa] Pressure of primary stream
Tp = 20 + 273.15; % [K] Temperature of primary stream
kp = 1.3589; % Primary stream's ratio of heat capacities 
rhop = 19.098; % [kg/m^3] Density of primary stream
% Secondary stream
Secondary_fluid = 'Air.ppf';
Ps = 0.1; % [MPa] Pressure of secondary stream
Ts = 20 + 273.15; % [K] Temperature of secondary stream
ks = 1.402; % Secondary stream's ratio of heat capacities 
rhos = 1.189; % [kg/m^3] Density of secondary stream
% Gas constant
R = 0.2871; % [kJ/(kg K)] 
% Parameters of ejector (same as in Simscape)
Spt = 1e-4; % [m^2] Throat area
Spn_to_Spt = 3; % Area ratio of nozzle exit to throat
Sm_to_Spt = 8; % Area ratio of mixing chamber exit to throat
eta_p = 0.95; % Efficiency of primary flow through nozzle
eta_s = 0.85; % Efficiency of secondary suction flow
eta_e = 0.88; % Efficiency for primary flow expansion
eta_m = 0.84; % Efficiency for mixing
% Initial guess for Pm
Pm_initial_guess = 0.001; % [MPa]
% Use fsolve to find Pm
options = optimoptions('fsolve', 'Display', 'iter');
Pm = fsolve(@(Pm) ejector_solver(Pm, Pp, Tp, kp, rhop, Ps, Ts, ks, rhos, R, Spt, Spn_to_Spt, Sm_to_Spt, eta_p, eta_s, eta_e, eta_m), Pm_initial_guess, options);
                                             Norm of      First-order   Trust-region
 Iteration  Func-count     ||f(x)||^2           step       optimality         radius
     0          2         2.30732e-09                        1.38e-05              1
     1          4         7.88222e-11    0.000167671         1.42e-06              1
     2          6          4.0154e-12    5.53742e-05         1.68e-07              1

Equation solved.

fsolve completed because the vector of function values is near zero
as measured by the value of the function tolerance, and
the problem appears regular as measured by the gradient.
disp(['Pm is: ', num2str(Pm), ' MPa']);
Pm is: 0.001223 MPa
% Calculate equations with Pm
Ejector_equations
% Critical mode operation
Pm_star = min(Pp,Ps) * (2 / (km+1))^(km / (km-1)); % Threshold that causes choking in both primary and secondary flow
Ejector_equations_star % Same equations to calculate the critical diffuser outlet pressure by replacing Pm with Pm_star
% Evaluate if critical mode operation
if Pd < Pd_star
    Pm = Pm_star;
else
    Ejector_equations % Repeat equations to save the correct values
end
disp(['Mm is: ', num2str(Mm), '']);
Mm is: 0.032822

3 个评论
显示 1更早的评论隐藏 1更早的评论

Yifeng Tang 2024-6-14

For the "cross-sectional areas of the three ports", are you talking about these ones in the red box? Not the throat area on the first row, right?

Please kindly clarify and I can provide an answer in the Answer section.

Paris 2024-6-14

Dear Yifeng,

Thank you very much for your reply and for the Simulink file. My code is wrong as my mass flows are thousands of times lower than those in the Simulink file. Also, I did not know that the outlet conditions of port B need to be fixed.

The cross-sectional areas are those in the red box, it seems to be their value does not appear in the equations written in https://www.mathworks.com/help/hydro/ref/ejectorg.html. If you could please tell me how they are used then I think I will be able to correct my solver.

Best regards,

Paris

请先登录，再进行评论。

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

Answer 1

Paris 2024-6-16

0
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/2128131-solving-the-equations-of-simscape-s-ejector-with-vsolve#answer_1472496

Dear @Yifeng Tang,

Thank you very much for your comment and for the Simulink file. I did not know that the outlet conditions of port B needed to be fixed.

The cross-sectional areas are those that you specified in the red box, it seems to me that their value do not appear in the equations written in https://www.mathworks.com/help/hydro/ref/ejectorg.html. I have varied these values from 0.1 to 1 m^2 in the SImulink file that you created, but I do not see any significant variations in the mass flows.

Also, I noticed that if I open the flow rate sensors under the model tree structure of the Simscape Results Explorer, the "M" mass flow values are 1000 times lower than those displayed by the Scope, is this related to the cross-sectional areas?

Best regards,

Paris

3 个评论
显示 1更早的评论隐藏 1更早的评论

Yifeng Tang 2024-6-17

The port area parameter appears in MANY components of fluids-related domains. They are mainly responsible for how the components interact with neighboring components.

(1) in the energy flow through the ports (one of the two domain through variables), there is usually a v^2/2 term for the kinetic energy. The velocity v is computed based on mass flow (the other domain through variable), density (function of pressure and temperature), and port area. For example, if your model contains a huge jump in the port area between two components, the local temperature may also change suddenly, as part of the internal energy has to convert to kinetic energy to conserve the total energy flow.

(2) Smoothing in the numerical scheme for fluids domains. Simscape uses so-called "upwind scheme", meaning the values of properties like temperature, concentration at ports are taken from the flow from the upstream of the port, i.e. the flow that goes "into" the port. When the flow is nearly zero, Simscape uses a diffusion-like term to smooth out any jump in properties at the port between two components, and the influence from the neighboring components has to do with the port area. The larger port are, more influence the properties inside the block has on the neighboring ports. You can see more details on this documentation page:

https://www.mathworks.com/help/hydro/ug/simscape-fluids-numerical-scheme.html

and you can see the Simscape language implementation in the "port_convection" code of each fluids domains (TL, MA, 2P, G).

Hope this helps.

Yifeng Tang 2024-6-17

Regarding your observation on the value of the results from Results Explorer vs. Scope, I suggest that you look in these places:

(1) PS-Simulink Converter: did you set a unit? Default unit of "1" will be interpreted as SI unit.

(2) in Simscape Results Explorer, the y-axis also contains a unit. You can click on the drop down arrow and change that.

Make sure the two units are consistent. I suspect you have g/s somewhere but kg/s at another place.

Paris 2024-6-17

Dear Yifeng,

Thank you very much for the explanation!

请先登录，再进行评论。

Solving the equations of Simscape's Ejector with vsolve

3 个评论
显示 1更早的评论隐藏 1更早的评论

回答（1 个）

3 个评论
显示 1更早的评论隐藏 1更早的评论

另请参阅

类别

标签

产品

版本

Community Treasure Hunt

Solving the equations of Simscape's Ejector with vsolve

3 个评论 显示 1更早的评论隐藏 1更早的评论

回答（1 个）

3 个评论 显示 1更早的评论隐藏 1更早的评论

另请参阅

类别

标签

产品

版本

Community Treasure Hunt

3 个评论
显示 1更早的评论隐藏 1更早的评论

3 个评论
显示 1更早的评论隐藏 1更早的评论