# Nonlinear dynamics and onset of chaos in a physical model of a damper
pressure relief valve

Lukas Schickhofer  KTM Research & Development, Mattighofen, Austria
<EMAIL_ADDRESS>Chris G. Antonopoulos  Department of
Mathematical Sciences, University of Essex, Wivenhoe Park, United Kingdom

###### Abstract

Hydraulic valves, for the purpose of targeted pressure relief and damping, are
a ubiquitous part of modern suspension systems. This paper deals with the
analytical computation of fixed points of the dynamical system of a single-
stage shock absorber valve, which can be applied for the direct computation of
its system variables at equilibrium without brute-force numerical integration.
The obtained analytical expressions are given for the original dimensional
version of the model derived from continuity and motion equations for a
realistic valve setup. Furthermore, a large part of the study is devoted to a
systematic sensitivity analysis of the model towards crucial parameter
changes. Numerical investigation of a potential loss of stability and
following nonlinear oscillations is performed in multi-dimensional parameter
spaces, which reveals sustained valve vibrations at increased valve mass and
applied pretension force. The dynamical behaviour is analysed by phase space
orbits, as well as Fourier-transformed valve displacement data to identify
dominant frequencies. Chaotic indicators, such as Lyapunov exponents and the
Smaller Alignment Index (SALI), are utilized to understand the nature of the
oscillatory motion and to distinguish between order and chaos.

###### keywords:

Non-smooth dynamical system, Impact oscillator, Stability analysis, Parametric
study

###### MSC:

[2010] 00-01, 99-00

##  1 Introduction

It is long known that shock absorbers and automotive dampers, which contain
various forms of mechanical valves subjected to fluid pressure, show complex
and sometimes unexpected behaviour, such as instabilities, or sustained
vibrations and noise during operation. This is mainly due to the dominant
physical nonlinearities, which are also present in the underlying equations of
motion. These nonlinear terms come in several forms, such as coupled pressure-
displacement relationships, impact laws, or variable spring stiffness, which
collectively render the dynamical system non-smooth and therefore difficult to
analyse with usual methods of nonlinear dynamics [Kunze and Küpper, 2001, Di
Bernardo et al., 2008].

Although early works, especially by Lang [1977], Reybrouck [1994], and Talbott
and Starkey [2002], introduced potent models for the prediction of shock
absorber variables, these approaches were either algebraic and linearized in
nature and therefore not suited for studies of its nonlinear dynamics, or they
applied a large set of semi-empirical fitting parameters that did not allow
investigation of the actual physical effects and co-dependencies of real model
properties. Duym et al. [1997a] also gives a critical evaluation of early
damper models, stating some of their weaknesses and shortcomings, such as the
problems in capturing hysteresis and comparable nonlinear effects.
Consequently, models based on nonlinear differential equations and actual
material properties and physical parameters are required for the analysis of
realistic cause-effect relationships between model properties and behaviour,
as well as its dynamics. More recently, Farjoud et al. [2012] devised a set of
nonlinear differential equations for the modelling of hydraulic damper valves
and successfully predicted their opening characteristics at different shim
stack and orifice settings. Boggs et al. [2010] and Sikora [2018] applied
similar analytic methods to model the suspension damping behaviour in real-
world vehicle applications. Furthermore, Liao et al. [2021] derived at
suitable two-degree-of-freedom model for fluid-structure coupling within a
relief valve for the parameter optimization of underwater applications.

As a result of ongoing research and development in pressure relief valve
dynamics, particularly in the automotive industry, it became necessary to
understand instability mechanisms and the loss of dynamic stability of
hydraulic valves in operation. In an early numerical investigation of
instabilities in hydraulic valve systems, Hayashi et al. [1997] could identify
the Feigenbaum route to chaos by power spectra, bifurcation diagrams and
Lyapunov exponents for a simplified poppet control valve. Benaziz et al.
[2012] applied a first-order differential equation system to predict the
dynamics of spring valves and studied the global stability of hydraulic valves
under realistic parameter changes. Thus, they could map out the occurrence of
instabilities under various flow rates. Additionally, they provided a
nonlinear model for the analysis of valve vibration and the resulting
structure-borne noise [Benaziz et al., 2015]. Using various approaches like
the brute-force numerical integration of long time series of trajectories with
implicit schemes, as well as bifurcation analyses of non-impacting periodic
solutions, Hős and Champneys [2012] and Bazsó et al. [2014] were able to
identify initial instabilities and various bifurcation types, such as fold,
torus, and grazing bifurcations, in a non-dimensional model of a pressure
relief valve. Recently, an experimental study by Ma et al. [2019] confirmed
these instability mechanisms of pressure relief valves and specified the
occurrence of chatter instabilities by general stability maps for a wide range
of volumetric flow rates and geometric properties such as inlet pipe length.
The valve types in these studies are similar to the one studied here, but lack
the energy dissipation of a viscous fluid surrounding the valve. Furthermore,
to date no systematic investigation of the effect of parameter variations on
the global dynamics of suspension valves has been attempted. One reason for
this was the non-existence of a validated closed system of differential
equations suitable for the application of numerical integration and fixed-
point analysis.

Recently, such a physical model was introduced and successfully applied to
analyse inertia and dynamic stability of a typical automotive suspension valve
system [Schickhofer and Wimmer, 2022, Schickhofer, 2022]. In this work, its
derivative in the form of a coupled, dimensional nonlinear system of equations
for a shock absorber pressure relief valve is studied with respect to its
equilibrium states and its sensitivity to parameter changes. The chosen
approach is therefore a mixture of analytical methods for the derivation of
formulae for direct computation, as well as numerical methods for the
integration of the first-order system. Moreover, to understand whether the
system undergoes changes in the nature of its dynamics from order to chaos, or
vice versa, chaotic indicators such as Lyapunov exponents [Benettin et al.,
1980a, b] and the Smaller Alignment Index (SALI) [Skokos et al., 2003, 2004]
are applied. Although the studied system is generally dissipative, as
explained in detail for variations of the damping parameter in Sec. 3.3.5, it
could well display transient chaotic motion with subsequent loss of function.
It is therefore important, both from a fundamental, as well as an applied
perspective on pressure relief valves and shock absorber systems, to
understand and accurately predict such transitions and sensitivities.

Various model parameters are known to have a significant impact on the shock
absorber valve function and dynamic behaviour. However, little is known on
their precise impact on equilibrium states, transition to chaos, and
instabilities. This study attempts to answer three important questions:

  * 1.

Is it possible to find analytic expressions for the equilibrium states of the
main shock absorber variables, such as valve displacement, velocity, and
pressure?

  * 2.

How are the nonlinear valve dynamics affected by realistic parameter changes?

  * 3.

Can we distinguish between order and chaos in the motion of the investigated
dissipative dynamical system?

In the following, we will address these these questions through analytical and
numerical methods and we will return to them in Sec. 4 for a detailed
discussion.

##  2 Theory

The dynamical system to be studied is the three-dimensional model of a valve
for the purpose of pressure relief in a hydraulic shock absorber. A shock
absorber acts as part of the suspension system of various types of vehicles,
such as cars, motorcycles, or trains, to absorb wheel displacements as a
result of ground irregularities or impact. It contains a viscous fluid,
typically mineral oil, which is compressed by the motion of a rod in a piston.
Its valve is further displaced by the acting pressure. The pressure drop
across the valve initiates flow rates through a piston with orifices of
various cross sections, most notably the flow rate $Q_{b}$ through a constant
bleed orifice and the flow rate $Q_{v}$ through a variable valve opening. This
entire process leads to an efficient relief of pressure by various flow losses
through the piston and valve gap flow. For a detailed description of the shock
absorber concept, the reader is referred to Dixon [2008] and Schickhofer
[2022]. Figure 1 shows the forces acting on the single valve and the
volumetric flow rates across the considered control volume. The underlying
system of equations have been previously validated and successfully
implemented for shock absorber development [Schickhofer and Wimmer, 2022].
Below, the dimensional version of the equation system is analysed with respect
to its critical points.

