Understanding and Improving Model Performance at Small Mass Flow Rates in Fluid System Models

This paper provides a detailed analysis of the reasons behind the poor simulation performance observed when mass flow rates become very small, commonly referred to as zero mass flow issues. By using simple example models, we effectively demonstrate the underlying causes of these simulation performance issues. We highlight various contributing factors that play a significant role in ex-acerbating the problem. Furthermore, we propose and examine countermeasures to mitigate these challenges. These countermeasures include modifications to the model itself, utilization of available settings in simulation tools, and adjustments to the solver. By implementing and evaluating these countermeasures, we illustrate their impact on improving simulation performance in scenarios involving low mass flow rates.


Introduction
Modelica is a great option to easily create models to simulate complex fluid systems.The created models can be simulated in one of the available tools.Thanks to advancements in algorithms and computational power, it is now possible to model these intricate systems within a short period of time.However, there is a challenge when it comes to systems that contain branches with no flow during certain periods or when the model is used to simulate rampup or shut-down sequences.In such cases, simulation time can drastically increase, resulting in what is commonly known as zero mass flow problems.From our experience in supporting several customers in creating and simulating models using different libraries, this issue causes large problems.Nevertheless, there is limited literature available on this subject.Dermont et al. (2016) presented an analysis of measures to improve robustness of models used to simulate air condition cycles.The demonstrated the impact of different measures including regularization of the mass flow pressure correlation at low mass flow rates, heat transfer modeling and choice of solver.The study (Li et al. 2020) highlighted the impact of an accurate and fast calculation of the system Jacobian matrix as part of the solution process.Qiao and Laughman (2022) analyzed the impact of different measures to improve perfor-mance of air conditioning models at low mass flow rates.They came up with a new regulation scheme for the pressure drop mass flow correlation around zero mass flows and an analysis of heat transfer handling methods.
In our opinion these article only scratch the surface of the numerical reasons which cause the slow simulations at zero mass flow rates.We strongly believe that gaining a deeper understanding of this phenomenon is crucial for developing effective strategies to improve simulation time in scenarios involving low mass flow rates.Therefore, the primary objective of this paper is to elucidate the underlying causes behind the sluggish performance observed during simulations with low mass flow rates.
To achieve this goal, we begin by examining and analyzing the solution process of a highly simplified fluid dynamic model.This analysis serves as a starting point to unveil the root causes of zero mass flow issues.Furthermore, we expand this model to incorporate additional complexities, thereby showcasing the impact of increased intricacy on simulation performance.By employing all these models, we illustrate the application of various countermeasures aimed at mitigating the challenges posed by zero mass flow problems.

