Asymptotic analysis of self-excited and forced vibrations of a self-regulating pressure control valve

Pressure vibrations in hydraulic systems are a widespread problem and can be caused by external excitation or self-exciting mechanisms. Although vibrations cannot be completely avoided in most cases, at least their frequencies must be known in order to prevent resonant excitation of adjacent components. While external excitation frequencies are known in most cases, the estimation of self-excited vibration amplitudes and frequencies is often difficult. Usually, numerical studies have to be executed in order to elaborate parameter influences, which is computationally expensive. The same holds true for the prediction of forced oscillation amplitudes. This contribution proposes asymptotic approximations of forced and self-excited oscillations in a simple hydraulic circuit consisting of a pump, an ideal consumer and a pressure control valve. Two excitation mechanisms of practical interest, namely pump pulsations (forced vibrations) and valve instability (self-excited vibrations), are analyzed. The system dynamics are described by a singularly perturbed third-order differential equation. By separating slow and fast variables in the system without external excitation, a first-order approximation of the slow manifold is computed. The flow on the slow manifold is approximated by an averaging procedure, whose piecewise defined zero-order solution maps the valve’s switching property. A modification of the procedure allows for the asymptotic approximation of the system’s forced response to an external excitation. The approximate solutions are validated within a realistic parameter range by comparison with numerical solutions of the full system equations.