\begin{overpic}[width=346.89731pt]{Figure1.pdf} \end{overpic} Figure 1:
Diagram of pressure force $F_{p}$, momentum force $F_{m}$, impact force
$F_{i}$, and pretension $F_{0}$ acting on the shock absorber valve (a), and
sketch of the control volume $V$ involving the total volumetric flow rate $Q$,
valve flow rate $Q_{v}$, and bleed flow rate $Q_{b}$ (b) used for the fully
integrated form of the continuity equation (2). A detailed description and
definition of forces and flow rates is given in Sec. 2.1. The boundaries at
$x=x_{1}$ and $x=x_{2}$ of the domain and valve motion are indicated.

###  2.1 Dimensional system

Essentially by applying Newton's second law on the valve shown in Fig. 1(a),
and the mass continuity equation on the control volume depicted in Fig. 1(b),
one arrives at the following equations:

$\displaystyle m\ddot{x}+c\dot{x}+kx$
$\displaystyle=F_{p}\left(p\right)+F_{m}\left(x,p\right)+F_{i}\left(x\right)-F_{0},$
(1) $\displaystyle\dot{p}$ $\displaystyle=\frac{1}{\beta
V}\left[Q-Q_{v}\left(x,p\right)-Q_{b}\left(p\right)\right].$ (2)

Equation (1)(a) models the valve as a driven oscillator of mass $m$, damping
$c$, and spring constant $k$, which is subjected to the forces given below.
The mineral oil used in the damper is assumed to be compressible with
effective compressibility $\beta$, and the integrated continuity equation
across the volume $V$ takes into account all flow rates introduced in Fig.
1(b).

##### Pressure force

The dominant force is given by the static pressure difference $\Delta
p=p-p_{0}$ across the valve between the variable upstream pressure $p$ and the
downstream pressure $p_{0}$ of the shock absorber piston. It depends crucially
on the imposed volumetric flow rate further described below and causes a
pressure force given by

$F_{p}\left(p\right)=A_{p}\Delta
p=A_{p}\left(p-p_{0}\right)\mathrel{\stackon[1.5pt]{=}{\stretchto{\scalerel*[width("=")]{\wedge}{\rule{4.30554pt}{12.91663pt}}}{0.5ex}}}A_{p}p.$
(3)

Here, $A_{p}$ is the total area on the valve on which the fluid pressure acts.
In many shock absorbers this area consists of the sum of several port
sections. If the reference pressure is considered zero ($p_{0}=0$), or if $p$
is assumed as gauge pressure, the expression in Eq. (3) simplifies to
$A_{p}p$.

##### Momentum force

Due to the redirection of the fluid jet by the valve surface, a momentum
change is taken into account by the force term

$F_{m}\left(x,p\right)=\frac{C_{f}\rho}{A_{p}}\left(\underbrace{\alpha\pi
d_{v}x}_{A_{v}\left(x\right)}C_{d,v}\sqrt{\frac{2p}{\rho}}\right)^{2},$ (4)

which is dependent on the valve opening $x$ and pressure $p$, but also on the
empirically determined momentum coefficient $C_{f}$, which controls the
relative amount of transferred momentum. Lang [1977] established this
parameter through experiments at $C_{f}\approx 0.3$, but it can vary
significantly between shock absorber geometries due to differences in piston
geometry design and the resulting flow paths, which directly affect the impact
angle of the jet on the valve surface and thus the effective area [Hős et al.,
2014]. The valve opening section $A_{v}\left(x\right)=\alpha Px=\alpha\pi
d_{v}x$ is a direct function of the product of perimeter and valve
displacement, which gives a generated surface for the flow. The flow
proportionality coefficient $\alpha$ is used to adjust the actual flow area
depending on the geometry. An alternative definition of the momentum force,
$F_{m}=\dot{m}u_{0}$, can be given as a function of mass flow rate $\dot{m}$
and inlet mean flow velocity $u_{0}$.

##### Impact force

The effect of impact at the lower and upper bounds $x_{1}$ and $x_{2}$ is
modelled by a force term

$F_{i}\left(x\right)=\begin{cases}-k_{i}\left(x-x_{1}\right),&\quad x\leq
x_{1},\\\ 0,&\quad x_{1}<x<x_{2},\\\ -k_{i}\left(x-x_{2}\right),&\quad x\geq
x_{2},\end{cases}$ (5)

which models the repulsion using a stiff spring constant
$k_{i}=$1\text{\times}{10}^{10}\text{\,}\mathrm{N}\text{\,}{\mathrm{m}}^{-1}$$
that is considerably larger than the system stiffness ($k_{i}\gg k$).

##### Pretension force

The constant force $F_{0}$ can be controlled via pretension of the spring
acting on the valve. Its purpose is the targeted delay of opening during
initial pressure increase.

##### Volumetric flow rates

The total volumetric flow rate $Q$ is introduced to the system by displacement
due to road irregularities. Furthermore, the definitions of the valve flow
rate

$Q_{v}\left(x,p\right)=C_{d,v}\underbrace{\alpha\pi
d_{v}x}_{A_{v}\left(x\right)}\sqrt{\frac{2\Delta p}{\rho}},$ (6)

and of the flow rate through the constant bleed orifice

$Q_{b}\left(p\right)=C_{d,b}A_{b}\sqrt{\frac{2\Delta p}{\rho}},$ (7)

follow directly from the Bernoulli equation, which is considered to accurately
describe the flow across the piston under the premise that the discharge
coefficients $C_{d,v}$ and $C_{d,b}$ are accurately chosen to reflect the
transition from laminar to turbulent flow, as suggested by Segel and Lang
[1981]. As introduced through Fig. 1, the bleed orifice is a constriction of
constant cross section in the centre of the shock absorber piston, while the
valve lift leads to a variable opening. The effect of both paths is a
cumulative flow rate and pressure drop across the shock absorber.

Using the above definitions, Eq. (1)-(2) can be rewritten as a three-
dimensional, first-order system

$\displaystyle\dot{y}_{1}$ $\displaystyle=y_{2},$ (8)
$\displaystyle\dot{y}_{2}$
$\displaystyle=\frac{1}{m}\left[-cy_{2}-ky_{1}+F_{p}\left(y_{3}\right)+F_{m}\left(y_{1},y_{3}\right)+F_{i}\left(y_{1}\right)-F_{0}\right],$
(9) $\displaystyle\dot{y}_{3}$ $\displaystyle=\frac{1}{\beta
V}\left[Q-Q_{v}\left(y_{1},y_{3}\right)-Q_{b}\left(y_{3}\right)\right],$ (10)

where $\vec{y}=\left(y_{1},y_{2},y_{3}\right)=\left(x,\dot{x},p\right)$.

The system (8)-(10) can subsequently be used for numerical integration within
a large parameter space. By varying the flow rate $Q$ and certain valve
properties, the dynamics and transitions from stable to unstable states are
investigated in Sec. 3.2 and Sec. 3.3. Numerical integration is carried out by
schemes capable of dealing with stiff equations and non-smooth systems, such
as ODE15s or ODE23s in MATLAB. For the data displayed in multi-dimensional
parameter spaces, such as the plots depicting valve displacement as a function
of flow rate and model parameters in Sec. 3.3, task farming on a computing
cluster is applied, which allows for the parallel integration of many single-
core cases for sufficiently long runtimes.

###  2.2 Chaotic indicators

####  2.2.1 Lyapunov exponents

