Simscape Fluids component behavior at low pressure (< 1 Pa)

I am modeling a gas flow system in Simscape Fluids (gasses) with pressures ranging from atmospheric down to vacuum conditions (~1e-5 Pa) at upstream and vacuum (~1e-5 Pa) at the downstream both using reservoirs. Then inbetween is a pipe basically but with a constant volumetric flow using a custom simscape component.

The system works fine above 1 Pa, but when the pressure drops below 1 Pa, my custom component seems to not get the constant volumetric flow rate which I set. While this does work for all the pressure above it (I can't really give the exact code which I know is a bit hard but I have nothing special (even more so) around 1 Pa), it does not have the constant volumetric flow rate over the pressure range going from 1Pa down. For this to work at all I changed the minimum pressure in gas properties to a very low value and have a perfect gas close to air (air has the same behavior the default gas properties variant).

I suspect the issue comes from the way Simscape handles density calculations using a log function. Specifically, I noticed that instead of using the standard equation:

Simscape applies a log-based approach like:

I also found that Simscape uses a protected log function:

function y = logProtected(x, x0)

% Returns log(x) for x >= x0 and a linear extension for x < x0

definitions

y = if ge(x, x0), log(x) else x/x0 - 1 + log(x0) end;

end

end

where x0=1 in components like a pipe (seemingly 1 Pa), which I suspect is causing numerical issues at very low pressures.

Observed Issue:

• When p>=1→ Model behaves as expected.
• When p<1 → Volumetric flow of the component shows still the same constant value which is used to calculate the mass flow. When using a flow sensor with volumetric flow at actual conditions (upstream) shows then decreasing volumetric flow as pressure of the upstream goes down towards the downstream pressure or 0.
• When trying solution or different approach to make the same happen. Steady-state fails to converge, especially in the equation:, seems to produce invalid values at low pressure. One of possible solutions was removing the log entirely and using the original equation which works and gives similar results for p>1 but does not converge below 1. Another try was removing logprotected and using log(x) even when x<1 for all the pressure calculations (and in turn density) but this also fails to converge at pressure below 1 Pa.
• The default pipe still does give results for pressures below 1 Pa, but I can't say a lot about those results as they are supposed to decrease with decreasing pressure. So any problems wouldn't necessarily show up as well as it does for my purpose/use.

Note that when I say 1 Pa, it mostly for convergence is often a bit below it but already gives weird/bad results below 1 and then break somewhere between 0.5 and 1 Pa.

Questions:

1. Why does Simscape use logProtected(x, x0=1) instead of a direct logarithm for density?
2. How can I ensure that my model works for vacuum conditions (<<1 - 1 Pa) without breaking at the log function transition?
3. Are there recommended settings or numerical stabilizations for low-pressure gas flow modeling in Simscape?
4. If these things are not possible is there any workaround for the use case of low pressure/vacuum and custom component (with constant volumetric flow over pressure ranges)?

My system is Windows 11 with Matlab 2024a: 24.1.0.2689473 (R2024a) Update 6 (but I would think that this shouldn't matter as it is all quite new and simscape does not have (many) changes towards newer versions )

It is a rather complex problem it seems and I tried to explain it as short as possible while being complete but let me know if anything is unclear.

Thanks in advance for any insights!

[Torsten]

I think you should address these questions directly to the developers - they are too specific to be answered by forum members.

https://uk.mathworks.com/support/contact_us.html

[Isaac Ito]

(Moved my comments to Answers.)

[Yifeng Tang]

Is it possible for you to reproduce this issue with blocks in the foundation library and the share the model? We'll need to rule out that the issue is due to (1) certain equations, or certain way the equations are written, in the custome component you mentioned, and (2) steady-state solve, which I assume you mean the "start simulation from steady-state" in solver configuration.

Regarding the log function, I'm not 100% sure but my assessment is that it's there to ensure the numerical accuracy when computing terms like needed for volume state calculation, and maybe also to help with the Newton iteration solver within each time step. For example, see line 122-124 of the pipe (G) code. When used as an input to a log function, the value of a variable needs to be non-dimensionalized, and the domain parameters like p_unit and T_unit are used to scale the variable. Such domain parameters can't be changed directly via the Gas Properties (G) in the foundation library, but I believe you can make a copy of it and implement as a custom component and attach to your network. See line 419-431 of the .ssc code of this block. Move some of the ?_unit into a parameter section that is not " (Access = private)". You should be able to change their values then. And you were right, the default scaling factor is 1Pa.

But again, let's rule out other possibilities before touching this more fundamental part, as it'll change the behavior of every component you have.

Answer 1

[Isaac Ito]

Try using Gas Properties (G) block in your model.

Depending on your choice of Gas specification parameter in the block, you can find parameters like Minimum valid temperature and Minimum valid pressure. By default, they are 1 K and 1 Pa, respectively. These parameters set thresholds for the solver to use linear approximation if simulation is configured to continue even outside of the specified temperature/pressure ranges. In your use case, you would have to use smaller values for these parameters.

[Paul Hermens]

I was already using gas properties block like said in my post. I already tried setting that to a low value as said in my post. That does however not give the problem itself. As the simulation still works below 1 Pa, however it does not work as expected as it changes suddenly from >1 Pa going to < 1 Pa. So this is not an answer.