Nonlinear dynamics and safety aspects of pressure relief valves

Direct spring operated pressure relief valves are regular elements of high-pressure hydraulic systems; however, valve vibrations may still occur endangering the valve and the protected system. The damping of these vibrations is a crucial issue, whereas the external viscous damping at the valve cannot be increased boundlessly based on safety aspects. The vibrations emerge via Hopf bifurcations of either super- or subcritical type, and the nonlinear characteristics have an important role in the valve dynamics. Thorough nonlinear analysis is needed to guarantee safe operation. The linear stability analysis is extended by local and global bifurcation analysis to reveal the stable domains bounded by subcritical Hopf points and the possible bistabilities that exclude the safe operation. Afterwards, the effect of a downstream pipe is analysed, which can also carry safety risks based on the literature; however, the time delay originated in the wave propagation in the pipe may also have stabilising effects. The nonlinear analysis of the presented system compares analytical, semi-analytical and numerical results and assigns the safe parameter domains of the valve with respect to a wide range of relevant parameter variations.


Introduction
Pressure relief valves are safety critical elements of high-pressure hydraulic and pneumatic systems. The essential elements of a direct spring operated pressure relief valve are a precompressed helical spring and a valve disc (see Fig. 1). In normal operation of the hydraulic system, the valve disc is pressed against the valve seat, so it closes the space of the pressurised fluid. The pressure level of the valve opening is called set pressure, and it can be adjusted with the precompression of the spring. The use of external viscous damping on the valve is optional, while it is forbidden by some standards in order to guarantee the fast valve opening. These valves open only in emergency events of the protected hydraulic system. Because of the possibly rare opening, the dashpot aiming to provide additional damping can get stuck, which could have catastrophic consequences; external viscous damping carries high safety risks and maintenance requirements. Still, the magnitude of damping becomes relevant when valve vibrations occur [1,2]. These vibrations can be detected by the harsh noise which often includes the sound of the impacts between the valve disc and the valve seat. The phenomenon can lead to mechanical damage of the valve or to insufficient discharging process.
These undesired vibrations are related to the loss of stability of the equilibrium state of the open valve. As it is already known from the literature [3,4], the loss of stability always takes place via Hopf bifurcations that could be either super-or subcritical. The global dynamics of the valve is strongly determined by the type of the Hopf bifurcation. The supercritical case was thoroughly studied in the literature [5]; the emerging stable oscillations can be so large that impacts occur between the valve disc and the valve seat, which is also called chatter. The corresponding grazing bifurcations and the cascade of period doubling bifurcations leading to chaotic oscillations were also identified [5] in the corresponding parameter regions. While the stable oscillation of the valve should be avoided, the impacts occurring after the grazing bifurcations are definitely not permissible. The study in [5] also indicates the sensitivity of the type of the Hopf bifurcation to parameter changes; however, there is no detailed analysis of the relevant system parameters on how they affect the stable operation of the valve, moreover, how they influence the type of the Hopf bifurcations.
While the possibility of the subcritical Hopf bifurcation is mentioned in the literature [5,6], the corresponding global dynamics of the valve has not been explored. In these cases, the hysteresis phenomenon can occur if, for example, the flow rate is increased and decreased, and the stable equilibrium state of the open valve may coexist with stable chatter. Consequently, the clarification of the effects of the corresponding initial conditions are relevant from the viewpoint of the operation of the valve. This requires the analysis of the full nonlinear system, since the linear stability analysis has limited practical relevance.
Further complication in practice is the industrial demand to categorise a valve whether its steady state is a stable equilibrium or a stable chatter. However, the environment of the valve has an essential effect on the steady states, so the valve behaviour can not be predicted without modelling the hydraulic elements (pipes, vessels) surrounding the valve. Increasing vessel volume can have a stabilising effect [7], while long upstream pipes usually destabilise the equilibrium [3]. This was also confirmed by the experimental work in [8]. Despite the common engineering approach that the pipes connected to the valve have a negative effect, the downstream pipe may stabilise the valve equilibrium [7,9].
In [7], a mathematical model is derived for a system consisting of a vessel with constant inlet fluid flow rate, a direct spring operated pressure relief valve, and a downstream pipe. In the limit case of the zero pipe length, the vessel-valve subsystem is a threedimensional system of nonlinear ordinary differential equations (ODE). The presence of the pipe extends the model with two partial differential equations (PDE), which can be transformed to delay differential algebraic equations (DDAEs). It was pointed out for given inlet flow rates, and for relatively low damping levels that in the parameter plane of the pipe length and vessel volume, stable and unstable regions follow each other with a certain periodicity of the pipe length, till an upper pipe length limit above which the system is always unstable.
In this study, the model presented in [7] is analysed further from nonlinear vibrations viewpoint. Section 2 briefly summarises the modelling of the vessel-valve system. The subsequent sections of the paper present the stability investigation and the nonlinear analysis of the system. Finally, Sect. 5 deals with the nonlinear dynamic effects of the downstream pipe; the conclusions are summarised in Sect. 6.