Lyapunov characteristic exponents (LCE) are a well-known tool to gauge the
rate of separation from an orbit point in $\mathbb{R}^{m}$ along its $m$
orthogonal directions. In three dimensions, this leads to a set of three
Lyapunov exponents,

$L_{k}=\lim_{t\to\infty}{\left(\frac{1}{t}\ln{\left[\frac{\vec{w}_{k}\left(t\right)}{\vec{w}_{k}\left(0\right)}\right]}\right)}=\lim_{t\to\infty}{\left(\frac{1}{t}\ln{\left[\vec{r}_{k}\left(t\right)\right]}\right)},$
(11)

indicating the shrinking or stretching of orthogonal axes
$\vec{r}_{k}\left(t\right)$ that span the state space with $k=1,\ldots,m$ at
time $t$. If the largest Lyapunov exponent settles to a non-zero, positive
value, this gives an indication of chaotic motion, since infinitesimally close
initial conditions result in orbits that diverge exponentially fast in time
[Alligood et al., 1996].

####  2.2.2 Smaller Alignment Index (SALI)

An alternative to Lyapunov exponents for the detection of chaos is the Smaller
Alignment Index (SALI) [Skokos et al., 2003, 2004], which is defined as

$\operatorname{SALI}\left(t\right)=\min\left\\{\left\lVert\frac{\vec{w}_{1}\left(t\right)}{\left\lVert\vec{w}_{1}\left(t\right)\right\rVert}+\frac{\vec{w}_{2}\left(t\right)}{\left\lVert\vec{w}_{2}\left(t\right)\right\rVert}\right\rVert,\left\lVert\frac{\vec{w}_{1}\left(t\right)}{\left\lVert\vec{w}_{1}\left(t\right)\right\rVert}-\frac{\vec{w}_{2}\left(t\right)}{\left\lVert\vec{w}_{2}\left(t\right)\right\rVert}\right\rVert\right\\}.$
(12)

It uses the fact that any two randomly chosen initial deviation vectors
$\vec{w}_{1}\left(0\right)$, $\vec{w}_{2}\left(0\right)$ will eventually
become aligned with the most unstable direction with the angle between them
tending to zero. Since one is only interested in the direction of the
deviation vectors, they are only normalized at each timestep and not
orthogonalized, which is a crucial difference to the computation of Lyapunov
exponents and saves computing time. In the case of chaotic orbits the
normalized vectors align and point to the same (or exactly opposite)
direction. Thus the norm of their sum (_antiparallel alignment index_) or
difference (_parallel alignment index_) tends to zero. It follows from Eq.
(12) that generally
$\operatorname{SALI}\left(t\right)\in\left[0,\sqrt{2}\right]$, with
$\lim_{t\to\infty}{\operatorname{SALI}\left(t\right)}=0$ being a clear
indication of chaos. SALI provides information on the instantaneous changes of
the dynamics and is therefore an important additional quantity for chaotic
indication. While Lyapunov exponents show the long-term persistent trend of a
dynamical system either towards order or chaos, in particular after
convergence and considerable integration time, SALI is able to identify also
transient chaos and is especially relevant for the analysis of short-lived
oscillations of realistic, dissipative systems, such as the one studied in
this paper.

###  2.3 Validation

Extensive validation of the model (8)-(10) has been carried out both against
fluid-structure interaction (FSI) simulation data in Schickhofer and Wimmer
[2022] as well as against experimental test bench data of damping curves in
Schickhofer [2022]. A typical approach in suspension engineering has been
carried out in these measurements: Using a time-resolved dynamometer with
controlled sinusoidal mass flow input of ground excitation $w_{g}=2\pi f_{g}$,

$\dot{m}\left(t\right)=\dot{m}_{max}\left(t\right)\sin\left(\omega_{g}t\right),$

the pressure drop across the hydraulic shock absorber valve has been measured
and compared to the pressure output of the model. Here, a reference frequency
of ground excitation of $f_{g}=$8\text{\,}\mathrm{Hz}$$ has been chosen. Since
the pressure reduction across the valve system is directly coupled to the
other model variables, such as valve opening displacement, this constitutes a
sufficient setup for valve model validation. Figure 2 shows a comparison of
numerically obtained pressure loss from the system (8)-(10) against test bench
data of a realistic shock absorber including a hydraulic pressure relief valve
of $0$, $15$, and $29\text{\,}\operatorname{Clicks}$, where
$\text{\,}\operatorname{Clicks}$ refers to a common configuration unit for the
adjustment of the bleed orifice area (${A}_{b}=$0\text{\,}{\mathrm{mm}}^{2}$$
at $0\text{\,}\operatorname{Clicks}$,
${A}_{b}=$1.4\text{\,}{\mathrm{mm}}^{2}$$ at
$15\text{\,}\operatorname{Clicks}$, ${A}_{b}=$2.7\text{\,}{\mathrm{mm}}^{2}$$
at $29\text{\,}\operatorname{Clicks}$). The other baseline model parameters of
the check valve are given by the values from A. It can be seen that the
agreement between the model's prediction and the test bench measurements is
satisfactory and lies consistently below $4\text{\,}\mathrm{\char 37\relax}$
for the compression stroke of positive pressure drop values and the rebound
stroke of negative values in Fig. 2.

\begin{overpic}[width=216.81pt]{Figure2.pdf} \end{overpic} Figure 2:
Numerically computed pressure loss of a realistic damper including a pressure-
relief valve modelled by the system (8)-(10) compared against test bench
measurements.

##  3 Results

In the following, both direct analytical approaches, as well as numerical
integration of the first-order system (8)-(10) are performed to investigate
the characteristics of the valve dynamics. While Sec. 3.1 offers a
mathematical formulation for the immediate computation of equilibrium states
in the valve motion, Sec. 3.2-3.2 explore the physical consequences of crucial
parameter changes of the system. The three-dimensional dynamical system, as
introduced by Eq. (8)-(10) is solved with the values for its model parameters
as presented in A. Wherever specific parameters are varied, this is noted in
text and figures.

###  3.1 Fixed points

The investigation of fixed points of a dynamical system is valuable for a
number of reasons: First of all, they give the long-term state to which the
system converges after all initial transients have disappeared and constitute
important stable solutions for a given set of initial and boundary conditions.
Additionally, if analytical expressions can be found for these critical
points, no numerical integration is necessary for their computation, which is
an important mathematical property.

In order to find the fixed points of the system (8)-(10), its gradients must
vanish such that

$\dot{\vec{y}}_{e}=\left.\frac{\partial\vec{y}}{\partial
t}\right|_{\vec{y}_{e}}=\vec{f}\left(\vec{y}_{e}\right)=\vec{0},$ (13)

which gives

$\displaystyle y_{2}$ $\displaystyle=0,$ (14) $\displaystyle-
ky_{1}+F_{p}\left(y_{3}\right)+F_{m}\left(y_{1},y_{3}\right)+F_{i}\left(y_{1}\right)-F_{0}$
$\displaystyle=0,$ (15) $\displaystyle
Q-Q_{v}\left(y_{1},y_{3}\right)-Q_{b}\left(y_{3}\right)$ $\displaystyle=0.$
(16)

Equation (14) simply enforces zero velocity of the valve at equilibrium, while
Eq. (16) recovers essentially the conservation of volume across the control
domain. Thus, there will be no compressibility effects at fixed points of the
system. As a result, Eq. (15) is analysed as the main equation of state
involving valve opening $y_{1}$ and acting pressure $y_{3}$:

$\displaystyle y_{1}$
$\displaystyle=\frac{1}{k}\left(F_{p}\left(y_{3}\right)+F_{m}\left(y_{1},y_{3}\right)+F_{i}\left(y_{1}\right)-F_{0}\right)$
$\displaystyle=\frac{1}{k}\left[A_{p}y_{3}+\frac{C_{f}\rho}{A_{p}}\left(\alpha\pi
d_{v}C_{d,v}y_{1}\sqrt{\frac{2y_{3}}{\rho}}\right)^{2}-\underbrace{k_{i}\left(y_{1}-x_{1/2}\right)}_{\forall
y_{1}\leq x_{1}\land y_{1}\geq x_{2}}-F_{0}\right].$ (17)

