modelica / ModelicaStandardLibrary

Free (standard conforming) library to model mechanical (1D/3D), electrical (analog, digital, machines), magnetic, thermal, fluid, control systems and hierarchical state machines. Also numerical functions and functions for strings, files and streams are included.
https://doc.modelica.org
BSD 3-Clause "New" or "Revised" License
471 stars 168 forks source link

Poor scaling of curve parameter s in Modelica.Fluid.Interfaces.PartialLumpedVolume #2717

Open ErikHenningsson opened 6 years ago

ErikHenningsson commented 6 years ago

Poor scaling of the curve parameter s in Modelica.Fluid.Interfaces.PartialLumpedVolume when portsDataheight[i] = 0 causes problems for mode switch detection. Especially, when switching from the mode // ports[i] is above fluidLevel, preventing outflow_

The curve parameter s is very small when in this mode but rather large when in other modes. This makes detection of zero-crossings difficult.

Problems can for example occur when simulating Modelica.Fluid.Examples.Tanks.ThreeTanks with tank levels starting at 10 m, 0 m, and 0 m, respectively. It may happen that the flow into tank2 starts very/too late.

The equation for s in the above mode is as follows: s[i] = (ports[i].p - vessel_ps_static[i])/Medium.p_default*(portsData_height[i] - fluidLevel); The pressure variables are normalized, but the height difference isn't. I suggest including a normalization of the height difference as well.

beutlich commented 6 years ago

with tank levels starting at 10 m, 0 m, and 0 m, respectively

-10m?

ErikHenningsson commented 6 years ago

with tank levels starting at 10 m, 0 m, and 0 m, respectively

-10m?

No, indeed +10 m. This will cause the other two tanks to fill up. However, the fill up of tank2 sometimes starts later than expected due to difficulties locating the zero-crossing.

beutlich commented 4 years ago

@casella Is this a shortcoming or rather an enhancement?

casella commented 4 years ago

The equation for s in the above mode is as follows: s[i] = (ports[i].p - vessel_ps_static[i])/Medium.p_default*(portsData_height[i] - fluidLevel); The pressure variables are normalized, but the height difference isn't. I suggest including a normalization of the height difference as well.

I think the real problem is slightly different. In one branch of the equation we have

s[i] = fluidLevel - portsData_height[i]

while in the other we have

s[i] = (ports[i].p - vessel_ps_static[i])/Medium.p_default*(portsData_height[i] - fluidLevel)

In the first branch s has the dimension of a level difference in meters. In the second branch it has the dimension of a pressure difference divided by the default medium pressure times a level difference, so it's again a level difference. The scaling is dimensionally correct, the problem is I guess in the scaling factors.

I need to figure out how to improve this

casella commented 4 years ago

@ErikHenningsson, I also see a problem starting at level = 0, since the energy balance equations becomes singular in that case.

casella commented 4 years ago

I performed some more experiments. If you start with level1 = 10, level2 = 1, level3 = 0, you see that the level in tank2 first goes down a bit (because there's more water flowing out to tank3 than in from tank2), until level3 increases enough, without getting to zero. This behaviour persists until level2 = 0.12, approximately. If you go below that, e.g. level2 = 0.11, the level hits zero after some time. At that point, if you look at the mass flow rate, you see a kind of instantaneous latch-up transient: the mass flow rate becomes zero, and stays there either indefintely or for a long time, even though the actual level of the two adjacent tanks (counting for the heigh differences of the two tanks) are much higher.

The only explanation I could come up with is that the branch equation

s[i] = (ports[i].p - vessel_ps_static[i])/Medium.p_default*(portsData_height[i] - fluidLevel)

which holds for negative s[i], tries to accommodate two phenomena that can cause the "check valve behaviour" of not allowing flow reversal: one is having the level below the port height, and the other one is having a negative delta-p across the port; the combination is handled by means of a product.

In the transient discussed above, as I understand, there is a combination of factors causing the two phenomena to take place at the same time. This somehow makes the zero crossing quadratic, instead of linear, and thus causes some kind of numerical lock-up around s[i] = 0, which leads to unphysical behaviour.

Maybe we should have two separate s coordinates to describe those two phenomena separately, and eventually combine them in some clever way which is not a product. I currently have no idea how to do that.

I wish I had more time to play with this very tricky and challenging problem, but unfortunately I have more pressing tasks and deadlines ahead.

I'll leave this ticket open, in case someone wants to pick it up, the issue is far from being resolved.