Analysis for simple model
For demonstration of the underlying phenomenon we use a simple model of an isothermal ideal gas in a volume with the variable pressure p connected over a flow resistance to a boundary with fixed pressure When including a quadratic flow model with the parameter c but neglecting the dynamics of the momentum flow as often done for gas flows we can derive the single ordinary differential equation In order simplify the equation we introduced ∆p = p b − p and τ = V R•T •c .Equation ( 4) could be integrated using an DOI 10.3384/ecp204189 Proceedings of the Modelica Conference 2023 October 9-11, 2023, Aachen, Germany implicit Euler scheme in which p t i denotes the solution at the current and p t i−1 the solution at the previous time step.Equation ( 6) is a nonlinear equation which has to be solved to determine the pressure at the current step t i .We have brought the equation to residual form.So we have to find the root of the residual function R(p).For the solution solvers usually apply a Newton method, in which the linearized equation is repeatedly solved Here the prime denotes the derivative and p j t i is the updated solution for pressure at time step t i in iteration j.The iteration is performed starting with an initial guess p 0 t i for the solution.Solving the equation requires solving the linear system R ′ p j−1 t i . For a single equation this can be done without much effort, but for a system of differential equations, the derivative of the residual equation system becomes a matrix (which is closely related to the Jacobian matrix of the system function F).Calculation of the derivative/Jacobian matrix and calculation of a decomposed form for solution comes at high computational effort.Additionally, in many systems the matrix does not change too much while advancing in time.Therefore, solvers usually never update the derivative during the iterative solution process for a time step and use (R ′ (p j t i )) −1 = (R ′ (p 0 t i )) −1 (known as Chord method (Kelley 1995)).Additionally, solvers try to use the same derivative (R ′ ( p)) −1 for solving multiple time steps.With this assumption the following iteration scheme is used with a constant R ′ ( p)) −1 Figure 1 visualizes an exemplary successful iteration for the example problem from equation ( 4).The slope of the linearized approximations (red lines) do not match the actual local slope of the residual function (blue curve).Nevertheless, the iteration schemes converges, but it has lost its quadratic convergence due to the constant derivative.But as shown in Figure 2 the iteration might fail.The iterations schemes drifts off from the actual root of the residual equation and finally circles around the solution.In the following we will analyze this phenomenon in more detail and demonstrate that the iteration with not-up-to-date derivative becomes more difficult to solve when mass flow rate becomes small.If the solution diverges, as shown in Figure 2 or if due to the slower convergence no solution is found within a given number of iterations the step will be rejected.As reaction the solver will request an update of the Jacobian matrix and repeat the solution process.If the solution still fails, it might be necessary (as we will show later on) to reduce the step size.Then a solution with the outlined method can only be achieved for very small time steps and if the Jacobian matrix is updated at every time step.In this setting the time step is not chosen to fulfill the chosen tolerance anymore, but to make the system solvable with the approximated Jacobian matrix from the initial guess.This slows down the solution process and is the cause of slow models performance what is named zero mass flow issues.
Let us take a closer look on the iterative solution process of Equation 4. The Banach fixed-point theorem states that an iteration scheme Φ(p j t i ) will converge if it is contractive.In that case a contraction factor λ < 1 exists and the iterative schemes fulfills For the Chord method we can calculate the contraction factor with By applying the mean value theorem there exists a p ξ ∈ (p j−1 t i , p j t i ) such that For the second step we included Equation 8.As convergence requires λ < 1 we can obtain with Equation 11 This criterion can be used for a more detailed analysis of the cause of zero mass flow issues and how it can be avoided, or its effects can be reduced.
But before applying this criterion to our demonstration problem Equation 1, we want to take a general look on the significance of Equation 12.For the reasons mentioned above a solver will try to reuse the derivative R ′ ( p)) for multiple steps.But the actual derivative in the region of the solution R ′ (p ξ ) might differ from the used derivative R ′ ( p)) by a factor of two as demonstrated in Figure 2. In that case the iterative solution of the nonlinear system will fail.The solver will reject the step and repeat the solution with a reduced step size and an updated value for R ′ ( p)).Unfortunately, the actual criterion for convergence is even stricter than that from Equation 12.If the condition is fulfilled, but with values close to two, the convergence is very slow.The integrator demand convergence within a few iteration, e.g.Hindmarsh et al. (2023), otherwise the step is rejected and the simulation performance decreases.But as the ratio |R ′ ( p)) −1 • R(p j t i )| represent somehow the convergence speed of the method, we can use it for an analysis.Smaller values for |R ′ ( p)) −1 • R(p j t i )| will lead to convergence within fewer steps.As convergence within few steps is demanded, a good initial guess is of great importance.The closer this initial guess is to the actual solution, the more likely is convergence of the method to the actual solution within the demanded solution tolerance.This has three consequences: 1. Methods which have a good approximation or forecast used for the initial guess will not be affected as much by the zero mass flow issue like method with no good forecast.
2. Reducing the time step improves the approximation of the initial guess.Therefore, reducing the step size has two positive aspects: it improves the initial guess and the ratio |R ′ (p 0 t i )) −1 •R(p j t i )| will be closer to one.But obviously this comes as the price of a slower simulation.
3. The stricter the tolerance demanded for the solution, the harder it becomes to reach the tolerance within the desired number of iteration.So decreasing the tolerance might cause more problems at low mass flow rates.
After that general analysis, we will take a closer look on that criterion for our demonstration problem.