This is a quadratic equation that can be rewritten to give an expression for
$y_{1}$. By rearranging and defining constant coefficients, we obtain

$\displaystyle\underbrace{\frac{2C_{f}}{kA_{p}}\left(\alpha\pi
d_{v}C_{d,v}\right)^{2}}_{a}y_{1}^{2}y_{3}-\left(1+\frac{k_{i}}{k}\right)y_{1}+\underbrace{\frac{A_{p}}{k}}_{b}y_{3}+\underbrace{\frac{k_{i}}{k}}_{c}\underbrace{x_{1/2}}_{d}-\underbrace{\frac{F_{0}}{k}}_{e}$
$\displaystyle=0,$

which can be rewritten as

$ay_{3}y_{1}^{2}-\left(1+c\right)y_{1}+by_{3}+cd-e=0.$ (18)

Equation (16) can be further applied to get a relationship between
displacement $y_{1}$ and pressure $y_{3}$:

$\displaystyle y_{3}$
$\displaystyle=\frac{\rho}{2}\left(\frac{Q}{C_{d,v}\alpha\pi
d_{v}y_{1}+C_{d,b}A_{b}}\right)^{2}$ $\displaystyle=\underbrace{\frac{\rho
Q^{2}}{2}}_{h}\frac{1}{\left(\underbrace{C_{d,v}\alpha\pi
d_{v}}_{f}y_{1}+\underbrace{C_{d,b}A_{b}}_{g}\right)^{2}},$

which results in

$y_{3}=\frac{h}{\left(fy_{1}+g\right)^{2}}.$ (19)

By using Eq. (19) for Eq. (18), one ends up with an equation of third order in
$y_{1}$:

$\displaystyle
ay_{1}^{2}\frac{h}{\left(fy_{1}+g\right)^{2}}-\left(1+c\right)y_{1}+b\frac{h}{\left(fy_{1}+g\right)^{2}}+cd-e$
$\displaystyle=0,$ $\displaystyle
ahy_{1}^{2}-\left(1+c\right)\left(fy_{1}+g\right)^{2}y_{1}+\left(cd-e\right)\left(fy_{1}+g\right)^{2}+bh$
$\displaystyle=0.$ (20)

In order to alleviate the analytical treatment, three cases are distinguished
that reflect typical realistic conditions of the shock absorber:

##### Case 1: $F_{i}=0$, $F_{0}=0$

For motion of the valve with vanishing impact force and pretension, this case
gives a reasonable approximation. It leads to the following expression:

$\displaystyle ahy_{1}^{2}-\left(fy_{1}+g\right)^{2}y_{1}+bh=0,$
$\displaystyle f^{2}y_{1}^{3}+\left(2fg-ah\right)y_{1}^{2}+g^{2}y_{1}-bh=0,$
$\displaystyle y_{1}^{3}+\underbrace{\frac{2fg-
ah}{f^{2}}}_{A}y_{1}^{2}+\underbrace{\frac{g^{2}}{f^{2}}}_{B}y_{1}-\underbrace{\frac{bh}{f^{2}}}_{C}=0.$
(21)

Thus, the resulting relationship can be written as a third-order equation in
$y_{1}$ of the type $y_{1}^{3}+Ay_{1}^{2}+By_{1}-C=0$, which can be solved by

$\displaystyle y_{1}\approx$ $\displaystyle
0.26457\sqrt[3]{-2A^{3}+5.1962\sqrt{-4A^{3}C-A^{2}B^{2}+18ABC+4B^{3}+27C^{2}}+9AB+27C}$
$\displaystyle-\frac{0.41997\left(3B-A^{2}\right)}{\sqrt[3]{-2A^{3}+5.1962\sqrt{-4A^{3}C-A^{2}B^{2}+18ABC+4B^{3}+27C^{2}}+9AB+27C}}$
$\displaystyle-0.33333A.$ (22)

##### Case 2: $F_{i}=0$, $F_{0}\neq 0$

If pretension force is applied and regular smooth motion without impact is
considered, we get

$\displaystyle
ahy_{1}^{2}-\left(fy_{1}+g\right)^{2}y_{1}-e\left(fy_{1}+g\right)^{2}+bh=0,$
$\displaystyle f^{2}y_{1}^{3}+\left(2fg-
ah+ef\right)y_{1}^{2}+\left(g^{2}+2efg\right)y_{1}-bh+eg^{2}=0,$
$\displaystyle y_{1}^{3}+\underbrace{\frac{2fg-
ah+ef}{f^{2}}}_{A}y_{1}^{2}+\underbrace{\frac{g^{2}+2efg}{f^{2}}}_{B}y_{1}-\underbrace{\frac{bh-
eg^{2}}{f^{2}}}_{C}=0,$ (23)

which can be solved in a similar fashion as above:

$\displaystyle y_{1}\approx$ $\displaystyle
0.26457\sqrt[3]{-2A^{3}+5.1962\sqrt{-4A^{3}C-A^{2}B^{2}+18ABC+4B^{3}+27C^{2}}+9AB+27C}$
$\displaystyle-\frac{0.41997\left(3B-A^{2}\right)}{\sqrt[3]{-2A^{3}+5.1962\sqrt{-4A^{3}C-A^{2}B^{2}+18ABC+4B^{3}+27C^{2}}+9AB+27C}}$
$\displaystyle-0.33333A.$ (24)

##### Case 3: $F_{i}\neq 0$, $F_{0}\neq 0$

For the case of acting impact and pretension forces, which involves all terms
of Eq. (20), we can attempt a comparison of the order of present forces. Since
the stiffness of the impact force term is several orders of magnitude larger
than the system's stiffness ($k_{i}/k\approx 10^{6}$), we can approximate Eq.
(20) by

$-c\left(y_{1}+d\right)=0,$ (25)

which leads to the fixed points

$y_{1}=d=x_{1/2},$ (26)

at the lower and upper bound of the impact oscillator (cf. Fig. 1). These are
pseudo-equilibria of the non-smooth dynamical system, which are defined
exclusively by the discontinuities of the domain.

With the analytically obtained expressions for equilibrium displacements of
Eq. (22) and Eq. (26), together with the definition of Eq. (19) for the
related pressure, one has the possibility of computing those valve parameters
without numerically integrating system (8)-(10). Figure 3 shows a comparison
of analytical and numerical solutions for relevant flow rates. After initial
transients have subsided, the time-dependent numerical result settles to an
equilibrium and the two solutions are identical.

\begin{overpic}[width=433.62pt]{Figure3.pdf} \end{overpic} Figure 3:
Analytically computed equilibrium solutions for displacement and pressure
compared to the results from numerical integration of the system (8)-(10) for
typical volumetric flow rates of
$[2,4,8,14,20]\,$\text{\,}\mathrm{L}\mathrm{/}\min$$.

###  3.2 Valve opening characteristics

An important aspect of valve dynamics is the opening behaviour during increase
of the flow rate. In this section, we investigate the possibility of popping-
off instabilities or valve chattering during opening, as well as the
occurrence of oscillations during impact at the boundaries of the non-smooth
system (8)-(10). The impact force (5) is acting at
$x_{1}=$0\text{\,}\mathrm{mm}$$ and $x_{2}=$1.7\text{\,}\mathrm{mm}$$ (cf. A)
in this setup, which are valve seat values for shock absorber applications.