Introduction
Hydraulic pressure control valves are widely used in hydraulic transmission systems. Their application represents the simplest way of maintaining a desired pressure level within a certain part of the overall hydraulic circuit independently of the current power demand of the consumer. This is a standard feature of many hydraulic circuits and enables to feed several consumers out of one sub-circuit at the same time. These circuits are also called constant-pressure systems.
Pressure control valves can be designed both as poppet valves or as spool valves. The pressure regulation is realized via an internal pressure feedback to a control surface of the valve piston, providing a feedback force proportional to the system pressure. On the other hand, the piston movement changes a variable control gap, leading to a piston position dependent oil flow out of the system. The pressure level is adjusted by applying a control force on the piston, acting against the pressure feedback force. This force can either be mechanic (typical for pressure relief valves or pressure reducing valves), hydraulic (pilot-operated pressure control valve) or magnetic (proportional valve). While the first variant merely enables one individual unchangeable set-pressure level, the other two valve types, capable of providing a variable control force, allow for more advanced applications, for example, in electro-hydraulic (e.g., brake-or clutch-) actuation mechanisms. From the above elucidations, it follows that pressure control valves are closed-loop control systems [18]. It can even be shown that the included pressure feedback has certain optimality properties in the sense of linear quadratic regulation problems [11]. However, like many closed-loop control systems, it is long known that pressure control valves exhibit stability problems.
Early discussions addressing instability phenomena in hydraulic control valve circuits have been published in the 1960s. In [7], the interaction between pipeline dynamics and valve piston is identified as a possible source of instability. In [19], the internal pressure feedback loop is detected to possibly destabilize the valve's equilibrium position and qualitative design recommendations are given to ensure stable operation conditions. A review of these and other early studies can be found in [8]. Recent investigations on local stability of pressure control valves are carried out analytically in [12,13,16] and experimentally in [1]. It is unanimously discovered that, e.g., high system damping or supply volume flow into the system stabilizes the valve, in case of both poppet and spool valves. Furthermore, local stability strongly depends on the valve's flow characteristic and hence on control edge geometry and leakage behavior.
Since these are purely local properties and given that nearly all hydraulic systems are inherently nonlinear, the question arises to what extent the investigation of local stability is capable of drawing conclusions about the valves suitability in practical use. In this context, two main questions arise: The first one is about the structure of the phase space, namely whether unstable limit cycles exist, which limit the domain of attraction of stable equilibrium positions. The second question is about the size of oscillation amplitudes and frequencies in case of an unstable equilibrium position. If the amplitudes remain small, it may be acceptable in some practical applications that the stability border is exceeded from time to time; however, in order to avoid resonant interactions with other vibratory elements, the frequency has to be known, too. These questions have been firstly addressed in [13], where simple analytical conditions for global stability in the presence of linear viscous damping are derived. For the case of Coulomb friction and quadratic damping, a harmonic balance approach is applied to show that an unstable limit cycle exists around a stable equilibrium position, limiting the domain of attraction of stable solutions. Nevertheless, this investigations remain in some sense local, as they do not take the valve's switching property into account.
Since their impact vibrations promise the more interesting dynamical behavior, poppet valves have so far received most attention in the field of self-excited and forced vibration dynamical analysis. Numerical investigations of circuits involving poppet valves, disclosing non-smooth bifurcation scenarios and even chaotic solutions, can be found in [2,3,9,10,16]. The selfexcited and forced oscillations of systems involving spool valves are numerically analyzed in [12,26,27].
Despite the efforts already made, a systematic investigation (either numerically or analytically) of pressure control valve dynamics is still pending. One reason may be that even the simplest circuits include a large quantity of parameters, rendering numerical investigations extremely complex and time consuming. Therefore, in order to gain a better understanding of the dynamical behavior, it would be helpful to complement these numerical investigations by analytical studies. To conduct these, knowledge of characteristic quantities (length, time, pressure) is essential, in order to properly scale system parameters and variables [17]. The hybrid nature of every hydraulic circuit (meaning that two physical domains are included) renders the search for these characteristic quantities very counterintuitive. For example, the undamped eigenfrequency of the mechanical part (valve piston) is often taken to compute a characteristic time scale, but, as will be shown later, has generally little to do with the valve's natural frequency inside a hydraulic circuit.
Among the first attempts to investigate nonlinear valve dynamics analytically, in [23] the method of multiple scales is exploited to analyze the nonlinear response of a relief valve with cubic spring stiffness under harmonic excitation. This analysis is built upon in [20], where saddle node bifurcations of the steadystate solutions are observed. However, the model, on which the investigations are founded, does not include the aforementioned closed-loop property. Considering hydraulic servo-systems, singular perturbation techniques are exploited for the dynamical analysis and control design in [21,22,25]. In [18], an approach using singular perturbation and averaging technique is proposed to obtain a reduced-order model for an electrohydraulic gearbox actuator featuring a hydraulic pressure control valve. The frequency of a dithering signal introduces the slow time scale, and it is observed that the dynamics of the pressure control sensing chamber and pipeline resonances are much faster compared to valve piston and line pressure dynamics. The reducedorder model is constructed by omitting the fast variables. Subsequently, the slow motion is averaged in order to obtain an autonomous system. In [16], a minimal model of a pressure control poppet valve is analyzed. By simplification of the equilibrium position, an analytical approximation of the stability border, on which a Hopf bifurcation occurs, is computed. Furthermore, by means of a center manifold reduction, analytical approximations for the amplitude and frequency of the arising limit cycle oscillations are found, which, however, are very rough and which do not account for impacting solutions. The method of direct separation of motion is proposed in [14,15] to compute the effect of high-frequency excitation in form of a dithering signal on the stationary operating point of a pressure control valve with linear and quadratic damping.
Most of these studies have in common that an external frequency (corresponding to a dither signal or to a control reference signal) is introduced to compute a characteristic time scale. However, both self-excited oscillation and resonance frequencies generally correspond to natural system frequencies, and dimensionless models should include them in order to portray the system's magnitude orders in the right way. Within this discussion, a systematic approach both for computing this characteristic time scale and for computing characteristic measures of the dependent variables is proposed. The resulting dimensionless system equations clearly reveal two different magnitude orders, which motivates the utilization of an asymptotic method. The system is represented as a singularly perturbed third-order differential equation. Herein, the valve's main nonlinearity, which is its non-smooth flow characteristic, is not small. This nonlinearity must be accounted for in the zero-order solution of the asymptotic approximation. By means of singular perturbation and averaging techniques, a first-order approximation of the slow mani-

Pressure control valve model
This section presents a simple hydraulic circuit, which is used for the investigation of pressure control valve dynamics. First, a physical representation is derived, followed by a dimensional analysis of dependent and independent variables. Based upon this analysis, a dimensionless model of the circuit is presented.