Analysis of Convergence for Demonstration Problem
If we apply Equation 12to our simple problem in Equation 1, we get First of all, one can see that as p − → p b results in p − → p ξ the method has only a chance to convergence for ∆t − → 0. Therefore, in almost all libraries the square-root in the mass flow pressure relation is replaced by an approximation which has a finite derivative when crossing zero.For the regulation named regRoot from the Modelica Standard Library (Modelica Standard Library 4.0.0 2023) eq. 13 becomes But still convergence for p − → p b can become problematic for small values of τ and if ∆p small is not chosen large enough.The solver will find a solution within the demanded number of iteration but only for small time steps.Additionally, the simple demonstration model can be used to explain some general phenomena which can be seen in more complex models concerning if zero mass flow issues occur or not.Therefore, the simple model has been implemented as Modelica model and was solved with Dymola 2022x (Dassault Systemes AB 2023) with the regulated form of the square root mentioned above.Two different solvers from Dymola were used: Dassl and Radau.Dassl as a general purpose solver with good performance and Radau, as it was recommend by Dermont et al. (2016), when encountering zero mass flow problems.Figure 3 shows the number of Jacobian matrix evaluation when solving the problem while varying the pressure level of the boundary p b .The pressure was initialized at the same level as the boundary.After one second the pressure in the boundary was increased or decreased by 10 Pa and the system was simulated for further 10 s to settle.One can see that after exceeding a certain threshold in pressure level the number of Jacobian matrix evaluation dramatically increases.This indicates presence of a zero mass flow problem as the integration steps are rejected, and a new Jacobian matrix must be calculated.This pressure level threshold depends on the chosen solver.Radau generally performs better in this case.But in general, it is surprising that the pressure level has an impact on the solution of the problem.The pressure level of p b is only a constant offset for the solution and therefore on first look it should not have an impact.Additionally, if Dassl is chosen as solver, the number of Jacobian matrix evaluations depends on the sign of the pressure step.When Radau is used, the same number of Jacobian matrix evaluation is required independent of the sign of the pressure step.Both phenomena can be explained if one focuses on how the derivative R ′ ( p) is calculated in the solution process.If no special settings are applied the derivative is calculated numerically presumably with a given relative variation ∆p ε = ε • p with a small number ε.Typical methods are Firstly, if the pressure level increases, the absolute variation for calculation of the derivative increases.But the physical meaning of a changed pressure difference has not changed.For a solution around the steady solution p ≈ p b the variation ∆p ε used to approximate the derivative might flip the sign of p − p b .In that case the absolute value of approximation of R ′ ( p) becomes smaller than the absolute value of the actual derivative (see Figure 4).As Equation 13shows this decreases the speed of convergence.Secondly, the problem itself is symmetric to an increasing or decreasing step of the pressure.But using the approximation Equation 15to calculate the derivative introduces an asymmetry.As only the Dassl solver is susceptible to the sign of the variation it suggests itself that Dassl uses a forward differencing scheme while Radau uses a central difference scheme like Equation 16.
We created an external external function which is called in the model to track the value of states during all models calls: by analyzing the state variation one can identify the times at which the Jacobian matrix is updated.When using Radau as solver the perturbation is applied in both directions while for Dassl only a perturbation in one direction is applied -therefore we can verify the assumption.The central scheme is is more accurate and symmetric, and therefore the results to not show a dependence on the sign of the pressure step.
There are two options to improve model performance for zero mass flow rate.If possible one could use an analytic calculation of the Jacobian matrix to get rid of the underestimation of the derivative.Additionally, or solely the problem could be reformulated to use ∆p = p− p b as state.This measure improves the accuracy of the numeric approximation of the Jacobian matrix.Pressure differences are driving the mass flow rates and therefore the changes in pressure.If the pressure becomes large, the perturbation applied ε • ∆p to calculated the derivative might be in the range of the driving pressure differences and the approximate Jacobian matrix is not accurately enough to solve the system with large steps.Choosing the pressure difference as state, leads to a perturbation much smaller used to calculate the derivative and the applied solution method is much more robust.In both cases the model performs much better as one can see in Figure 3.If one of the tweaks is applied, no issues occur even if pressure level is increased.Additionally, the number of evaluations is independent of the step sign.