Modelling
The mechanical model in question is shown in Fig. 1. The valve is modelled with the modal mass m of the moving parts, the spring stiffness k and damping coefficient c. The angular natural frequency of the valve and the damping ratio are respectively. The set pressure p set of the valve is proportional to the spring precompression z 0 . The closed valve opens if the fluid force acting on the valve disc through the outflow orifice cross-section A becomes equal to the spring force. The mathematical formulation of this operation rule leads to the following expression for the set pressure: Custom-designed pressure relief valves can be chosen from catalogues for different operation conditions considering mainly the set pressure and operating temperature ranges. The modal mass and the system stiffness are determined by the manufacturer that means the angular natural frequency is given, and the set pressure is prescribed by the desired pressure level in the protected system. In this context, the damping ratio primarily depends on the damping coefficient c, which has to be estimated during modelling. The vessel of volume V is charged by a constant inlet flow rate Q in . The fluid has its density ρ and bulk modulus K . The state variables of the vessel-valve system are the valve lift z(t), the valve disc velocity v(t) = dz/ dt and the gauge pressure p(t) in the vessel. The valve can be modelled as a one degree of freedom oscillator where the force F(t) acting on the valve disc depends on the pressure drop Δ p at the valve: In case of the vessel-valve system, the valve discharges the medium directly to the atmospheric pressure p 0 corresponding to zero gauge pressure, which means that the pressure drop at the valve is Δ p = p. The differential equation describing the pressure variation in the vessel can be derived from the mass balance equation: where the outflow Q out through the valve can be calculated based on the Bernoulli Law: The outlfow orifice of the valve is a cylindrical surface of area Dπ z(t) depending on the valve lift z. The expression √ 2Δ p(t)/ρ is the fluid velocity in the orifice. The discharge coefficient C d is an empirical factor, which can be measured or calculated by CFD simulations [10,11]. It is assumed to be constant during the calculations.
Thus, the model of the vessel-valve system is a threedimensional nonlinear system of ordinary differential equations. The dimensionless time T and the new state variables y 1,2,3 are introduced in the following way: where the reference quantities are: With these state variables, the dimensionless system of ordinary differential equations can be derived: where dot refers to the derivative with respect to the dimensionless time T , and the dimensionless parameters are expressed in the form of The dimensionless flow rate q is proportional to the constant inlet flow rate Q in into the vessel. The dimensionless set pressure δ is proportional to the set pressure adjusted with the spring precompression. The dimensionless system stiffness β is proportional to the bulk modulus K of the fluid and inversely proportional to the vessel volume V . During the analysis, the vessel volume is associated with this dimensionless parameter. In case of an infinitely large vessel (β → 0), the rate of change of the vessel pressure p is zero, see (5) and (27). This leads to "infinite softness" of the hydraulic system before the valve. The infinitely small vessel (β → ∞) means direct delivery of the fluid. This is the case of "infinit stiffness" of the hydraulic system. Note that this way, a solid vessel is considered, which does not influence the hydraulic stiffness. The elasticity of the walls or other elastic elements added to the vessel can have influence on the parameter β, which are now neglected for the sake of simplicity, while they could be used to manipulate the system stiffness in a useful way.