Physical system representation
The hydraulic circuit investigated in this paper is depicted in Fig. 1. The principal components of the system are a pump, an ideal consumer and the valve itself. The consumer is modeled as a simple hydraulic resistance. Following the usual modeling approaches for the investigation of stability and dynamics of hydraulic systems [12,16,27], the governing equations are given by Herein, parameters m, k and r denote mass, stiffness and radius of the valve piston and the control force f c represents either a spring preload, a solenoid force or a static pressure force applied by a pilot stage. Spool position x and pressure p represent the dependent variables. Parameters q P and C h are the constant pump volume flow and system's hydraulic capacitance and the function q h , representing the system's nonlinear flow characteristic including non-smooth control edge flow and nonlinear volume flow to consumer, reads with γ V the turbulent flow coefficient, u the valve overlap, A C the area of the orifice representing the ideal consumer and with σ (x) and ( p) defined as Expression dẋ can be seen as a surrogate for various frictional effects in the system. Secondary physical effects such as hydraulic inductance, Bernoulli forces, position-dependent capacitance and Coulomb friction are not taken into consideration. The working principle of the valve is demonstrated in Fig. 2, showing numerical simulations of the transient responses of spool position x(t) and system pressure p(t) to a change in consumer demand (change of orifice area A C ) at t = 0.2 s and in power supply (change of pump volume flow q P ) at t = 0.4 s. As expected, stationary system pressure p 0 is maintained, whereas the stationary spool position x 0 varies depending on the difference between supply and load volume flow. Through enlargement of spool mass m at t = 0.6 s, the valve's equilibrium position is destabilized and limit cycle oscillations emerge. Here, Piston oscillation amplitudes are approximately 0.15 mm and pressure amplitudes are about 5.5 bar, which are rather high values. This is due to the fact that the simulation involves a high pump volume flow of q P = 20 l/m. The simulations are carried out with the parameter ranges given in Table 1, whose values are orientated on real valves. All numerical simulations are computed with Matlab, using its ode23t-solver including eventdetection of the switching edge.

Scaling of dependent and independent variables
In order to find the "right" scaling for a non-dimensional system representation, it is noted that the parameter values, with which the numerical simulations of the limit cycle oscillations in Fig. 2 at t > t 3 are performed, are chosen so that the valve is operated in an unstable parameter domain, but close to the stability border. It strikes that the piston oscillations x(t) crosses the switching edge, but only to a small extent. As a matter of fact, further numerical simulations show that this is always the case in proximity to the stability border. Hence, under the assumption that the piston oscillates approximately symmetrical around its equilibrium position, the distance x 0 − u between equilibrium position x = x 0 of valve spool and switching edge x = u seems to be a characteristic distance measure, at least in proximity to the bifurcation point, which is also the parameter domain of the most practical interest.
For the computation of this measure, it is shown in the sequel how to derive approximations for the stationary valve position x 0 and, as a by-product, set-pressure p 0 . Along with p T = 0Pa, the equilibrium conditions of Eqs. (1) and (2) read At this point, it is useful to assume that kx 0 f c . This is a typical case and is always strived for in pressure control valve design, because it corresponds to a precise preservation of the set-pressure level p 0 , as can be directly seen from Eq. (3). Following this assumption, the approximate stationary equilibrium takes a relatively simple form. The parameters in Table 1 correspond to a stationary set-pressure level p 0 ≈ 30bar. For the exact computation of the equilibrium position, the reader is referred to [16] or [27]. A characteristic time scale can be found by assuming that the piston inertia is small compared to pressure and damping forces, which is also a rather typical case [14]. Neglecting it in Eq. (1) allows to express the system pressure and its time derivative as Along with Eq. (2) and by performing a Taylor series expansion of the nonlinear volume flow function q h around the approximate equilibrium position, the representation of a linear damped oscillator is derived, whose approximate (considering, that k ∂q h ∂ x |x 0 ,p 0 ) undamped angular frequency is a natural choice for the introduction of a characteristic time scale. With standard parameter values given in Table 1, a characteristic angular frequency of ω ≈ 934rad/s is obtained. For the parameters of the limit cycle oscillations in Fig. 2, a value of ω ≈ 810rad/s is computed, which is an acceptable rough estimation for the actual oscillation angular frequency of ≈ 730rad/s.

Dimensionless system representation
The non-dimensional variables and parameters are introduced by using the characteristic quantities from Eqs. (5) and (6). This leads to the piecewise smooth (continuous, but not continuously differentiable) non-dimensional system representation Parameters M, K and F c are the dimensionless spool mass, spring stiffness and control force and Q P and with α V and α C the dimensionless spool diameter and orifice area. The values of the newly introduced dimensionless parameters are given in Table 2. With values in the third column, the valve's equilibrium position is unstable; it is therefore used in Sect. 3. Parameter values in the fourth column lead to a stable equilibrium position, and the simulation results in Sect. 4 are carried out with them.