Non-linearity as root of the problem
Furthermore, we want to use the simple model to emphasize the importance of the non-linearity of the problem.Dermont et al. (2016) gave the time constant of the problem which increases when decreasing the mass flow rate as reason for the poor model performance for small mass flow rates.We want to demonstrate that the time constant itself is not the root of the problem.For this test we run the model with Radau as solver and the pressure level of 1×10 6 Pa with a step of 10 Pa.For the second run the regRoot mass flow pressure relation was replaced with a linear relation which has the same slope for vanishing mass flow rates (Figure 5).Due to the linearity this slope remains the same for all mass flow rates.The slope of the mass flow pressure relation determines the time constant of the problem.The comparison is designed in a way that the time constant of the non-linear problem is at maximum the same as the time constant of the linear problem. Figure 6 shows the solver step when both models are simulated with the same settings.Though the linear problem should be tougher to solver from a simple perspective of the time constant, it is actually the other way round.A linear model cannot face zero mass flow  issues as the Jacobian matrix calculated at any position is valid for the whole domain.As a result the linear problem can be solved with much higher time steps, though it actually has the higher dynamics.Actually, the presented findings could have been derived from Equation 12.But the results are presented as they impressively show which solver steps are possible coming from the pure dynamic of the problem.The non-linearity of the problem which results in a change of the time constant, causes the Chord method to fail and as a result the possible time steps of the integration method cannot be exploited.

Extended model with transported scalar
Obviously, models are normally much more complex.One additional level of complexity are additional equations to be solved e.g.balances for balance, mass, composition etc.The form of these equations has a crucial impact on the resilience of the model against zero mass flow issues.In order to demonstrate and analyze this impact we extend our simple models with an additional equation The solution of equation R(ψ t ) = 0 depends on G(∆p) = H(∆p) • F(∆p) and hence on the flow direction through H(∆p) and on the actual mass flow rate through F(∆p).It is given by Though in our model it would be possible to solve the equation for pressure and the equation for the passive scalar sequentially, in a Modelica tool these equations are solved simultaneously.So Equation 17is added to the model from the previous section and solved again while tracking the number of Jacobian matrix evaluations as it was plotted in Figure 3.This time only Radau and a step of 10 Pa was used.The results are show in Figure 7.For reference purposes the results for only the pressure equations are shown additionally.Adding the equation for the transported scalar causes the zero mass flow issue to occur at lower pressure levels and from that pressure level additional Jacobian matrix evaluations are caused compared to the reference case.When integrating the solution of Equation 17 depends on the actual mass flow rate.Therefore, the iteration can only converge after the Equation 1 has converged reasonably close to the actual solution.If Equation 1 alone converges only just within the desired number of iteration, the system will not converge, causing the zero mass flow issue occurring at lower pressure levels.Figure 8 shows an exemplary course of the residuals of the two equations over the iterations.One can clearly see that the residual of the scalar (Equation 17) even increases before it starts to converge.But convergence rate is slow, while the mass flow rates oscillates around its final value.Figure 7 contains measure which are intended to improve convergence around zero mass flow rate.One idea is to introduce diffusion into Equation 17 by adding a transport term which is independent of the mass flow rate We use a artificial diffusion constant µ.
The second idea is to introduce a (nonphysical) regulation of the convective transport as depicted in Figure 9.The convective transport with is set to be zero if the mass flow is regulated.After the mass flow has left the regulation the convective flow is ramp up to its actual value.
From Figure 7 we can see that the pure diffusion does not improve the system convergence around zero.No improvement was observed for reasonable big values of the artificial diffusion constant µ.The problem is that the solution for ψ still depends on the mass flow rate, though it is shifted.Nevertheless, the scalar converges only after the mass flow has converged.For an improvement the solution should be independent of the mass flow rate for very small mass flow rates.That is the aim of the second measure, with the additional regulation of the convective flow.This method has to be implemented carefully to avoid a violation of conservation.When this measure is applied the system is more robust against zero mass flow issues, but it still performs worse than the model with only mass conservation.