Linear stability analysis
The equilibrium of the vessel-valve system is at the trivial solution y 1,2,3 (T ) ≡ y * 1,2,3 of the mathematical model, which can be calculated from the equations in a flow rate dependent form. It has only one positive solution for the equilibrium valve lift, that is, y * 1 ≥ 0. The dimensionless setpressure δ > 0 is considered as a fixed parameter because it is usually not varied continuously. In the followings, δ = 3 is considered, which means that the set pressure is 3 bar gauge pressure compared to 1 bar atmospheric pressure. It is known from the literature that increasing set pressure has negative effect on valve stability [5,7].
The perturbation vector η = [η 1 η 2 η 3 ] T can be introduced around the equilibrium as y(T ) = y * + η(T ). The linearised system at the equilibrium assumes the form: where J is the Jacobian matrix: The corresponding characteristic equation assumes the form The linear stability analysis is carried out by means of the Routh-Hurwitz criterion. According to the physical meaning of the parameters in case of an open valve, the following conditions are considered: β, ζ, y * 1 , y * 3 > 0, which guarantees that all the coefficients are positive in (18). The criterion of the second Hurwitz determinant has to be fulfilled in order to achieve stability. The limit of dynamic loss of stability is at H 2 (q, β, ζ ) = 0, which has one positive root for the damping ratio: It presents the lower limit of stability for a given value of the inlet flow rate q and system stiffness β based on the stability criterion H 2 > 0. The ζ min (q, β) defines a surface in the parameter space of the inlet flow rate represented by q, the vessel volume given by β and the damping ζ of the valve. This boundary surface of the dynamic loss of stability is shown in Fig. 2a. Note that in case of an undamped valve, stable operation is not possible, because ζ = 0 substitution in (19) leads to the condition of stability −β y * Two damping levels are represented in Fig. 2a by horizontal planes, which have intersection lines with the boundary surface. These intersection lines are the boundary curves in the stability charts constructed for fixed values of the damping ratio ζ . This way, the stability charts give information about possibilities of intervention, because the size of the vessel can be designed in a proper way and the inlet flow rate can be regulated during operation. The damping ratio is hardly tunable. The damping consists of the possible external damping and the damping effect of the oiled and/or moistened surfaces sliding on each other during valve vibrations. Accurate measurements are necessary for estimating the actual level of damping at the valve.
The stability charts for the two different damping levels in Fig. 2a are presented in panel (b). For a strongly damped valve (ζ = 0.6), the stable domain is wide, and a turning point of the boundary curve appears inside the investigated parameter window. However, the slightly damped valve (ζ = 0.2) has only a narrow stable domain for small values of β, that is, for large vessel volumes V . For clear visualisation, a logarithmic scale of β is applied in the second case. The turning point of the stability boundary is out of the practically relevant parameter window (for q > 100), so the stability chart suggests that practically, there are two separated stable domains. The upper stable domain Because of the uncertainty of the damping, the minimal damping ζ min has an important role for a given vessel size and operation flow rate. The equilibrium can be stable only for damping values larger than ζ min , which implies the question how large external damping is feasible. In the followings, the ζ = 1 is considered as an upper limit of the feasible damping ratio. Figure 3 shows a map for the minimal damping for different values of the dimensionless inlet flow rate q and system stiffness β. The damping values are saturated at 1; thus, white colour refers to the domains with minimal damping ratios larger than or equal to 1.
The damping levels in the range ζ < 0.4 might be feasible without additional external dashpot. The location of the slightly damped stable domains, for example, the one bounded by ζ = 0.3, is shown in Fig. 3.
The results show that β → 0 (V → ∞) has a stabilising effect (see the stable domain S V ∞ in Fig. 2b); for any damping, there exists a vessel volume above which the stability is flow rate independent. On the other hand, for relatively large β, and for relatively large inlet flow rates q, another beneficial stable domain appears. (See the stable domain S V 0 in Fig. 2b.) Thus, both the small and large vessels can have a stabilising effect. The gradient of the surface in Fig. 2a implies that large inlet flow rates also stabilise the system.
Design aspects serve limitations on the utilisation of these stabilising engineering solutions. Large vessel cannot be an ultimate solution for valve vibration issues because the size of the protected hydraulic circuit and the space in the engine may restrict the vessel size. The valve must operate around its nominal flow rate determined by the set pressure and the outflow orifice of the vessel, so the flow rate cannot be increased boundlessly either.
Although the above linear stability analysis presents a guide for the design of these systems, it does not provide information about the basins of attraction of the stable valve operations. Detailed nonlinear analysis is needed, for example, to explore the sensitivity for perturbations in the initial conditions.

Hopf bifurcation calculations
The linear analysis in Sect. 2 pointed out that the damping not only influences the decay of the transient vibrations, but also determines the contour levels in Fig. 3, which are stability boundaries corresponding to the dynamic loss of stability. The actual self-excited vibrations can emerge through sub-or supercritical Hopf bifurcations depending on the system parameters.
In engineering sense, the loss of stability is safer through supercritical Hopf bifurcation than through a subcritical one. In supercritical case, the equilibrium becomes unstable, and stable periodic orbits grow in the neighbourhood of the unstable equilibrium. In case of the subcritical Hopf bifurcation, an unstable limit cycle exists around the stable equilibrium hedging in the basin of attraction of the stable equilibrium. This dynamical behaviour usually leads to coexisting stable solutions in the bistable regions resulting in hysteresis phenomena, which cannot be predicted by the linear analysis.