Approximate solutions of self-excited oscillations
Within this section, the procedure for computing asymptotic approximations for the limit cycle amplitude and frequency of the pressure control valve model is described. The procedure combines classical singular perturbation with an averaging technique, as suggested, for example, in [28]. This procedure is closely related to the asymptotic analysis of partially damped systems proposed in [6] and generalized in [5]. Firstly, the system is reformulated as a third-order ordinary differential equation, which, as it appears, is singularly perturbed. Secondly, the singularly perturbed (fast) variable is approximated with an asymptotic expansion up to magnitude order ε. This approximation converts the singularly perturbed system into a regularly perturbed one, on which first-order averaging is applied.

Asymptotic approximation of the slow manifold
In order to clarify the magnitude orders, the new variable ξ = X − X 0 , representing deviations from the equilibrium position, is introduced and Eqs. (7) and (8) are reformulated as a single ordinary differential equation of third order in which P 0 = K X 0 + F c is the dimensionless stationary set-pressure. With the parameters given in Table 1, the dimensionless mass is small, i.e., M 1. Therefore, Eq. (10) is a singularly perturbed equation. The usual procedure in handling such equations is to separate slow variables ξ, ξ and the fast variable η = ξ , cf. [4,28]. Furthermore, to enable the application of asymptotic methods, Taylor series expansions of the two switching states of the nonlinear volume flow characteristic Q h , Eq. (9), are executed around the equilibrium position [ξ, ξ , ξ ] = [0, 0, 0] . By introducing a small parameter ε = M and neglecting terms of magnitude order ε 2 and higher, one obtains the piecewise smooth system representation in which all newly introduced parameters Table 3 Two sets of dimensionless parameter values for the system representation defined by Eqs. (11) and (12). The two sets based on the same physical parameter values as the sets in Table 2 Symbol are of magnitude order one, see Table 3.
In fact, parameter value N o is at least significantly higher than magnitude order ε and therefore considered as magnitude order 1. Parameters stem from the Taylor series expansions of the two switching states of Q h . The expressions εη and εK i ξ , i = {o, c}, associated with inertia and spring forces, are of magnitude order ε, which is in agreement with the assumptions stated in Sect. 2.2. The subscripts {o, c} are associated with the two switching states indicating whether the valve is open (ξ + X 0 > 0) or closed (ξ + X 0 ≤ 0). It is noted that the exact same system representation as in Eqs. (11) and (12) with the same magnitude orders can be derived for valves with pressure supply instead of volume flow supply, as studied in [12][13][14], but only if the magnitude orders of the difference between supply pressure and system pressure match the magnitude order of the system pressure itself. In terms of singular perturbation theory, an asymptotic approximation for the limit cycle oscillations corresponds to a regular or outer expansion of the solution [28]. For autonomous systems, this outer expansion can be associated with a so-called slow manifold M ε in phase space. By setting ε = 0 in the boundary layer Eq. (12), a continuous but non-smooth first-order approximation M 0 of the slow manifold M ε is obtained. This is a globally asymptotically stable (eigenvalue −1) solution of the equation in which ξ is a parameter and not a variable. Hence, all solutions of Eqs. (11) and (12) are attracted by the slow manifold. Figure 3 shows the three-dimensional phase space of the full system, in which the first-order approximation of the slow manifold M 0 is represented as a light gray, bent surface. One can clearly see that the trajectory, starting at the point indicated as IC, is strongly attracted by it and remains in an ε-vicinity for all time.

Asymptotic approximation of the system dynamics on the slow manifold
The flow on the slow manifold is described by Eq. (11). Herein, the fast variable η appears with magnitude order one, why an asymptotic expansion η ≈η = η 0 + εη 1 up to magnitude order ε has to be computed first. Substituting it into Eq. (12) and balancing terms with the same magnitude orders, the asymptotic approximation up to magnitude order ε can be directly computed to S. Schröders, A. Fidliñ by solving linear algebraic equations. It is noted that this approximation, in contrast to Eq. (13), is not continuous at ξ = −X 0 anymore. By substituting Eq. (14) into Eq. (11), the singularly perturbed system, Eqs. (11) and (12), is transformed into the regularly perturbed equation describing the flow on the slow manifold, on which averaging [24] is applied subsequently. Substituting the asymptotic expansions ξ ≈ξ = ξ 0 +εξ 1 , ξ ≈ξ = ξ 0 +εξ 1 into Eq. (15) and balancing terms with magnitude order ε 0 delivers the zero-order solutions of the slow variables with integration constants A, α, c 1 and c 2 . The constants c 1 and c 2 as well as the unknown switching times τ 0 , τ 1 and τ 2 can be expressed in dependence of A and α by solving transition conditions cf. Fig. 4. Note that the zero-order solution ξ 0 (τ ) is continuously differentiable also at τ = τ 2 , although this condition is not explicitly formulated. Solving the above conditions gives expressions Fig. 4 Structure of zero-order solutions ξ 0 (τ ) and ξ 0 (τ ) with according switching times, A = 1.5