Extended model with multiple volumes
In this step we want to extend the model from section 2 to a model with multiple connected control volumes.

Pressures
For the pressures we obtain where we have extended the notation from Equation 4to the control volume label k = 1...N with ∆p . Following the procedure in Equation 14, we now approximate In a numerical experiment with n = 9 the control volume V [k] was chosen uniformly V [k] = V except for volume five, in which it was chosen to be a tenth of the other values V [5] = 0.1 • V , such that τ[5] = τ/10.The other parameter varied is ∆p small [k].Two strategies are applied: in the first strategy all regularization intervals have the same size and in the second strategy all have the same size except for that of k = 5.The value is chosen to be ∆p small [5] = 10 • ∆p small .The results are given in Figure 10.It is important to mention that the given level of ∆p small corresponds to the value applied for most of the equations.The first important thing which one can see is that increasing the regularization interval will improve the performance until a certain threshold.After that one could say that the zero mass flow problem is not present anymore.But the increase in ∆p small comes at a certain price, the pressure drop at low mass flow rates will be computed incorrectly.If this is an issue, the second strategy might be an interesting option.It is only required to apply a larger linearization for flow models which are connected to pressure state with the large values of τ.So the overall accuracy of the model will be better.

Passive Scalars
Using the same notation as in the previous section we can extend Equation 18for the passive scalar to obtain , where we have used the symmetry properties of the Heaviside step function H and the function F due to Equation 4.

Generalized convergence criterium
In this section we would like to extend the computation of λ in Equation 11 and the resulting convergence criterium in Equation 12to general state space systems with N states For a given time step t i we have and for an implicit integration algorithm it holds that with the non linear residuum equation which can be solved iteratively by a Newton-Raphson scheme Here J ⃗ R denotes the Jacobian of the residuum vector ⃗ R(⃗ x j t i ).In analogy to Equation 8the Chord method attempts to solve Equation 29 with the inverse Jacobian fixed at some state ⃗ x: Now in analogy to Equation 9 this iteration converges if it is contractive, that is for some λ < 1 it holds that In order to further develop this expression, notice that we re-write the difference on the left as the result of an integration along a straight line where J ⃗ Φ ⃗ x(s) denotes the Jacobian of ⃗ Φ at position⃗ x(s).Now clearly the tangent vector d⃗ x(s)/ds is constant along the line and given by d⃗ x(s) Moreover from Equation 30 it follows that with 1 denoting the N-dimensional identity matrix.So we can formally re-write Equation 32as with the matrix A given by With these ingredients we can replace Equation 31 similarly to Equation 11by Here ∥A∥ = ∥A∥ 2 denotes the spectral norm of the matrix A as induced from the L 2 vector norm ∥⃗ x∥.The spectral norm is defined as the square root of the largest eigenvalue of A T A, with A T the conjugate transpose of A. We can further estimate ∥A∥ as follows Although Equation 38holds for any matrix norm, for practical application in the context of Equation 31 one has to use the spectral norm ∥J∥ 2 .

Explicit application
The bound obtained in Equation 38 can be used in order to give conditions for convergence and also recommendations for the regularization parameters of the regRoot and SM functions contained in F(∆p) and G(∆p).

General case
The matrix elements J ⃗ R (⃗ x)[I, J]of the Jacobian J ⃗ R (⃗ x) can be computed from Equation 28as follows: Here δ [I, J] denotes the Kronecker delta.
For the coupled pressure scalar system we have the time derivatives f [I] given from Equation 22and Equation 24. Let Hence the resulting incidence matrix is tri-diagonal for the pressures and tri-diagonal for the scalars with an additional tri-band between pressures and scalars.This property may be used in order to arrive at an estimate to the spectral norms of Equation 38.An explicit computation of the eigenvalues can in principle be avoided by approximating the spectral norm by the general property: where ∥A∥ 1 = max j ∑ N i=1 |a i j | is the maximum absolute column sum norm of A and ∥A∥ ∞ = max i ∑ N j=1 |a i j | is the the maximum absolute row sum norm of A. This also avoids computation of A T A. An explicit symbolic analysis for the pressure system in the fashion of section 7 in J. Brunnemann (2008) will be subject to future work.