Methods, algorithms
The bistable region means that for certain parameter combinations, the equilibrium position of the open valve is stable, while a stable oscillation also exists for large perturbations in the initial conditions. This scenario must be avoided in case of a safety device. The bistable regions can be revealed by constructing the bifurcation diagram. The inlet flow rate q is an appropriate choice for the bifurcation parameter because it can be modified during operation.
With the help of the vibration amplitudes, the location of grazing bifurcation can be predicted where the valve touches the valve seat. The grazing bifurcation separates the regions of the chatter and the vibrations without impacts. In case of the unstable limit cycles, the basin of attraction of the stable equilibrium is also defined by the radius of the limit cycle. The type of the Hopf bifurcations and the amplitudes of the emerging self-excited vibrations can be determined with the analysis of the nonlinear system.
Analytically, the type of the Hopf bifurcation and the radius of the limit cycle are predicted by means of the Hopf bifurcation calculations that include the centre manifold reduction. The details of the calculations can be found in "Appendix A". The sign of the so-called Poincaré-Lyapunov coefficient Δ indicates the type of the Hopf bifurcation: Δ > 0 corresponds to the subcritical cases, Δ < 0 to the supercritical ones. The radius of the limit cycle is approximated by the following formula: The appearing Re(λ )| q=q cr derivative of the characteristic root with respect to the bifurcation parameter q can be calculated by the implicit differentiation of the characteristic equation P(λ) = 0 in (18). This nontrivial step includes the derivation of the inlet flow rate dependent equilibrium state y * 1 (q), 0, y * 3 (q) given in (15). Note that the nonlinearities are strongly asymmetric due to the square root appearing in (27). The asymmetric nonlinearities alter the shape of the emerging limit cycles compared to the standard symmetric parabolic approximation.
The DDE-Biftool MATLAB package is a welltested tool using continuation method for tracing the equilibrium state and the limit cycles of a time delay system. The mathematical model of the vessel-valve system does not contain time delay, so it is set to be zero in the inputs. The critical point and the corresponding Poincaré-Lyapunov constant can be calculated and compared to the analytical results. The software package allows to follow the possible folds in the shape of the limit cycle branches depending on the bifurcation parameter q, which is not possible with the analytical Hopf calculations in case of third-order approximations; fifth-order power series would be needed for the analytical prediction of the folds in the bifurcation diagrams.
The coexisting stable solutions can be explored by analysing the steady-state solutions of the nonlinear system. The time series of the valve lift can be produced by an appropriate Runge-Kutta (RK) method, for example, the ODE45 built-in solver of MATLAB. The event handler of MATLAB can be used to detect the time instant t imp of the impact between the valve disc and the valve seat. In the case when the condition y 1 = 0 fulfils, the simulation stops and the new initial conditions are served by the Newtonian impact law: where t − imp and t + imp are the time instants before and after the impact, respectively, and e = 0.9 is considered for the coefficient of restitution. In case of chatter, this method is applied at each impact. The continuation method and the analytical calculations do not handle this discontinuity, but they are able to predict the location of the grazing bifurcation assigning the starting point of the region of the impacting oscillations. For this relatively high value of the coefficient of restitution, similar chaotic oscillations can be observed as in [5].

Local nonlinear analysis
The nonlinear analysis is carried out by combining the above analytical, semi-analytical and numerical methods. The type of the Hopf bifurcations is determined in the surface shown in Fig. 2, which presents an overall picture in the parameter space about the local nonlinear dynamic behaviour. Along the stability boundaries in Fig. 3, the unsafe sections will be selected by the subcritical Hopf bifurcations. The global dynamics is to be explored by constructing bifurcation diagrams.
First, the type of the Hopf bifurcation is determined with the help of the DDE Biftool on a (q, β) mesh. At the mesh points, the minimal damping ratios ζ min (q, β) are known from (20). The inputs for the DDE-Biftool are the system stiffness values and minimal damping ratios at these mesh points, where the bifurcation parameter is the inlet flow rate q. The bifurcation point  Fig. 4. Red plus symbols refer to the subcritical Hopf bifurcation points and blue ones to the supercritical points. The border line separating the domains of the different bifurcation types, which consists of the Bautin points [12], does not coincide with any contour level. While the supercritical parameter domain is much larger than the subcritical one, the practically relevant low minimal damping levels fall in the subcritical domain.
The contour level of ζ = 0.39 runs close to the limit between the sub-and supercritical domains, but it clearly intersects the line of the Bautin points. Along this contour level, the Poincaré-Lyapunov constants were checked analytically. Red dots mark the subcritical cases and blue dots the supercritical points. The Bautin point along the contour level ζ = 0.39 is at β ≈ 9.7, q ≈ 5.4. The analytically calculated critical vibration frequencies along this contour level are also shown in Fig. 4 depending on the inlet flow rate q. These frequencies correspond to the critical complex conjugate roots ±iω of the characteristic polynomial (18). According to the definition of the dimensionless time T in (7), the real vibration frequencies are = ωω n , where ω n is the angular natural frequency of the valve.

