Stability of dynamic fluid transport simulations

Pipeline transport is an efficient method for transporting fluids in energy supply and other technical applications. While natural gas is the classical example, the transport of hydrogen is becoming more and more important; both are transmitted under high pressure in a gaseous state. Also relevant is the transport of carbon dioxide, captured in the places of formation, transferred under high pressure in a liquid or supercritical state and pumped into underground reservoirs for storage. The transport of other fluids is also required in technical applications. Meanwhile, the transport equations for different fluids are essentially the same, and the simulation can be performed using the same methods. In this paper, the effect of control elements such as compressors, regulators and flaptraps on the stability of fluid transport simulations is studied. It is shown that modeling of these elements can lead to instabilities, both in stationary and dynamic simulations. Special regularization methods were developed to overcome these problems. Their functionality also for dynamic simulations is demonstrated for a number of numerical experiments.


Introduction
For stability of numerical simulation of the networks used for the transport of fluids, not the pipes play an important role.Purely pipe networks are known to be unconditionally stable, while the other elements, used for control of pressures and flows, can bring instability to the system.Such elements used to raise pressure (compressors), to drop pressure (regulators) or to define the direction of the flow (flaptraps).A study of the influence of such elements to stability of dynamical simulation of fluid transport networks is the main purpose of this work.
In general terms, dynamic simulation of network problems is reduced to the numerical solution of a system of differential equations.If the time derivatives are set to zero, the system becomes algebraic and describes a stationary solution: dynamic: F (x) = dx/dt, stationary: F (x) = 0. ( Although the stationary system can be solved directly, e.g., using the stabilized Newton method, solving it by integrating until the stationarity regime may represent a more efficient alternative. The stationary system corresponds to the singular points of the vector field dx/dt, and only points of the attracting type represent stable solutions.Points of repulsive or saddle type must be excluded, which is not easy to do in a stationary solution, while a dynamic solution will automatically select stable points.Dynamical systems can have a more complex type of asymptotic behavior, have limit cycles, strange attractors (random behavior), runaway solutions, which also cannot be found in the framework of stationary simulation.In addition, sometimes stationary problems do not have a single solution, a prototype example is an indeterminate pressure in a pipe between two closed valves.Dynamic problems define such solutions trivially, by the history of integration.In particular, the pressure in isolated sections of the network can simply remain equal to the starting value.
A small clarification, the systems we are considering can be non-autonomous, their l.h.s. may explicitly depend on time.This type of system includes time-varying problem parameters such as boundary pressures and flows.In addition, the systems do not belong to ordinary differential equations (ODE), but to the differential-algebraic ones (DAE).Some variables and equations may not have time derivatives.Systems are assumed to be discretized in space.Initially, they belong to 2D partial differential equations (PDE), contain derivatives with respect to time and spatial coordinate along the pipe.Such equations are constrained by boundary conditions corresponding to the network topology.After discretization in time, a system of algebraic equations is obtained: continuous: F (x, t, dx/dt) = 0, discretized: F (x, t, (x − x prev )/dt) = 0. (2 This system should be solved at each time step dt, with respect to the variables x, given the value of these variables in the previous step x prev .Here we have chosen a fully implicit scheme for time discretization.It gives the most stable integration, although it introduces numerical viscosity effects into the solution.These effects require small integration step if the system dynamics is to be studied in detail.At the same time, integration only for the purpose of determining the stationary solution can be carried out in large steps, and here the effects of numerical viscosity only help the convergence of the solution.More advanced integration schemes may be used if needed, such as Runge-Kutta or Rosenbrock-Wanner.In this paper, we draw attention to the fact that a discretized system can degenerate at some step, which leads to a divergence of the solution or a hang of the integrator.Such problems arise in realistic fluid transport simulation problems, and they are directly related to the simulation of the control elements used there.Physical modeling and basic algorithms for solving equations are implemented in our system MYNTS (Multi-phYsics NeTwork Simulator, [1]).The system provides an open and userfriendly modeling in the form of a list of variables and equations for each node and edge of the network graph.A linguistic approach is used, in which the description of the network and the mathematical description of the system of equations are considered as two domain specific languages.A fast translation algorithm maps one representation into another, using the userspecified modeling as a dictionary.The resulting system of equations is then given to a Newtonian solver with a globally convergent Armijo line search stabilizer.
In our earlier works, modeling relevant to fluid transport problems was described: for natural gas [2][3][4], for hydrogen [5], for carbon dioxide [6].A common problem for stationary setting is the degeneration of systems due to marginally stable modeling of control elements, which was investigated in [7,8] using principal component analysis (PCA).In this paper, we will show that some of these problems can be solved in a dynamic approach, although they are still a limiting factor for the simulation speed and require careful adjustment of the integration parameters.
The physical modeling of fluid transport is based on the Darcy-Weisbach friction formula in pipes, in specific implementation with Nikuradse or more advanced Hofer approximation [9,10].The equation of state (EOS) relates the fluid density to pressure, temperature and fluid composition.In the simplest case, analytical formulas are used: for natural gas Papay [11], for hydrogen and carbon dioxide Peng-Robinson and Soave-Redlich-Kwong [12,13].More accurate calculations use complex EOS based on a large number of empirical factors.These EOS include AGA8-DC92 and GERG2008, which are currently ISO standards [14,15].For problems with phase transitions, which include the transport of carbon dioxide, the homogeneous equilibrium model (HEM, [16]) can be used as a first approximation.In complex problems such as pipe depressurization, more sophisticated modeling can be used [17].
Globally convergent methods for solving stationary problems are based on the stabilized Newton method, for piecewise linear problems [18][19][20], for general nonlinear problems [21,22].There are also well-tested standard integration methods [23] for solving dynamic problems.
The study of the stability of gas transport simulations was carried out in the works [24][25][26][27][28][29][30].In particular, in [24] the stability of solutions was considered in the mathematical sense, irrespective of the discretizations performed.In [25], the stability of two numerical solution methods was considered: the method of characteristics and the discrete-space continuous-time approach.[26] considered the modeling of the simplest gas transport networks, with a practical indication of numerically unstable situations.In [27], in application to high pressure steam networks, an iterative fully implicit finite-difference method was developed and a review of other methods with stability considerations was carried out.[28-30] consider a wide class of finite difference schemes for solving gas transport problems, with a deep analysis of their stability.As applications, in the cited works, a single pipe segment is most often considered, less often -networks containing only pipes, and even less often -networks that include control elements in simplified modeling.In our work, we will consider a more realistic modeling of control elements, show that certain combinations of them can destabilize the simulation, and construct regularization methods to overcome these problems.
In Section 2, the modeling of dynamical fluid transport is described in details.In Section 3, numerical experiments are presented and discussed, with a main focus on stability of the integration procedure.Section 4 summarizes the results.

Dynamical equations of fluid transport
Fluid transport networks usually contain a large number of pipes and a relatively small number of other elements.Pipes have stable modeling.Stationary problems for networks consisting only of pipes usually converge in several iterations.The exception is cycles of short thick pipe segments, to stabilize the solver they should be detected and contracted.The main influence on stability is related with the control elements, the modeling of which is discussed in detail below.
Variables: per node -pressure P , mass density ρ, compressibility z; per edge -flow in different normalizations: mass flow Q m in kg/s, molar flow Q ν in mol/s, power flow Q p , normal volume flow Q N .The power flow is defined only for combustible fluids, measured in (mega)Watts, with the reference to a given calorific value of the fluid.The normal volume flow is measured in normal cubic meters per second (or thousands of such units per hour), with the reference to the fluid volume at normal conditions P norm = 1 atm= 1.01325 bar, T norm = 273.15K.In this paper we will use abbreviations where ρ norm is fluid density at normal conditions.
Equations: per node -gas law P = ρRT z/µ, with universal gas constant R, absolute temperature T , molar mass µ; EOS usually defined as a function z = z(T, P ).In this paper we use ideal gas z = 1 and analytical Papay formula [11] relevant for natural gas.For stationary simulations we have also implemented more complex EOS: AGA8-DC92 [14] and GERG2008 [15], they are not considered in this work.The specific form of EOS does not affect the stability of the solver as long as the dependence of ρ(P ) at a fixed temperature is monotonically increasing [2].
Fluid transport through pipes is described by differential equations [5]: where t is time, x is coordinate along the pipe, v is fluid velocity, D is pipe diameter, g is gravitational acceleration, h is height.The first term on r.h.s.describes frictional force contribution, with λ(k/D, Re) defined by Nikuradse [9] or Hofer [10] formula, where Re = 4|m|/(πµ visc D) is Reynolds number, µ visc is the dynamic viscosity, k is the pipe roughness.Mass flow is related with velocity by a formula m = ρvS, where S = πD 2 /4 is the pipe crosssection area.
The physical meaning of these equations: the first equation describes the law of conservation of mass, the second -the law of conservation of momentum.In the general case, there is a third relation that expresses the law of conservation of energy, in this paper we replace it with the isothermal condition T = T s , where T s is the temperature of soil or other environment.This mode is realized at a high heat transfer coefficient between the pipe and the environment.
Discretization: in this paper, we will use the following physically-based collocation scheme.Per node, dynamical Kirchhoff equation is formulated: where V n is the volume assigned to the node, I ne is the incidence matrix for the network graph.
The physical meaning of this equation: the change of the mass in the node is related with mass flows coming into the node (taken with the + sign) and going out of the node (taken with the − sign).In addition, in supply/consumer nodes of Qset type in r.h.s. an external mass flow is also added, with the corresponding sign.In the entry node of the Pset-type, the Kirchhoff equation is replaced by the condition P = P set .Volume V n is defined as a half of the volume of elements adjacent to the given node.For non-pipe elements, a given positive constant volume 2V 0 is assigned, providing the volume V 0 per every adjacent node.Such an algorithm serves as a regularization of the equations.Further, the parameter V 0 will play an important role in the analysis of the stability of the solver.If necessary, pipes are subdivided into sufficiently small segments, in which the discretization used reproduces continuous equations (3)(4).
Further, resistors are defined as short pipes; valves have the equation P in − P out = 0 in open state, m = 0 in closed state, dynamic valve can be defined by the equation where α ∈ [0, 1].Shortcuts are equivalent to open valves.Static valves and shortcuts can be removed/contracted from the graph in a pre-processing stage.
Control elements: the behavior of these elements is described by diagrams in Fig. 1, while the analytical description has a form [2]: max(min(P in − P L , defining polyhedral surfaces of Fig. 1 in the piecewise linear, so-called min-max representation [20]. For flaptraps, the equation defines the well-known linear complementarity problem (LCP), for other control elements, the representation can be constructed as a sequence of LCPs.
Singularities of modeling: the described control element modeling has a problem with regularity [2].The problem is not even that the equations are non-smooth.It is known that Newton's method converges for piecewise linear systems, moreover, in a finite number of steps [18][19][20].The problem is that for the stationary system to be nondegenerate, the equations of the elements of the form f (P in , P out , Q) = 0 must satisfy the generalized resistivity conditions: For control elements, this condition is met marginally, which geometrically corresponds to the normals to the surfaces Fig. 1, directed strictly along the coordinate axes.Under certain conditions on the network topology, such a configuration can lead to degeneracy of the Jacobi matrix of the stationary system.This corresponds to the appearance of zero eigenvalues in the Jacobi matrix, which can be proved, for example, using PCA [7].
Physically, such a property can be demonstrated on the P/Q-undefined configurations shown in Fig. 1.For two sequential regulators located on the QH-face of the control equation, the condition Q = Q H is imposed twice.This leads, firstly, to the degeneracy of the system, and secondly, to the fact that the variable P in the node between the regulators does not enter into any equation.Thus, the intermediate pressure can be arbitrary, as long as the solution is on the QH-face.However, the solution algorithm based on Newtonian iterations will diverge for a degenerate Jacobi matrix, and no solution will be found.For two parallel regulators located on the PH-face, the condition P = P H is imposed twice.Here the system also degenerates, and the balance of flows δQ between the regulators is uncertain.The work [7] also contains examples where similar conflicts are not limited to the nearest neighborhood of control elements, but spread over wide areas of the network.
We also note that the described modeling of compressors belongs to the simplest, so-called 'free' model.The work [4] describes a more complex, 'advanced' simulation that takes into account the individually calibrated characteristics of the compressors and their drives.However, such modeling only includes additional curvilinear faces in compressor diagrams.The normals to these faces satisfy the resistivity conditions strictly, so the modeling of the advanced faces themselves is stable.However, in the full diagram there are pieces of the free model with a marginal position of the normals.The final answer or intermediate iterations may belong to these faces, thus the problem of degeneracy is not excluded by the advanced compressor modeling.
For stationary problems, in [2] we proposed the following regularization algorithm.In l.h.s. of the control equation, a resistive term + s (P in − P out − R s Q) is added, with a small positive factor s and a fixed positive R s (one can choose, for example, R s = 1 bar/(10 3 m 3 /h)).Now the generalized resistivity conditions are satisfied, the Jacobi matrix becomes nondegenerate.PCA detects positive eigenvalues proportional to s .The system has a unique solution, however, adding a regularizing term leads to small deviations from the target faces.In particular, the conditions Q = Q H and P = P H get small deviations.Thus, a tradeoff arises between the stability of the solution (large s ) and the accuracy of the answer (small s ).
Dynamic regularization: there is another possibility for dynamic problems.The undefined P -value in the intermediate node is determined by the integration history.In particular, in completely isolated sections of the network, it can simply remain equal to the starting value.However, for the applicability of this consideration, it is required that the node has a non-zero volume, proportional to V 0 .Then, according to (5), the density remains constant, and due to EOS, so does the pressure.Similarly, an indeterminate balance between parallel elements can be fixed by the kinetic term − d ∂m/∂t introduced in l.h.s. of the control equation.Here the positive quantities V 0 and d are new regularizing parameters.In contrast to the stationary regularization, the influence of the dynamic regularizing terms tends to zero when the stationary point is reached.Thus, the target values in the control equations in the new approach are asymptotically satisfied.
In more detail, dynamic regularization requires special parameter tuning for integration stability.Indeed, after discretizing − d (m − m prev )/dt, the kinetic term gives a contribution of − d m/dt to the equation comparable to the − s R s Q of the stationary regularization.Equating them, we obtain the relation d = s dt(R s /ρ norm ).Choosing here the value s = 10 −2 , which in practice leads to a stable solution of the stationary problem, we obtain the value d , which gives a similar stability of dynamic integration.A similar argument for the V 0 term leads to the formula V 0 = s dt(ρ norm /R s )c 2 s , where c s = (∂P/∂ρ) T is an isothermal speed of sound.Note that the values of V 0 and d in these formulas are proportional to the step dt, which is usually chosen large for solving stationary problems by the integration method.This gives rise to upper bounds from the rate of filling the network.Indeed, if the volume of all V 0 -elements is significantly less than the total volume of pipes in the network, then the network filling time from the initial pressure (for example, P norm ) to the final one (typically 50-100 bar) is determined by the volume of the pipes.In this case, it makes sense to increase V 0 and the dt-step proportional to it in order to speed up the integration process from the starting point to stationarity.If the volume of V 0 -elements begins to exceed the volume of pipes, then the filling time increases in proportion to V 0 , as does the dt-step.In this case, the number of steps stops changing and a further increase of V 0 does not make sense.Similarly, when the kinetic term proportional to d in the control elements starts to dominate in the equations, then the control elements become a bottleneck for filling the network, and further increase of d should be stopped.
There is another theoretical issue related to the choice of dynamic regularizing parameters.Since the regularizing terms are switched off at the stationary point, the position of the stationary point does not depend on V 0 and d .However, it is not known whether the stationary point remains attractive for any choice of these parameters.Further research is required to answer this question.
Infeasible problems are a different kind of irregularity.Simulations are usually carried out in order to test the feasibility of certain scenarios, for example, in planning problems.This implies the existence of impossible problem statements.Such examples are easy to construct for both stationary and dynamic problems.Indeed, let there be a low-capacity compressor at the beginning of a long pipe, and a high-flow consumer at the end.Obviously, the parameters of this stationary problem can be chosen so that it does not have a solution.Another example is a fixed flow consumer drawing gas from a pipe that is closed at the other end.After a finite time, all the gas will be drawn, and the solution of the dynamic problem cannot be continued beyond this point.
For stationary problems, in [2] we have constructed a special unfolding procedure that extends the equations beyond the physical domain, which for natural gas transport usually has the form P ∈ [1, 150] bar, |m| ∈ [0, 1000] kg/s.The extensions are chosen in such a way that the Jacobi matrix of the stationary system is globally nondegenerate.In this case, under certain conditions on the behavior of the equations at infinity [3], the solution of the system exists and is unique.We emphasize that the equations are continued also into a completely non-physical region of negative pressures.The need for such a continuation consists in the fact that intermediate iterations can lead to this region, and even there the system must be non-degenerate, for the convergence of the iterative procedure.In the case that the found unique solution of the system is located outside the physical domain, then the problem statement itself is infeasible.
For dynamic problems, we perform a similar unfolding to prevent the integration procedure from stopping due to the loss of regularity of the equations.If the solution goes beyond the physical domain, it means that certain dynamic scenarios cannot be implemented.However, if the integration is carried out only to find a stationary solution, then the exit of the solution outside the physical region at the intermediate stage can be ignored.It is only important the final stationary solution is located in the physical domain.At the same time, our numerical experiments reveal one more possibility, the solution of a dynamic problem may not have a stationary limit at all, going into a runaway mode.This mode is also an indication of the infeasibility of the problem.
Other types of instabilities.In the work [26], instabilities are noted in the presence of both large and small flows in different parts of the network.Our analysis [8] confirms this conclusion, however, we found such an effect only in the presence of cycles with small flows, while Qsets of different scales do not present problems.The reason is the appearance of circulating flows as independent degrees of freedom and the quadratic behavior of the friction term (∼ m|m|).In the Jacobi matrix, such degrees of freedom enter with a small coefficient (∼ |m|), leading to degeneracy of the matrix.The cycles of thick and short pipes present a similar problem, due to their small hydraulic resistance.Possible solutions to these problems are to introduce a laminar term (∼ m) into the friction force and/or to contract the problematic cycles.In dynamic simulation, the available dynamic terms (∼ dm/dt) in pipes provide a natural regularization.However, when implicit integration with large steps is applied to find a stationary solution, the dynamic terms are effectively switched off, and the degeneracy of Jacobi matrix can reappear.In that case one can (i) integrate with smaller steps; (ii) increase the magnitude of dynamic terms proportionally to the time step (similar to d -factor); (iii) use stationary regularizers.Here one method or a combination of those can be selected.

Numerical experiments
To test dynamic integration algorithms, we use a variety of network scenarios described in Table 1.The simplest of them, dyn1-3, also shown in Fig. 2, are used to test basic integration characteristics such as time resolution, response to abrupt changes in equations, etc.A more complex network N1 is shown in Fig. 3.There is also a set of 85 highly complex realistic networks provided by our industry partner for benchmarking.On these networks, we test our dynamic regularization methods.The results of numerical experiments are shown in Fig. 4.
Scenario dyn1 demonstrates the dependence of the solution on the integration step and the choice of the starting point.In the dyn1a scenario, integration is carried out with a small step, while the fine structure of the solutions, oscillations, is visible.Such setup can be used for a detailed study of the dynamics.In the dyn1b scenario, the integration step is large, the oscillations are smoothed out by the integrator.These settings are more convenient for finding a stationary solution by the integration method.In the dyn1c scenario, the same step is chosen as in dyn1a, but the initial ρ-value is deliberately chosen to not satisfy the equations in the initial state.This scenario illustrates the general problem, if not all DAEs are satisfied in the initial state, then in the first step there is a jump in the solution, forcing the system to fulfill all the equations.This jump is transmitted into kinetic terms, resulting in a high amplitude shockwave in the network.Although the effects of this shockwave eventually fade out, the echo of the initial shock continues to roam the network for a long time, making the approximation to a stationary solution inefficient.This property makes it necessary to choose a starting point that exactly satisfies all equations.It distinguishes the approach from an arbitrary choice of a starting point for a direct solution of a stationary system by Newton's method.We solve the problem of choosing a starting point in the following way.The solution starts from m = 0, P = P 0 = Const throughout the network.It can be seen from the diagrams in Fig. 1 that the regulator equation is satisfied for an arbitrary choice of P 0 , and the compressor equation -for P 0 ≤ P L .In practice, if there are compressors in the system, one can choose P 0 = P norm , since P L never drops below P norm .In this case, all P set = P norm , Q set = 0.The values ρ, z are calculated from the given EOS.This initialization depends on the composition of the fluid, and in this work we have fixed it to pure methane.For stationary fluid transport simulation, we have already implemented variable composition, mixing equations, non-isothermal thermodynamics and complex EOS.We are working on an appropriate upgrade for dynamic simulation.
As a result of the initialization, the starting point satisfies all equations with zeroing of time derivatives, being the exact solution of an auxiliary trivial stationary problem.Next (phase I, network filling) we start to continuously increase Psets towards the required values.Also, all terms in the equation that may conflict with the choice of the starting point (for example, gravity pressure drop ρg∆h) should be initially omitted, then gradually turned on.Then (phase II, start of consumption) we continuously change Qsets to their values.After that (phase III), the actual change in parameters begins according to the dynamic scenario under study, or the period of relaxation of the system, if the goal is to find a stationary solution.Finally (phase IV), the system reaches stationarity, if entering to this regime is possible.
Scenario dyn2 presents the solution of a dynamic problem with two supplies, of Pset and Qset types.Qset supply is initially inactive, then gradually takes over from Pset supply.The problem can be considered as quasi-stationary, that is, a sequence of stationary problems with changing parameters.The dynamic terms in this setup can be neglected.The purpose of this test is to test the functionality of dynamically changing parameters in equations.On the flow profiles through the 3 pipes in this network, all 4 the described solution phases can be seen, with the actual switching of supplies in phase III.
Scenario dyn3 explores the dynamic opening of valve, which we count as the control element in Table 1.Note that the goal was not to describe accurately the physical phenomena occurring in this case, that is similar to the process of pipe depressurization.Such exact description generally requires more complex modeling [17].The goal was to study the stability of the integrator in response to a sharp change in the equations in the system.The dynamic valve was modeled by the equation ( 7), the state of the valve was described by the parameter α ∈ [0, 1].In the piecewise constant dyn3-pwc scenario, α jumped from 0 to 1. Consequently, the pressure at the valve output jumped, leading to a jump in density through EOS.The mass flow required for this has the form of a delta function, which can be seen in the corresponding image.Also of interest are the oscillations after the main peak shown in the inset.In the piecewise linear dyn3-pwl scenario, α went from 0 to 1 linearly in time.In this case, the main peak was smoothed out, while the oscillations following the peak are similar to the pwc case.
Scenario N1 integrates a realistic network of moderate complexity.The corresponding graphs show 4 phases of the dynamic solution, ending with stationarity in phase IV.Scenarios N85 with realistic networks of high complexity were used to test the stability of the integrator.As a result, 1 scenario out of 85 turned out to be divergent.For comparison, the direct solution of the stationary system by the Newton method under the same conditions gave 6 divergent scenarios out of 85.Thus, the solution of the stationary problem by the integration method is more stable than the direct method, although this stability is not absolute.The possibility of universal adjustment of the integration parameters in order to achieve convergence is the subject of our further research.
In greater detail, when solving a discretized system (2) at each time step, it is possible to solve it with a stabilized Newton method up to the required accuracy.One can also solve a linearized system per time step, which is equivalent to limiting Newton's method to one step and turning off the line search stabilizer.Although the second option is faster, in practice it is much less stable for large realistic networks.When degeneracy occurs, the second option diverges exponentially in time steps, while the first option, at the cost of a large number of iterations, can get out of the degeneracy region.Even if the method hangs or cycles, it has a chance to get out of there at the next time step.Therefore, we choose the first option in our tests.
Infeasible scenario N85.1 is deliberately chosen as such to study the behavior of algorithms on problems of this type.The stationary solution with s -regularization gives an answer in the nonphysical region of negative pressures.Interestingly, the dynamic solution is localized in the region of positive pressures, however, it exceeds the practically admissible boundary P max = 150 bar.In this case, infeasibility manifests itself in a different way.The part of the network adjacent to node n1 goes into stationarity.The other part, adjacent to the node n2205, goes not into stationarity, but into the runaway inflatory mode.This node possesses a maximal speed |dx/dt| taken as a criterion of non-stationarity.A detailed consideration of the solution shows that some control elements turn out to be closed, cutting off a section of the network with an unbalanced sum of Qsets.In such regions, a stationary state is impossible, either inflation or deflation is realized there, depending on the sign of the sum of Qsets.The study of the relation between the occurrence of such regions and the infeasibility of a stationary problem is also in our future plans.

Conclusion
This paper studies the effect of control elements such as compressors, regulators and flaptraps on the stability of fluid transport simulations.It is shown that the modeling of these elements can lead to instabilities, both in stationary and dynamic simulations.Special regularization methods have been developed to overcome these problems, their functionality has been demonstrated in a number of numerical experiments.The methods assign a positive volume to the nodes where only control elements are adjacent, serving as a buffer for dynamic pressure; introduce a kinetic term in equations of control elements, stabilizing the flow balance; choose a large integration step to suppress the oscillations; start solving the problem from zero flow and atmospheric pressure, to satisfy initially all equations in the system.
As numerical experiments show, the integrator withstands different modes, with a small step, with a large step, a sharp change in the initial state or equations (although not recommended), with the appropriate setting -integration of complex networks with a large number of control elements, both for feasible and for infeasible problems.
Our ongoing work consists in extending dynamic modeling to non-trivial fluid composition, temperature dynamics, complex equations of state and phase transitions.The study of the stability of stationary points, the automatic choice of integration parameters and the adaptive choice of the integration step are also in our future plans.

Figure 1 .
Figure 1.Working diagrams for regulator, compressor, flaptrap.Connections leading to undefined pressure and flow.

Figure 4 .
Figure 4. Results of numerical experiments (see text for details).

Table 1 .
Parameters of dynamic test networks.