Single pressure
For the case of a single pressure p[1] and boundary pressure p b the Jacobian J ⃗ R (⃗ x) has only one single element: For simplicity of notation we have dropped the discretization index on the right hand side.Plugging this into Equation 38 we obtain max In order to fulfill this inequality it must hold that This re-produces the result of Equation 12.

Single pressure and passive scalar
For the case of a single pressure p[1] with boundary pressure p b acoupled to a single passive scalar ψ[1] with boundary ψ b the Jacobian J ⃗ R (⃗ x) is a 2x2 matrix with three non-zero elements: Here we have left out the discretization indices for p, ψ on the left hand side for simplicity of notation.
In the sequel we will suppress the (s)-dependence of (a, c, d) for simplicity of notation.The above setting implies From the symmetry of M T M it follows that the two eigenvalues λ 1 , λ 2 are real.For the 2D case we can explicitly compute them.However one may also apply Gershgorins circle theorem (Gershgorin (1931) This is of particular use for higher dimensional cases of Equation 38.

Solver Modification
The adaptions shown which should lead to an improved simulation of the model where all on the model side.These modification can be applied by the simulation engineer or library developer and they work independently of the used tool.But these modifications cause deviations between model results and the physical correct or the expected results.These deviations might be tolerable in many cases, but in some cases they are not.Therefore, it would be interesting to have a solution process which completely avoids zero mass flow issues or reduces the impact.Simple Mass Flow (section 2) 7.1×10 1 2.0×10 4 7.0×10 3 1.4×10 5 2.1×10 3 3.4×10 4 Mass Flow + Scalar (section 3) 5.3×10 2 6.7×10 3 3.0×10 4 7.7×10 5 1.0×10 4 3.0×10 4 Multi volume (section 4) 2.2×10 1 1.1×10 4 2.7×10 3 1.3×10 5 1.1×10 3 4.6×10 4 BranchingDynamicPipes 7.3×10 1 7.1×10 1 4.3×10 3 4.3×10 3 3.4×10 3 3.4×10 3 PID 1.9×10 1 1.9×10 1 9.5×10 2 9.5×10 2 8.4×10 2 8.4×10 2 BatchPlant_StandardWater 1.7×10 2 1.9×10 2 6.8×10 3 6.7×10 3 5.1×10 3 5.2×10 3 Overshooting during the iteration process is a known problem even for normal Newton methods, if the starting point is far off from the actual solution and/or if the equation to solve is highly nonlinear.In that case damping of the solution can reduce the required number of iteration or even avoid divergence of the method.As shown above overshooting is as well the problem which causes the zero mass flow issue.Therefore, we tried to apply a damping strategy to the solution process.For this test we used OpenModelica and the Cvode solver as for both the full code is available and can easily be modified.The algorithm is taken from Dahmen and Reusken (2008) and modified to be usable for the Chord method.For a general residual equation a Chord step is done to calculate an update the solution for all states The update is not applied directly.Instead the condition is checked until it is fulfilled with the series λ = 1, 0.5, 0.25, ... using C λ = 1 − λ 4 .When the condition is fulfilled the step is applied and the same procedure is repeated for the next step ∆x j+1 t i .The right hand side of Equation 48 is (a fraction) of the norm of the full step calculated from Equation 47.The left hand side is the next full step which is calculated if the current damped step is applied.Therefore, by using this method we damp the calculated step until the norm of the next step is smaller than the current.As the size of the step depends of the residual R(x j t i )), we enforce the residual to reduce.Checking the condition comes with no relevant extra computational effort.The price of a recalculation after a rejected Chord update costs the same as a normal step.So if the rejected steps are included in the total number of Chord updates, one could design a method which comes at almost no extra costs.
The method was applied to the example problems described above in configuration in which the zero mass flow issue occurs.The performance of the damped strategy is compared to the normal algorithm.The results are summarized in Table 1.Damping the steps causes a significant reduction in Jacobian matrix evaluations and functions calls.The number of integrator time steps is reduced as well, though in case of the problem from section 3 the reduction of integrator steps is not a high as in the other.But function calls and the decomposition of the Jacobian matrix are the main expenses during integration.Therefore, using the damped solving approach significantly reduces the solution time in all the cases.Beside the some models from the Modelica Standard Library were simulated with both methods, though zero mass flow rate was not a particular issues for these models.The performance indicator numbers are very similar.The slight difference can result from differences in the solution and the fact that with the used settings damping steps were not counted in the total step count for an iteration.Therefore, the number of function with damping exceed the number of calls without damping though less steps were taken and less Jacobian updates were required.These minor deviation could be tolerated.All in all the method looks promising to reduce the impact of zero mass flow issues.