Depending on the flow velocity magnitudes, the valve might show delays in
settling to an equilibrium state. This effect is demonstrated in Fig. 4 and
occurs due to initial vibrations and further violent chattering during
closing, especially at intermediate to high flow rates. Initially, there are
persistent instabilities at low volume flow due to the fact that the forces
acting on the valve are of similar order of magnitude (cf. Eq. (1)). After
recovery to a constant delay of equilibrium above approximately
$50\text{\,}\mathrm{L}\mathrm{/}\min$, it increases linearly with the applied
volume flow at higher flow rates at above
$200\text{\,}\mathrm{L}\mathrm{/}\min$. Thus, there is a clear indication of
an ideal range of operation of the shock absorber valve with respect to the
fluid flow it can handle. Predictable and largely constant opening
characteristics are obtained for volumetric flow rates of
$60$–$190\text{\,}\mathrm{L}\mathrm{/}\min$.

\begin{overpic}[width=433.62pt]{Figure4.pdf} \end{overpic} Figure 4: Dynamics
of valve opening for a range of $0$–$1000\text{\,}\mathrm{L}\mathrm{/}\min$
(a) and $0$–$100\text{\,}\mathrm{L}\mathrm{/}\min$ (b). The time needed for
reaching the opening amplitude at $1.7\text{\,}\mathrm{mm}$ is shown together
with a black dashed line indicating the point at which the solution has
settled to equilibrium and all gradients of the first-order system (8)-(10)
have vanished.

###  3.3 Parametric studies

The impact of various parameter changes is studied below, where crucial
quantities such as the bleed orifice area, mass, or pretension are varied (cf.
Sec. 2.1). All parameter ranges are chosen based on a realistic check valve
model as it would occur in a hydraulic shock absorber and using the baseline
properties from A.

####  3.3.1 Bleed orifice area

The cross-sectional area of the bleed orifice regulates the flow through the
constant orifice of the valve, as detailed in Eq. (7). It thereby generates a
pressure loss due to flow constriction and subsequent expansion, which occurs
in parallel to the flow losses through the spring-loaded valve. Fig. 5
reflects the fact that changes in the bleed orifice area have a linear effect
on the valve dynamics. The orbits change only slightly in magnitude, but the
underlying motion pattern stays the same. Therefore, all dynamic motion
captured in Fig. 5 can be considered stable.

\begin{overpic}[width=433.62pt]{Figure5.pdf} \end{overpic} Figure 5: Orbits in
phase space of displacement, velocity, and pressure for specific values of the
bleed orifice area and increasing flow rates of
$Q=$10\text{\,}\mathrm{L}\text{\,}{\mathrm{min}}^{-1}$$ (a),
$Q=$50\text{\,}\mathrm{L}\text{\,}{\mathrm{min}}^{-1}$$ (b),
$Q=$100\text{\,}\mathrm{L}\text{\,}{\mathrm{min}}^{-1}$$ (c), and
$Q=$500\text{\,}\mathrm{L}\text{\,}{\mathrm{min}}^{-1}$$ (d).

####  3.3.2 Mass

A crucial parameter for the dynamic behaviour of a pressure relief valve is
its mass. It can be seen from Fig. 6 that an increase of mass has a
potentially destabilizing effect due to larger inertia. This leads to violent
_chattering_ , a boundary effect of non-smooth dynamical systems and impact
oscillators, which entails a series of impacts that cause vibrations of the
valve in real applications [Mora, 2013]. Table 1 indicates the fundamental
stability properties of the valve at different mass and flow rates. This
nonlinear phenomenon intensifies at higher flow rates and causes sustained
loops, or even surfaces, along which the trajectories oscillate (cf. Fig.
6(b)-(d)). The tendency of the valve to show stable oscillations thereby
clearly increases at higher mass, as further demonstrated by the Fourier
transformation of the displacement time series data shown in Fig. 7.
Additionally, the sections through the parameter space in Fig. 8 reveal ranges
of volumetric flow rates of sustained oscillations visualized in black that
extend to lower flow rates within the sections of higher masses. Furthermore,
initial instabilities during valve opening are bleeding into lower flow rates
at intermediate mass, as indicated by the isolated patches of higher
displacement in Fig. 8(b).

\begin{overpic}[width=433.62pt]{Figure6.pdf} \end{overpic} Figure 6: Phase space with trajectories for various valve masses and increasing flow rates of $Q=$10\text{\,}\mathrm{L}\text{\,}{\mathrm{min}}^{-1}$$ (a), $Q=$50\text{\,}\mathrm{L}\text{\,}{\mathrm{min}}^{-1}$$ (b), $Q=$100\text{\,}\mathrm{L}\text{\,}{\mathrm{min}}^{-1}$$ (c), and $Q=$500\text{\,}\mathrm{L}\text{\,}{\mathrm{min}}^{-1}$$ (d). Table 1: Stability map summarizing the results of Fig. 6 for various mass values and flow rates, where the symbol $\times$ means _stable_ and $\circ$ means _unstable_. Flow rate $Q$ Mass $m$ | $0.0005\text{\,}\mathrm{kg}$ | $0.0105\text{\,}\mathrm{kg}$ | $0.0205\text{\,}\mathrm{kg}$ | $0.0305\text{\,}\mathrm{kg}$
---|---|---|---|---
$10\text{\,}\mathrm{L}\text{\,}{\mathrm{min}}^{-1}$ | $\times$ | $\times$ | $\times$ | $\times$
$50\text{\,}\mathrm{L}\text{\,}{\mathrm{min}}^{-1}$ | $\times$ | $\circ$ | $\circ$ | $\circ$
$100\text{\,}\mathrm{L}\text{\,}{\mathrm{min}}^{-1}$ | $\times$ | $\circ$ | $\circ$ | $\circ$
$500\text{\,}\mathrm{L}\text{\,}{\mathrm{min}}^{-1}$ | $\times$ | $\circ$ | $\circ$ | $\circ$
\begin{overpic}[width=433.62pt]{Figure7.pdf} \end{overpic} Figure 7:
Displacement time series (a) and Fourier transformation (b) for realistic
volumetric flow rates and the chosen mass values.
\begin{overpic}[width=433.62pt]{Figure8.pdf} \end{overpic} Figure 8: Parameter
space of valve displacement as a function of mass, time, and flow rate in a
range of $0$–$100\text{\,}\mathrm{L}\mathrm{/}\min$ (a) and
$0$–$50\text{\,}\mathrm{L}\mathrm{/}\min$ (b).

####  3.3.3 Pretension

The pretension force applied to a shock absorber valve delays opening and
causes an overall higher pressure drop. However, there is a problematic side
effect of using a loaded spring for pretension: Particularly at lower
volumetric flow rate and during the initial phase of opening there is a
tendency to trigger instability and oscillations, which is visible in Fig. 9
and Tab. 2. In Fig. 9(a) the orbits in phase space exist on geometric
structures similar in appearance to a chaotic attractor. Also Fig. 10 gives
indications to aperiodic motion: While strong oscillatory behaviour can be
recognized for all considered pretension force values, the motion undergoes a
transition from relatively unordered vibrations without distinct frequencies
at lower flow rates to regular periodic oscillations at frequencies of
$2100\text{\,}$–$2900\text{\,}\mathrm{Hz}$. The transient character of this
unwanted valve response becomes obvious by investigating the parameter spaces
of Fig. 11, where the black areas again indicate strong vibrations that are
dominant when the terms for opposing forces in Eq. (1), such as the pressure
force and momentum force on the one hand, and the pretension force on the
other hand, become similar in magnitude. As soon as the volumetric flow rate
is larger than a threshold value for the considered pretension, valve
oscillations subside. This effect is also underlined by Tab. 2, which shows
stable motion at intermediate to high flow rates above
$100\text{\,}\mathrm{L}\text{\,}{\mathrm{min}}^{-1}$ for all pretension
values.