Global bifurcation diagrams
The numerical bifurcation diagrams are produced with the help of the MATLAB ODE45 solver. Parameter sweep is carried out for increasing and decreasing series of the bifurcation parameter q and the steady-state solutions are represented with the equilibrium valve lift and the maximal and minimal displacements of the valve. At each parameter step, the slightly perturbed steady state of the previous step serves as the initial conditions of the simulation. For the parameter step size Δq = 0.2, 15% perturbation was applied in the dimensionless valve lift y 1 compared to the steady state of the previous step. The simulation time in each parameter step is relatively long: ΔT = 1000. This way, the stable solutions are followed and the bistable regions are identified in the numerical bifurcation diagram.
In Figs. 5 and 6, examples are given for the different types of the Hopf bifurcations. The parameter combinations were selected from Fig. 4 at the contour level ζ = 0.39. The local minima and maxima of the steady-state times series are denoted by dots and circles, respectively. Blue refers to the supercritical cases and red to the subcritical cases. The black dashed line represents the analytically determined equilibrium states. The branches of the stable and unstable periodic orbits emerging through the Hopf bifurcations are calculated both analytically and by continuation. The analytical results are produced by transforming back the radius of the limit cycle given by the Hopf calculations. Figure 5 shows a bifurcation diagram at q = 5 for the lower value of the system stiffness β ≈ 0.41. This is the safest case from engineering point of view, because by increasing the inlet flow rate q, the equilibrium is stable until the bifurcation point, after which the amplitudes increase relatively slowly, without producing immediate impacts with the valve seat. These vibrations correspond to stable limit cycles and the amplitudes calculated by the three applied methods are in good correspondence. The result of the Hopf calculations is shown in yellow dotted line, and the stable limit cycles calcu- lated by continuation is represented in cyan continuous line. Both the analytical limit cycle calculation and the continuation method can be applied to predict the grazing bifurcation, which is the lower flow rate limit of the chatter.
A subcritical case is represented with the dimensionless stiffness value β = 12 in Fig. 6a and another supercritical case is shown for β = 6 in panel (b). In both cases, hysteresis appears, however, at β = 6, there is supercritical Hopf bifurcation (see also Fig. 4). In the subcritical case, the unstable limit cycle calculated by continuation method is shown in light red dashed line fitting perfectly to the numerically produced stable solutions. The analyitically predicted unstable limit cycle does not meet the result of the continuation method well because of the asymmetric nonlinearities in the system, but it still shows the wide range of coexisting stable solutions.
The subcritical behaviour is dominant in the neighbourhood of the Bautin point; as Fig. 6b shows, a fold bifurcation changes the shape of the limit cycle branch in the neighbourhood of the supercritical Hopf point leading to a bistable region. This fold bifurcation can be detected with the help of the continuation method. The analytically predicted limit cycle is valid only locally in a narrow vicinity of the bifurcation point, but it cannot be applied for the prediction of the grazing bifurcation, which assigns the upper flow rate limit of the chatter. Panel (b) illustrates how the third-order approximation Consequently, the local bifurcation analysis visualised in Fig. 4 is still not satisfactory to predict the unsafe domains; the global behaviour has to be taken into account to predict the bistable regions in order to guarantee the safe operation of the valve.
Although supercritical Hopf bifurcation (see Fig. 4) can occur at feasible low damping values and small vessel volumes (see the domain S V 0 if Fig. 2b), it is an important observation that these domains are covered by bistable zones leading to unsafe valve functioning.
In case of the pressure relief valves, it must be emphasised that the emergency case of the hydraulic system prescribes the initial conditions for the valve dynamics: the valve is closed with zero disc velocity when the increasing pressure level just reaches the set pressure. This means that the initial conditions are y 1 (0) = 0, y 2 (0) = 0, y 3 (0) = δ, which are outside the domain of attraction of the stable equilibrium of the valve in the bistable zone (see Fig. 6); thus, the large amplitude oscillations are unavoidable. In this respect, this situation is different from other dynamical systems with subcritical Hopf bifurcations like the shimmy problem of rolling wheels [13], the machine tool vibrations [14] or the human balancing [15], where the steady states of the normal operation are the equilibrium states, and the stable large amplitude vibrations occur only for large enough random perturbations. In case of the pressure relief valve, the valve opening determines the perturbation so the bistable regions are generally not utilisable.
The other supercritical boundaries for low values of β in Fig. 4 (see also S V ∞ in Fig. 2) correspond to a classical supercritical Hopf bifurcation without fold bifurcations and bistability. Although these stable domains are safe here, the values of the system stiffness refer to large vessel volumes V → ∞. If the enlargement of the vessel volume is possible, it serves as a rule of thumb for the design of the hydraulic system. The large vessel cannot be an ultimate solution for the chatter problem of valves because of the limited space in many hydraulic circuits. All in all, we have to compromise between the need of low external damping level and the safe operation in a wide range of the parameters.