Summary and Outlook
In this paper we demonstrated the reason for the slow integration process of fluid system close to vanishing mass flow rates: Due to the high non-linearity of the problem the solution of the non-linear system for calculating an integration step fails, as the simplified chord method is used.The system can only be solved for steps smaller than the steps required from the dynamics of the system itself.And additionally the Jacobian has to be updated almost every step.Therefore, more step which are itself even more computational expensive are required.The problem can be aggravated due to inaccuracies in the calculation of the Jacobian matrix depending on the choice of the states.Though we focused on models of fluid systems, this phenomenon is not restricted to this kind of problem.It occurs for any model with high non-linearity in state where a system should come to rest.Some approaches which help to avoid the issues were given: If possible the non linearity of the model should be reduced e.g. by linearization close to the zero mass flow rate.The linearization interval can be adapted to the local model behavior.Additionally, it is helpful to make other variables independent of the solution of the mass flow rate if the mass flow rate becomes small.Furthermore, it is advisable to improve the accuracy of the Jacobian matrix calculation e.g. by modified state choice or by a analytic calculation of the Jacobian matrix.With the knowledge of the exact cause of the problem it might be possible to develop additional improvements or enhance the existing once.For example it might be possible to develop an automatic calculation of the linearization interval for the relation of mass flow rate and pressure drop.
As the zero mass flow issue comes from the non-linear solution process, a modification of that process was suggested as well.By including a damping procedure in the solution process, the solution process for small mass flow rates can be improved.Unfortunately, this method only improves solution process if the iterative solution "overshoots" the actual solution.Non-linearity might also lead to a slow iteration progress without overshooting: for example when calculating the root of ṁ2 = 0.If a dynamic momentum balance is used, the problem to solve is in that nature and a damping of the steps is not helpful anymore.This problem should be analyzed as well as it is done here for models without mass flow rate state.It might be possible to improve the solution process as well, e.g. by including factors or solution processes which do an approximated update of the (inverted) Jacobian matrix during iteration (e.g.Broyden method).The usage of these methods might be useful in the presented cases as well.

Figure 1 .
Figure 1.Example of a converging iterative solution for the example problem.

Figure 3 .Figure 4 .
Figure 3. Number of Jacobian matrix evaluations for the example problem for different settings.

Figure 5 .
Figure 5. Linearized and actual non-linear mass flow pressure relation used in the example model to demonstrate the impact of the non-linearity.

Figure 6 .
Figure 6.Solver step size for example model with linear and non-linear mass flow pressure relation as shown in Figure 5.

)
Now we introduce the abbreviation ∆ψ = ψ b − ψ and observe that sign(m flow ) = sign(∆p) in order to use the Heaviside step function H(∆p).Moreover from Equa-tion 4 and Equation 2 we can replace m flow = c • F(∆p) and m = V R•T • p. introduced G(∆p) = H(∆p) • F(∆p)) and again use τ = V R•T •c .The residual equation becomes the form

Figure 7 .Figure 8 .Figure 9 .
Figure 7. Number of Jacobian matrix evaluations for example model with (w/) and without (w/o) transported scalar and different settings.

Figure 10 .
Figure 10.Number of Jacobian matrix evaluations for a multivolume model with and without adjusted linearization interval.

Table 1 .
Performance key figures for the solution process for the presented models with and without damping strategy.