\begin{overpic}[width=433.62pt]{Figure9.pdf} \end{overpic} Figure 9: Orbits in phase space of displacement, velocity, and pressure for specific values of pretension and increasing flow rates of $Q=$10\text{\,}\mathrm{L}\text{\,}{\mathrm{min}}^{-1}$$ (a), $Q=$50\text{\,}\mathrm{L}\text{\,}{\mathrm{min}}^{-1}$$ (b), $Q=$100\text{\,}\mathrm{L}\text{\,}{\mathrm{min}}^{-1}$$ (c), and $Q=$500\text{\,}\mathrm{L}\text{\,}{\mathrm{min}}^{-1}$$ (d). Table 2: Stability map for the considered pretension values and flow rates of Fig. 9, where $\times$ indicates a _stable_ setup and $\circ$ an _unstable_ one. Flow rate $Q$ Pretension $F_{0}$ | $0\text{\,}\mathrm{N}$ | $30\text{\,}\mathrm{N}$ | $60\text{\,}\mathrm{N}$ | $90\text{\,}\mathrm{N}$
---|---|---|---|---
$10\text{\,}\mathrm{L}\text{\,}{\mathrm{min}}^{-1}$ | $\times$ | $\circ$ | $\circ$ | $\circ$
$50\text{\,}\mathrm{L}\text{\,}{\mathrm{min}}^{-1}$ | $\times$ | $\circ$ | $\circ$ | $\circ$
$100\text{\,}\mathrm{L}\text{\,}{\mathrm{min}}^{-1}$ | $\times$ | $\times$ | $\circ$ | $\circ$
$500\text{\,}\mathrm{L}\text{\,}{\mathrm{min}}^{-1}$ | $\times$ | $\times$ | $\times$ | $\times$
\begin{overpic}[width=433.62pt]{Figure10.pdf} \end{overpic} Figure 10:
Displacement in the time (a) and frequency domain (b) for a range of
volumetric flow rates. \begin{overpic}[width=433.62pt]{Figure11.pdf}
\end{overpic} Figure 11: Parameter space with displacement as a function of
pretension, time, and volumetric flow rate in the range of
$0$–$100\text{\,}\mathrm{L}\mathrm{/}\min$ (a) and
$0$–$50\text{\,}\mathrm{L}\mathrm{/}\min$ (b).

In order to further investigate the nature of the oscillatory motion described
above, chaotic indicators are applied to the cases with pretension leading to
orbits resembling chaotic trajectories. As visible in Fig. 9, these are in
particular pretensions of $F_{0}=$30\text{\,}\mathrm{N}$$ and
$F_{0}=$60\text{\,}\mathrm{N}$$. Figure 12 shows the evolution of the Lyapunov
exponents in time. It becomes clear from the double-logarithmic plot that it
tends to zero. Nevertheless, intermediate ranges around $t=10^{-4}$ show a
highly unsteady behaviour of the Lyapunov exponents, which are possibly
related to transient chaos. Moreover, the time series of SALI in Fig. 13 shows
values that never fully settle below $\operatorname{SALI}\left(t\right)\leq
10^{-12}$. This confirms the notion of the consistently decreasing Lyapunov
exponents from above and gives convincing reason to believe that the observed
orbits of Fig. 9 are indeed ordered and periodic, despite their Fourier-
transformed trajectories of Fig. 10(b) showing only weak frequency
information. However, the initial SALI time evolution in Fig. 13(a)-(b) for
the pretension values of $F_{0}=$30\text{\,}\mathrm{N}$$ and
$F_{0}=$60\text{\,}\mathrm{N}$$ has a sharp exponential decline hinting once
again to transient chaotic behaviour.

\begin{overpic}[width=433.62pt]{Figure12.pdf} \end{overpic} Figure 12:
Lyapunov exponents of the system at a pretension of
$F_{0}=$30\text{\,}\mathrm{N}$$ (a) and $F_{0}=$60\text{\,}\mathrm{N}$$ (b)
for $t\to$1\text{\,}\mathrm{s}$$ in double-logarithmic plot.
\begin{overpic}[width=433.62pt]{Figure13.pdf} \end{overpic} Figure 13: Time
evolution of SALI at a pretension of $F_{0}=$30\text{\,}\mathrm{N}$$ (a) and
$F_{0}=$60\text{\,}\mathrm{N}$$ (b) for $t\to$1\text{\,}\mathrm{s}$$ in
logarithmic plot.

####  3.3.4 Spring stiffness

An increase in stiffness of the spring loading the valve shows a strong
influence on the equilibrium solution to which the motion settles. However,
the basic dynamics of the blowing-off or popping-off process is not affected
(cf. Fig. 14). In particular, there is a clear initial overshooting of
trajectories at intermediate to high stiffnesses of around
$5100\text{\,}\mathrm{N}\text{\,}{\mathrm{m}}^{-1}$ to
$7600\text{\,}\mathrm{N}\text{\,}{\mathrm{m}}^{-1}$, as visible in Fig. 15.
This is probably due to a delayed reaction of the valve to the stronger spring
forcing. Furthermore, Fig. 16 shows a gradual diminution of the opening
amplitude at increasing values of spring stiffness, especially at small to
intermediate flow rates up to approximately
$70\text{\,}\mathrm{L}\text{\,}{\mathrm{min}}^{-1}$.

\begin{overpic}[width=433.62pt]{Figure14.pdf} \end{overpic} Figure 14: Phase
space with trajectories for the system's variables at various values of the
spring stiffness and at flow rates of
$Q=$10\text{\,}\mathrm{L}\text{\,}{\mathrm{min}}^{-1}$$ (a),
$Q=$50\text{\,}\mathrm{L}\text{\,}{\mathrm{min}}^{-1}$$ (b),
$Q=$100\text{\,}\mathrm{L}\text{\,}{\mathrm{min}}^{-1}$$ (c), and
$Q=$500\text{\,}\mathrm{L}\text{\,}{\mathrm{min}}^{-1}$$ (d).
\begin{overpic}[width=433.62pt]{Figure15.pdf} \end{overpic} Figure 15: Time
series of displacement data for various stiffness values and flow rate ranges
of $0$–$100\text{\,}\mathrm{L}\mathrm{/}\min$ (a) and
$0$–$50\text{\,}\mathrm{L}\mathrm{/}\min$ (b).
\begin{overpic}[width=433.62pt]{Figure16.pdf} \end{overpic} Figure 16:
Parameter space of displacement as a function of stiffness, time, and flow
rate ranges of $0$–$100\text{\,}\mathrm{L}\mathrm{/}\min$ (a) and
$0$–$50\text{\,}\mathrm{L}\mathrm{/}\min$ (b).

####  3.3.5 Damping

Both structural and fluid damping is taken into account in the resulting
damping coefficient of the system. This makes it difficult to quantify its
exact value through experiments. In the following, the impact of its typical
parameter ranges on the valve dynamics is investigated. Figure 17 shows a
significant reduction of velocity already at intermediate values, pushing the
orbits into displacement-pressure planes. Moreover, any chattering motion,
which is significant at low damping of
$c<$100\text{\,}\mathrm{N}\text{\,}\mathrm{s}\text{\,}{\mathrm{m}}^{-1}$$,
effectively vanishes at
$c>$200\text{\,}\mathrm{N}\text{\,}\mathrm{s}\text{\,}{\mathrm{m}}^{-1}$$,
even for large volumetric flow rates. It can be recognized in Fig. 18 that the
valve opening is also significantly delayed at higher damping, causing a shift
in displacement amplitudes and opening times. This effect is occurring at all
flow rates, but more severely at the lower range (cf. Fig. 19).