Downstream pipe as a damper
The model presented in Sect. 1 can be extended for the case of attached downstream pipes by means of the Navier-Stokes and continuity equations. These two first-order, linear partial differential equations (PDE) have the distributed state variables of the gauge pressurep(x, t) along the pipe and the fluid velocityṽ(x, t) in the pipe; the space coordinate along the pipe is x ∈ (0, l), where l is the pipe length. There are two boundary conditions expressing that the mass flow through the valve must be equal to the inlet mass flow into the pipe, while on the other side, the pipe delivers to the atmospheric pressure. The connection between the ODEs and PDEs appears in the pressure drop Δ p at the valve, which changes to Δ p(t) = p(t) −p(0, t) compared to p(t) introduced in Sect. 1. The PDEs can be transformed into a dimensionless form where the state variablesṽ andp are related to y 4 (X, T ) and y 5 (X, T ), respectively. Here, T is the dimensionless time similarly to the one in (7), and X = x/l is the dimensionless space coordinate along the pipe. The detailed derivation of the system of coupled ODEs and PDEs is given in [7] with the final form: By means of the travelling wave solutions, these partial differential equations can be transformed to delay differential algebraic equations (DDAEs), where the delayed state variable is denoted by f that gives the dimensionless gauge pressure y 4 and also the dimensionless flow velocity y 5 in the pipe [7]: The dimensionless time delay τ = 2lω n /a is proportional to the pipe length; it is the dimensionless time needed for the wave propagation in the pipe to and fro with the speed of sound a. The resulting nonlinear DDAE model consists of the following equations: Besides the time delay τ , another new parameter is the dimensionless pipe inlet parameter φ that is inversely proportional to the pipe cross-section area A p : Note that various terminologies appear in the literature for such kind of systems that involve a delay algebraic equation (see [16]) similar to (37), like shift map [17] or renewal equation [18,19]. According to the latter terminology, the algebraic or difference equation in (37) is a neutral renewal equation.
The linear stability analysis at the equilibrium results in closed form vibration frequencies of the emerging self-excited vibrations [7]: which is delay-independent. Since these values do not depend on the pipe length, (39) presents the same vibration frequencies which were calculated from the pure imaginary zeros of the characteristic polynomial in (18) for the ODE system describing the vessel-valve subsystem only. The critical vibration frequencies for the system without and with downstream pipe are presented together in Fig. 7 against the inlet flow rate q.
The stability boundaries are also given in closed form expressed for the critical system stiffness β cr (τ ; ω) (see the formula (56) in "Appendix B"), where ω must be substituted as ω 1 or ω 2 leading to two types of stability losses in the plane of the system stiffness β and time delay τ . The stability charts in this plane characterise the environment of the valve because the time delay includes the downstream pipe length, and the system stiffness β refers to the vessel volume V . In case of various flow rates, typical stability charts are shown in Fig. 8 for the damping level ζ = 0.39. As the panels in Fig. 8 show, a wide stable domain appears for increasing time delay with an upper system stiffness boundary of stability. These stable domains have an ultimate upper time delay limit corresponding to an asymptote of the expression for the stability boundary in (56) or at an intersection point of the two types of the stability boundaries. In case of τ → 0 at a given inlet flow rate q, the critical system stiffnesses β cr,1 and β cr,2 converge to the points of the contour level ζ = 0.39 in Fig. 4. This means that the lower stability boundaries corresponding to the safe supercritical cases illustrated in Fig. 5 can be shifted in the direction of increasing β values, so the beneficial S V ∞ stable region (see Figs. 2b and 4) can be enlarged. The results for increasing time delays are shown in Fig. 9, where the stability boundaries for the different time delays are upper limits of stability in the S V ∞ domain. Thus, for increasing time delays, this stable domain increases significantly.
The identificaton of the type of the Hopf bifurcation is not straightforward in the delayed model because of the implicit form of the DDAEs in (34-37). Development of numerical methods dealing with renewal equa-  [20]. The dots on the stability boundaries refer to the predicted type of the Hopf bifurcation based on (40-43) and Fig. 10. Blue refers to super-, red to subcritical cases, while filling refers to the local bifurcation type and the edgecolour to the global characteristics due to existing fold bifurcations. The white dot shows the location of the Bautin point. (Color figure online) tions is in progress [21], but it is not yet directly applicable for the system of (34-37). However, an artificial stiffness S can be introduced in order to allow the application of the continuation method by means of the DDE Biftool package. The original implicit system of DDAEs is replaced with approximate DDEs in the following way: When S → ∞, the above system becomes equivalent to (34-37). Theoretically, the value of S must be as high as possible, although in case of numerical calculations, the cumulation of the numerical errors has to be taken into account: it causes the breakdown of the code limiting the maximal applicable value of this artificial stiffness.
The continuation method works with the bifurcation parameter choice of the system stiffness β instead of the inlet flow rate q, presumably because the equilibrium state is independent of β, while the inlet flow rate q influences the equilibrium state step by step during continuations. Still, the local bifurcation analysis determines the type of the Hopf bifurcation independently of the parameter choice [12]; thus, the subsequent results can be compared to the results in Fig. 4. By increasing the artificial stiffness level S, fast convergence can be observed for the vibration frequencies ω and the critical bifurcation parameters β cr towards their analytical values. On the other hand, the Poincaré-Lyapunov constant is highly sensitive to the values of the artificial stiffness S and to the starting point of the continuation, that is, its sign can alternate even when the vibration frequencies and the critical bifurcation parameters reached their exact values. Consequently, rather the graphical results (like the ones in Fig. 10) of the continuation method are used to characterise the local and global nonlinear behaviour.
The continuation method was carried out at some chosen critical points of the stability boundaries in Fig. 9, where blue dots show the supercritical cases, red dots the subcritical ones, while the Bautin point is marked by white dot. In some cases, despite the locally subcritical bifurcation, the global characteristics is supercritical due to a fold bifurcation of the unstable limit cycle branch located close to the Hopf point (see Fig. 10b). The edgecolour of the dots in Fig. 9 refers to the type of the global behaviour. The results show that the stability boundary of the beneficial stable domain S V ∞ remains supercritical up to an upper limit of the inlet flow rate q while shifting upward with increasing τ ; this means that the extended stable S V ∞ domain still has a practically safe region.