Fig. 5 Switching times in dependence of amplitude
in dependence of A and α only, see Fig. 5.
A system in standard form for averaging can be obtained by varying constants A = A(τ ) and α = α(τ ) and using Eqs. (16) and (17) as a variable transformation with ξ = ξ 0 , ξ = ξ 0 . Inserting this transformation into Eq. (15), a non-smooth system of differential equations in the form for the variables A and α with right hand sides of magnitude order ε is derived. It is subsequently averaged with respect to dimensionless time τ providing a smooth right hand side. Variable A(τ ) represents the maxima of the limit cycle oscillations. The corresponding minima are computed along with Eq. (16). The dimensionless oscillation frequency can be computed from f = (τ 2 − τ 0 ) −1 , which does not depend on variable α. Therefore, the relevant equation for variable α is of no interest. The computation of the steady-state amplitude can be performed very efficiently numerically, either by calculating the root of the right hand side or by numerical time integration of Eq. (22). So far, extensive parameter studies have not yielded coexisting or unstable solutions corresponding to periodic solutions of the full system. This is an important result of practical interest, meaning that attraction domains of stable equilibrium solutions of the original system are not limited by unstable limit cycles. The existence condition for the solution of Eq. (22) results from the argument of C = arccos(X 0 /A), namely A > X 0 . That means that Eq. (22) is only valid for oscillations crossing the switching edge.
In case of oscillations not crossing the switching edge, A < X 0 , the averaging takes the form As a by-product of this result, the stability condition for the equilibrium position of the valve, corresponding to a stable trivial solution A = 0 of Eq. (23), can directly be computed, which with the dimensionless parameters from Sect. 3.1 and upon simplifying assumptions P 0 ≈ F c , X 0 ≈ 1 and K 1 < F leads to An increase of spring stiffness K or a decrease of the slope of the flow characteristic, corresponding to dimensionless parameter α V , stabilizes the equilibrium position, which is in accordance with the results presented in [16] and [13]. Because the right hand side of Eq. (20) is not continuous, the classical proofs for the validity of the asymptotic approximations do not hold here; however, they are validated by comparison with full time numerical simulations, see Fig. 6. Both the transient decay and the steady-state amplitude are reproduced excellently. Figure 7 shows a comparison between dimensionless oscillation frequencies f computed from original and averaged equations under variation of the small parameter ε = M. The results show a reasonable concordance, even for relatively high values of the small parameter ε. The same holds true for the oscillation minima and maxima, see Fig. 8. With these results, the question of what the limitations of the proposed method are, can be addressed. These are related to the small parameter ε. As an example, it is assumed that for ε < 0.2 the approximation is within an allowable error margin. In regard to the allowable physical parameters, e.g., system damping-starting from the standard values in Table 1-can be reduced to d ≈ 100 N s/m, or the hydraulic capacitance can be reduced to C h ≈ 10 −13 m 3 /Pa without the method failing.