\begin{overpic}[width=433.62pt]{Figure17.pdf} \end{overpic} Figure 17: Orbits
in phase space for various damping coefficients at volumetric flow rates of
$Q=$10\text{\,}\mathrm{L}\text{\,}{\mathrm{min}}^{-1}$$ (a),
$Q=$50\text{\,}\mathrm{L}\text{\,}{\mathrm{min}}^{-1}$$ (b),
$Q=$100\text{\,}\mathrm{L}\text{\,}{\mathrm{min}}^{-1}$$ (c), and
$Q=$500\text{\,}\mathrm{L}\text{\,}{\mathrm{min}}^{-1}$$ (d).
\begin{overpic}[width=433.62pt]{Figure18.pdf} \end{overpic} Figure 18:
Displacement time series for increasing values of the damping coefficient and
flow rates of $0$–$100\text{\,}\mathrm{L}\mathrm{/}\min$ (a) and
$0$–$50\text{\,}\mathrm{L}\mathrm{/}\min$ (b).
\begin{overpic}[width=433.62pt]{Figure19.pdf} \end{overpic} Figure 19:
Displacement in a parameter space of damping coefficient, time, and volumetric
flow rate of $0$–$100\text{\,}\mathrm{L}\mathrm{/}\min$ (a) and
$0$–$50\text{\,}\mathrm{L}\mathrm{/}\min$ (b).

####  3.3.6 Compressibility

A change in compressibility is the most likely reason for strong hysteresis
effects of pressure relief valves in shock absorbers, as described in some
detail by Duym et al. [1997b]. The compressibility of hydraulic oils can
increase due to several reasons such as oil quality, thermal expansion, or
dissolved gas phase in the form of air bubbles. Particularly the presence of
air in the mineral oil can alter the nominal resistance to compression. As
Fig. 20 shows, higher compressibility clearly leads to a damping effect on all
variables, introducing unpredictability and irreversibility to the system.
This is one of the reasons why resulting damping force characteristics of the
shock absorber do not coincide with compression and rebound strokes, despite
comparable valves for both phases.

\begin{overpic}[width=433.62pt]{Figure20.pdf} \end{overpic} Figure 20: Orbits
in phase space for discrete values of compressibility at volumetric flow rates
of $Q=$10\text{\,}\mathrm{L}\text{\,}{\mathrm{min}}^{-1}$$ (a),
$Q=$50\text{\,}\mathrm{L}\text{\,}{\mathrm{min}}^{-1}$$ (b),
$Q=$100\text{\,}\mathrm{L}\text{\,}{\mathrm{min}}^{-1}$$ (c), and
$Q=$500\text{\,}\mathrm{L}\text{\,}{\mathrm{min}}^{-1}$$ (d).

##  4 Conclusions

The basic function of the shim valves within a shock absorber resembles that
of an impact oscillator, a non-smooth dynamical system that shows a multitude
of nonlinear effects such as grazing bifurcations and chaotic motion [Di
Bernardo et al., 2008]. Such effects can be detrimental and can compromise the
regular function of the shock absorber by drastically increased response times
or even complete loss of damping at high excitation frequencies and flow rates
through the valve system. These effects are further complicated by the
presence of multiple phases in the shock absorber (mineral oil in the main
chamber, nitrogen in the gas reservoir). Due to this inherently nonlinear
behaviour, it is of vital importance to understand the impact that parameter
changes have on the dynamic features of a shock absorber pressure relief valve
and on the nature of its motion (i.e. steady, periodic, aperiodic, or
chaotic).

Regarding the main research questions outlined in Sec. 1, the following
findings are documented:

  * 1.

The quantitative value of shock absorber variables, such as valve
displacement, velocity, and pressure, can be computed at its critical points
for any set of boundary conditions and system properties by using the
analytical expressions of Eq. (22) and Eq. (19) for the case without impact
($F_{i}=0$) and Eq. (26) for the case with impact ($F_{i}\neq 0$), where the
fixed points are pseudo-equilibria defined by the bounds of the domain.

  * 2.

It is shown in Sec. 3.3 that while certain parameters, such as the bleed
orifice cross-sectional area have a predictable, linear effect of the motion
of the valve, other parameters, such as the mass and pretension of the spring-
loaded valve show a tendency to destabilize the system under the current
setup. This effect of the pretension force has also been noted by Hős et al.
[2015] with respect to instability mechanisms of spring-loaded pressure relief
valves in gas service. Although shock absorbers operate with a viscous fluid
and are therefore highly dissipative, we could detect instability towards
transient chaos in analogy to this study. An increase in valve mass leads to a
notable occurrence of low to medium frequency tonality in the range of
$400$–$1200\text{\,}\mathrm{Hz}$ and at medium to high volumetric flow rates
of $40$–$500\text{\,}\mathrm{Hz}$ (cf. Sec. 3.3.2). In contrast, the
application of a pretension force leads to sustained periodic and non-periodic
motion (cf. Sec. 3.3.3).

  * 3.

The oscillatory behaviour summarized above is further analysed by using
chaotic indicators, which quantify the level of order and chaos in the system
and provide a measure of the time needed for transient chaos to subside.
Moreover, Fourier transformations alone already give a good indication of the
periodicity of the valve motion. It is shown that intermediate values of
pretension, as mentioned above, lead initially to aperiodic, transient chaotic
valve dynamics, which are eventually damped out due to the high dissipation of
the investigated system. The chaotic indicator SALI is shown to be a potent
method for the detection of transitory chaos in dissipative systems, such as
the one studied here, due to its ability to gauge the immediate dynamical
state, as opposed to the indication of long-term behaviour by the Lyapunov
exponents.

Based on the above conclusions, it becomes clear which quantities and
parameter ranges lead to unsteady, vibratory valve response and to potentially
dangerous loss of function due to instabilities. The results presented in Sec.
3 are a direct extension to the study of the non-dimensional valve model by
Schickhofer and Wimmer [2022], which investigated the exact onset of
instabilities and Hopf bifurcations for realistic flow rates by application of
the Routh-Hurwitz criterion.

Direct conclusions for the effective design of common shock absorber valves
can be drawn from the data presented in this paper: The valve mass should be
kept as low as possible to ensure dynamic stability throughout the entire flow
rate regime. Additionally, if pretension (e.g. by springs) is applied, this
should be carefully adjusted for the operating flow rate, such that no strong
valve oscillations are excited due to pretension and impact rebound force
being of the same order of magnitude. Consequently, by taking into account the
global stability results from Tab. 1 and Tab. 2, one can choose a balanced
setup between mass, pretension, and target flow rate. Finally, the approach
described in Sec. 2 can be applied to other hydraulic valve systems to gain
valuable information on their stability and nonlinear dynamics.

Acknowledgements
This work has been performed under the Project HPC-EUROPA
(INFRAIA-2016-1-730897), with the support of the European Commission for
Research and Innovation under the H2020 Programme (Grant: HPC17ZLPYQ). In
particular, the authors gratefully acknowledge the support of the Department
of Mathematical Sciences at the University of Essex and the computer resources
and technical support provided by the Irish Centre for High-End Computing
(ICHEC).

Conflict of interest statement
The authors have nothing to disclose that would have biased this work.

##  Appendix A Baseline system properties

Below the values for the model quantities of the shock absorber valve are
given. They are based on a realistic pressure relief valve setup for the
purpose of flow control in a suspension system. Wherever different values are
chosen, this is explicitly stated in the text.