Conclusions
A simplified model of a pressure relief valve and its environment was presented and analysed in order to explore the complex valve dynamics in view of the parameters describing the valve and its environment. The main goal is to guarantee safety, which means fast valve opening without vibrations. Due to the demand of the immediate opening in an emergency case, the level of the allowable external damping is limited considering the outlined disadvantages of a dashpot with the possibility of getting stuck.
The undesired vibrations of the valve emerge through Hopf bifurcations and correspond to stable periodic solutions. The linear stability analysis shows the locations of the Hopf bifurcation and the stable regions of the equilibrium states, while the nonlinear analysis identifies the safe parameter domains. In case of subcritical Hopf bifurcations, the basin of attraction is usually restricted by the radius of the unstable limit cycle causing bistability. In addition, the initial conditions are prescribed by the emergency case of the protected hydraulic system; the valve is closed until the vessel gauge pressure reaches the set pressure. Based on this prescribed large initial perturbation, the stable domains bounded by subcritical Hopf bifurcations are pronounced unsafe, while the ones with supercritical boundaries are associated to safe valve operation.
The global dynamics, however, has to be also checked in the neighbourhood of the supercritical Hopf bifurcations since fold bifurcations can change the shape and stability of the limit cycle branches leading to bistable behaviour similarly to the locally subcritical cases. The nonlinear analysis compares three different approaches: analytical Hopf calculations, semianalytical investigation by DDE Biftool using continuation method, and numerical integration by Runge-Kutta method considering also impacts between the valve disc and valve seat. The results of these analyses of the vessel-valve system can be summarised in the following steps: -The linear stability analysis was carried out in the parameter space of the dimensionless inlet flow rate q, dimensionless system stiffness β and damping ratio ζ . As a conclusion, a lower damping limit ζ min exists for each (q, β) combination, above which the equilibrium is stable; the ζ min (q, β) boundary surface represents the limit of the dynamic loss of stability. The stable parameter domains of the low damping levels are restricted for two regions: one for large vessels (β → 0), and another for relatively small vessels (that means relatively large β) with large inlet flow rates. -The results of the local bifurcation analyses prove that a wide range of the critical points are supercritical including the β → 0 stable domain; however, in the linearly favourable cases with low minimal damping levels for small vessel volumes and large flow rates, the stable domains are mostly bounded by subcritical Hopf bifurcations. The line consisting of the so-called Bautin points is also identified, which separates the subcritical and supercritical parameter regions. -The global dynamics was explored by constructing bifurcation diagrams with respect to the bifurcation parameter of the inlet flow rate q. The boundary of the beneficial stable domains for small minimal damping levels and for large vessels (β → 0) has the safe supercritical characteristics. For small minimal damping levels and for relatively small vessels, the global nonlinear analysis showed subcritical characteristics even in the locally supercritical cases. In brief, for the vessel-valve system, only one safe domain exists: the case of large vessels (β → 0).
The assigned safe region for β → 0 cannot be utilised if the necessary vessel volume is larger than the available space in the protected hydraulic system. Nevertheless, the environment of the vessel may provide an engineering solution to achieve stability and safe operation for smaller vessels. The presence of a downstream pipe involves the time delay τ in the mathematical model originated in wave propagation; thus, it is proportional to the pipe length l. The main conclusions regarding the vessel-valve-pipe system are: -In case of an appropriately chosen pipe length, the time delay has a powerful stabilising effect, which enlarges the single safe stable domain in the (q, β) plane of the vessel-valve system. Stability charts of the delayed cases were presented to compare these extended stable domains to the delay-free case. -Global nonlinear analysis was carried out via continuation method of DDE Bibftool. The β cr boundary is supercritical for τ = 0, and the investigated stabilising τ > 0 cases inherit this property up to an upper inlet flow rate limit corresponding to the Bautin point. Still, a wide safe stable domain exists which can be utilised instead of increasing the vessel size or using external viscous damper. -The delayed model of the vessel-valve-pipe system is a system of implicit DDEAs. An approximate system of DDEs was presented by introducing an artificial stiffness, which makes the continuation method possible by means of DDE Biftool. The vibration frequencies and the Hopf point converge fast towards their exact values by increasing the artificial stiffness values, and the limit cycle branches can be followed graphically. The applicable stiffness values were bounded by the breakdown of the code, while the equilibrium state had to be independent of the bifurcation parameter. The possible implementation of DDE Biftool requires further study for nonlinear implicit DDAEs.
In summary, the design of pressure relief valves is not possible without the investigation of its global nonlinear dynamics including the hydraulic environment surrounding the valve.
Author contributions Both authors contributed to the study conception and design. Material preparation, data collection and analysis were performed by Fanni Kadar. The first draft of the manuscript was written by Fanni Kadar, and Gabor Stepan corrected previous versions of the manuscript. Both authors read and approved the final manuscript.
Funding Open access funding provided by Budapest University of Technology and Economics. The research reported in this paper has been supported by the Hungarian National Research, Development and Innovation Office, Grant Numbers NKFI-KKP-133846, NKFI-K-132477 and NKFI-FK-134496.