Approximate solutions of forced oscillations
Apart from self-excited vibrations, another important scenario of practical interest in connection with nonlinear valve dynamics is the case of a valve being operated under stable conditions, but subject to periodic excitation. In real-world applications, these excitation can be evoked from pump pulsations, unstable neighboring valves or a dithering signal. It is expected that forced oscillation amplitudes remain not small in case of a resonant excitation, meaning that the excitation frequency matches the valve's natural frequency ω. Fig. 8 Comparison between non-dimensional stationary oscillation minima A min and maxima A calculated from full system and averaged equation The above described procedure can easily be modified in order to obtain asymptotic approximations of the system response to a harmonic excitation. The case of pump pulsations is considered here. For this case, the non-dimensional pump volume flow is assumed to be of the form in which parameter A Q is assumed to be of magnitude order one. Hence, the excitation amplitude is of magnitude order ε. Following the procedure described in Sect. 3, the governing equation for the slow variables, Eq. (15), takes the form It is noteworthy that in the present case it would be possible not to approximate the stationary solution of the fast dynamics by a series approach, but to compute it directly. This can be done by introducing the so-called strongly damped variable y = η − η 0 (compare [5]). Through this transformation, all expressions containing slow variables ξ , ξ are of magnitude order ε, which is why they can be seen as constant in this case. Since the equation of the fast dynamics is linear in η (and y, respectively), it can be explicitly solved.
In this exact solution, the harmonic excitation would cause an additional phase shift, which is, however, of magnitude order ε 2 . As expected, numerical simulations have shown for the present case that this does not result in a significant improvement of the approximate solution.
It is clear that the zero-order solutions of the slow variables, Eqs. (16) and (17), must be modified to ensure that the period of the zero-order solution matches the excitation period. There are several modification possibilities, and one has to keep in mind that the analytical manipulations necessary for the later averaging procedure must remain possible. From here on, it is assumed that the excitation frequency is nearby the resonance frequency; that means | Q − ω o | 1. To keep the equations as simple as possible, the zero-order solutions are proposed with phase variable ψ. Parameter κ has to be computed so that the period of the zero-order solution matches the period of the excitation. Solving the same conditions as in (18) and (19) but with τ 2 − τ 0 = 2π/ Q , expressions are derived. The further procedure is identical to that of the autonomous case, except that the cosine term in Eq. (24) for ξ + X 0 < 0 has to be approximated by a quadratic polynomial, which is computed by a Taylor series expansion around τ = (τ 0 + 2π/ Q + τ 1 )/2. This approximation is done in order to enable the analytical computation of the averaged equations, which are, however, too lengthy to write them here. The numerical simulations in Fig. 9 are executed with a stable parameter set with increased damping compared to that used in the previous section, see According to Fig. 10, also the stationary maxima and minima of the forced vibrations show a reasonable concordance. Stationary solution curves for A < X 0 even lie almost on top of each other, so that no difference is visible. Larger deviations can be seen for the case A > X 0 in proximity to the highest resonance peak. In particular, the first-order approximation underestimates the region, in which multiple solutions exist. The dashed lines correspond to unstable solutions, which most likely also exist in the full system, but are not computed here.

Conclusions
Within this contribution, the nonlinear dynamics of a hydraulic pressure control valve have been investigated by means of asymptotic approximations. After introducing the hydraulic circuit including pump, control valve and ideal consumer, the working principle has been demonstrated by numerical simulations of the transient system behavior, showing that the stationary set-pressure is maintained independently of the power supply and consumer demand. It has also been shown that the valve's internal feedback loop can cause stability problems and the occurrence of limit cycle oscillations. In order to prepare asymptotic analysis, considerations with regard to magnitude orders of system parameters and states have clarified that pressure and damping forces dominate inertia and spring forces. By neglecting the non-dominant forces, simple expressions for the equilibrium position and a characteristic system frequency have been derived, which are suitable for the scaling of dependent and independent variables of the system.
Because of the consideration of a very light piston, the dimensionless system representation is a singularly perturbed third-order differential equation. It has been shown that the limit cycle oscillations can be associated with a slow manifold in phase space, whose attraction properties have been demonstrated numerically and theoretically, the latter by showing that the corresponding first order approximation is a globally asymptotically stable solution of the unperturbed equation of the fast variable. The singularly perturbed problem has been transformed into a regularly perturbed equation by computing asymptotic approximations of the fast variable up to magnitude order ε.
The scaling considerations have also revealed that the switching nonlinearity of the valve's control edge flow is not small. Consequently, to capture the system's dominant nonlinearity, a non-smooth zero-order solution of the slow variable has been constructed, which subsequently has been used as a variable transformation to obtain a system in standard form for averaging. With this procedure, the non-smooth third-order system representation has been reduced to a first-order differential equation. The asymptotic approximations have been validated by comparison with numerical simulations of the original system equations and have shown to give a good estimate of limit cycle oscillation amplitude and frequency.
For the approximation of the forced system response due to pump pulsations, the zero-order solution of the slow variable has been modified in order to match its period and that of the excitation. The analytical expressions obtained from averaging are more complex; however, both transient and stationary solutions are well represented.
Funding Open Access funding enabled and organized by Projekt DEAL.

Compliance with ethical standards
Conflict of interest The authors declare that they have no conflict of interest.
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/.