Table 1: Properties of the valve system. Parameter | Symbol | Value | Unit
---|---|---|---
Mass | $m$ | $6.8768\times 10^{-4}$ | $\text{\,}\mathrm{kg}$
Spring constant | $k$ | $200$ | $\text{\,}\mathrm{N}\text{\,}{\mathrm{m}}^{-1}$
Impact stiffness | $k_{i}$ | $10^{10}$ | $\text{\,}\mathrm{N}\text{\,}{\mathrm{m}}^{-1}$
Damping coefficient | $c$ | $10$ | $\text{\,}\mathrm{N}\text{\,}\mathrm{s}\text{\,}{\mathrm{m}}^{-1}$
Pretension force | $F_{0}$ | $0.35$ | $\text{\,}\mathrm{N}$
Pressure area | $A_{p}$ | $7.2150\times 10^{-5}$ | $\text{\,}{\mathrm{m}}^{2}$
Bleed orifice section | $A_{b}$ | $2.7\times 10^{-6}$ | $\text{\,}{\mathrm{m}}^{2}$
Discharge coefficient of valve opening | $C_{d,max,v}$ | $0.6$ | $-$
Discharge coefficient of bleed orifice | $C_{d,max,b}$ | $0.6$ | $-$
Momentum coefficient | $C_{f}$ | $0.3$ | $-$
Table 2: Properties of the mineral oil. Parameter | Symbol | Value | Unit
---|---|---|---
Reference density | $\rho_{0}$ | $830$ | $\text{\,}\mathrm{kg}\text{\,}{\mathrm{m}}^{-3}$
Reference compressibility | $\beta_{0}$ | $7.6\times 10^{-10}$ | $\text{\,}{\mathrm{Pa}}^{-1}$
Kinematic viscosity | $\nu$ | $2.3\times 10^{-5}$ | $\text{\,}{\mathrm{m}}^{2}\text{\,}{\mathrm{s}}^{-1}$
Control volume | $V$ | $1\times 10^{-4}$ | $\text{\,}{\mathrm{m}}^{3}$

## References

  * Alligood et al. [1996] Alligood KT, Sauer TD, Yorke JA.  Chaos: An Introduction to Dynamical Systems.  Springer, 1996.
  * Bazsó et al. [2014] Bazsó C, Champneys AR, Hos CJ.  Bifurcation analysis of a simplified model of a pressure relief valve attached to a pipe.  SIAM Journal on Applied Dynamical Systems 2014;13(2):704–21.
  * Benaziz et al. [2012] Benaziz M, Nacivet S, Thouverez F.  Nonlinear dynamic analysis of a shock absorber hydraulic spring valve.  In: International Conference on Noise and Vibration Engineering. 2012. .
  * Benaziz et al. [2015] Benaziz M, Nacivet S, Thouverez F.  A shock absorber model for structure-borne noise analyses.  Journal of Sound and Vibration 2015;349:177–94.
  * Benettin et al. [1980a] Benettin G, Galgani L, Giorgilli A, Strelcyn JM.  Lyapunov characteristic exponents for smooth dynamical systems and for Hamiltonian systems; a method for computing all of them. Part I: Theory.  Meccanica 1980a;15(1):9–20.
  * Benettin et al. [1980b] Benettin G, Galgani L, Giorgilli A, Strelcyn JM.  Lyapunov characteristic exponents for smooth dynamical systems and for Hamiltonian systems; a method for computing all of them. Part II: Numerical application.  Meccanica 1980b;15(1):9–20.
  * Boggs et al. [2010] Boggs C, Ahmadian M, Southward S.  Efficient empirical modelling of a high-performance shock absorber for vehicle dynamics studies.  Vehicle System Dynamics 2010;48(4):481–505.
  * Di Bernardo et al. [2008] Di Bernardo M, Budd CJ, Champneys AR, Kowalczyk P, Nordmark AB, Tost GO, Piiroinen PT.  Bifurcations in nonsmooth dynamical systems.  SIAM Review 2008;50(4):629–701.
  * Dixon [2008] Dixon JC.  The Shock Absorber Handbook.  West Sussex, United Kingdom: John Wiley & Sons, 2008.
  * Duym et al. [1997a] Duym S, Stiens R, Reybrouck K.  Evaluation of shock absorber models.  Vehicle System Dynamics 1997a;27(2):109–27.
  * Duym et al. [1997b] Duym SW, Stiens R, Baron GV, Reybrouck KG.  Physical modeling of the hysteretic behaviour of automotive shock absorbers.  Technical Report; SAE Technical Paper; 1997b.
  * Farjoud et al. [2012] Farjoud A, Ahmadian M, Craft M, Burke W.  Nonlinear modeling and experimental characterization of hydraulic dampers: effects of shim stack and orifice parameters on damper performance.  Nonlinear Dynamics 2012;67(2):1437–56.
  * Hayashi et al. [1997] Hayashi S, Hayase T, Kurahashi T.  Chaos in a hydraulic control valve.  Journal of fluids and structures 1997;11(6):693–716.
  * Hős et al. [2014] Hős C, Champneys A, Paul K, McNeely M.  Dynamic behaviour of direct spring loaded pressure relief valves in gas service: Model development, measurements and instability mechanisms.  Journal of Loss Prevention in the Process Industries 2014;31:70–81.
  * Hős et al. [2015] Hős C, Champneys A, Paul K, McNeely M.  Dynamic behaviour of direct spring loaded pressure relief valves in gas service: Reduced order modelling.  Journal of Loss Prevention in the Process Industries 2015;36:1–12.
  * Hős and Champneys [2012] Hős C, Champneys AR.  Grazing bifurcations and chatter in a pressure relief valve model.  Physica D: Nonlinear Phenomena 2012;241(22):2068–76.
  * Kunze and Küpper [2001] Kunze M, Küpper T.  Non-smooth dynamical systems: an overview.  Ergodic theory, analysis, and efficient simulation of dynamical systems 2001;:431–52.
  * Lang [1977] Lang HH.  A Study of the Characteristics of Automotive Hydraulic Dampers at High Stroking Frequencies.  Ann Arbor, Michigan: University of Michigan, 1977.
  * Liao et al. [2021] Liao M, Zheng Y, Gao Z, Song W.  Fluid–structure coupling modelling and parameter optimization of a direct-acting relief valve for underwater application.  Nonlinear Dynamics 2021;105(4):2935–58.
  * Ma et al. [2019] Ma W, Ma F, Guo R.  Experimental research on the dynamic instability characteristic of a pressure relief valve.  Advances in Mechanical Engineering 2019;11(3):1687814019833531.
  * Mora [2013] Mora K.  Non-Smooth Dynamical Systems and Applications.  Ph.D. thesis; University of Bath; 2013.
  * Reybrouck [1994] Reybrouck K.  A non linear parametric model of an automotive shock absorber.  SAE Transactions 1994;:1170–7.
  * Schickhofer [2022] Schickhofer L.  A universal nonlinear model for the dynamic behaviour of shock absorbers.  Nonlinear Dynamics 2022;:1–23.
  * Schickhofer and Wimmer [2022] Schickhofer L, Wimmer J.  Fluid–structure interaction and dynamic stability of shock absorber check valves.  Journal of Fluids and Structures 2022;110:103536.
  * Segel and Lang [1981] Segel L, Lang H.  The mechanics of automotive hydraulic dampers at high stroking frequencies.  Vehicle System Dynamics 1981;10(2-3):82–5.
  * Sikora [2018] Sikora M.  Modeling and operational analysis of an automotive shock absorber with a tuned mass damper.  Acta Mechanica et Automatica 2018;12(3):243–51.
  * Skokos et al. [2004] Skokos C, Antonopoulos C, Bountis T, Vrahatis M.  Detecting order and chaos in Hamiltonian systems by the SALI method.  Journal of Physics A: Mathematical and General 2004;37(24):6269.
  * Skokos et al. [2003] Skokos C, Antonopoulos C, Bountis TC, Vrahatis MN.  How does the Smaller Alignment Index (SALI) distinguish order from chaos?  Progress of Theoretical Physics Supplement 2003;150:439–43.
  * Talbott and Starkey [2002] Talbott MS, Starkey J.  An experimentally validated physical model of a high-performance mono-tube damper.  SAE Transactions 2002;:2422–39.