Data Availability Statement
The manuscript has no associated data.

Conflict of interest
The authors have no relevant financial or non-financial interests to disclose.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/ by/4.0/.

Appendix A
The Hopf calculations for the three-dimensional system of ordinary differential equations (25-27) are carried out at the critical point q cr , where the characteristic Eq. (18) has one complex conjugate pair of roots λ 1,2 = ±iω and a negative root λ 3 < 0. These are the eigenvalues of the matrix J in (17), and the corresponding eigenvectors are the complex conjugate s 1 = s 2 ∈ C 3 , and the real s 3 ∈ R 3 , respectively. The original system has to be expanded in the forṁ η ∼ = Jη + g(η), where g(η) contains the nonlinearities until third order. The origin of the nonlinearity is the square root of the dimensionless vessel pressure multiplied by the dimensionless valve lift (see 27); thus, the nonlinearity vector g has nonzero element only in its third coordinate g 3 : (45) First, the system has to be transformed to the eigenbasis to reach the first-order normal form. The corresponding transformation matrix T is: In case of the vessel-valve model, the transformation matrix is based on the following general form of the eigenvectors: where u is the new vector of the state variables defined by η = Tu, and T −1 JT serves the Jordan block form at the critical point: where the vector of the nonlinearities is denoted by G(u). The next step is the centre manifold reduction in order to reduce the number of dimensions to two, where the Hopf theorem can be applied directly. The invariant centre manifold h(u 1 , u 2 ) is approximated with second degree terms: If the initial perturbation is outside this manifold, the trajectories tend to the centre manifold 'along' the direction of the third eigenvector s 3 . Formula (51) of the centre manifold is substituted back to (50) as u 3 = h(u 1 , u 2 ): where the coefficients of the second degree terms of G 3 (u) will be denoted by c 20 , c 11 , c 02 . For the comparison of the coefficients in the polynomial balance in (52), the derivatives u 1 , u 2 have to be substituted up to first degree: u 1 ∼ = ωu 2 , u 2 ∼ = −ωu 1 . This leads to the linear nonhomogeneous system of algebraic equations which has a unique solution for h 20 , h 11 , h 02 . By substituting this formula of the centre manifold into (50) as u 3 = h(u 1 , u 2 ), the system is reduced to a twodimensional system of ordinary differential equations with the following notations: which is used in the determination of the Hopf bifurcation and the vibration amplitudes of the emerging stable/unstable periodic motions.