Solving the equations of Simscape's Ejector with vsolve

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

[Yifeng Tang]

Haven't got a chance to review your code yet, but I saw your comment in another MATLAB Answer thread, too. Attached is a Simscape model that may be useful.

With the default ejector parameters, primary at 5MPa, 293.15K, secondary at 0.3MPa, 293.15K, and outlet at 0.5MPa, I'm seeting the p_mix, in Simscape Result Explorer, to be ~0.158MPa.

When I set Pp=1MPa, Ps=0.1MPa, as you suggested, and assume Pout=0.1MPa, I see p_mix = 0.0528MPa.

Simscape Onramp is a good resource for getting started with Simscape. It doesn't discuss ejector, of course.

https://matlabacademy.mathworks.com/details/simscape-onramp/simscape

[Yifeng Tang]

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]

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]

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

[Yifeng Tang]

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]

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]

Dear Yifeng,

Thank you very much for the explanation!