## Abstract

A novel modeling approach for viscoelastic hydraulic cylinders, with negligible inertial forces, is proposed, based on the extended fractional-order Jeffreys model. Analysis and physical reasoning for the parameter constraints and order of the fractional derivatives are provided. Comparison between the measured and computed frequency response functions and time domain transient response argues in favor of the proposed four-parameter fractional-order model.

## 1 Introduction

In viscoelasticity, fractional-order models have since long been recognized as being more accurate and flexible than other approaches, in addition to having a lower dimension of parameter space than their integer-order counterparts. Analysis of the various properties of viscoelastic behavior, such as creep, relaxation, viscosity and initial conditions, can be found, e.g., in Refs. [1] and [2]. It is also worth recalling that one of the widespread formulations of the fractional-order differential operator, namely, the Caputo derivative, has been elaborated in the context of a viscoelastic stress–strain relation and memory mechanism associated with the initial conditions (see, e.g., [3] for details).

It is not only viscoelastic solids that have benefited from fractional-order modeling. Soft biological tissues have also been studied experimentally and identified by means of the fractional calculus in viscoelasticity [4]; a more recent review of fractional calculus in modeling biological phenomena can be found in Ref. [5]. Other examples of the complex impact and retardation dynamics, addressed with the help of fractional-order modeling, can be found when analyzing, e.g., backlash [6] or piezoelectric creep [7]. Fractional-order modeling has been shown to be efficient and has also been identified in experiments describing the transient behavior of supercapacitors [8]. Besides the above-mentioned examples, significant efforts in fractional-order calculus of dynamic systems should be credited to fractional-order controllers, though not directly related to the recent study. For exemplary overviews of the tuning of fractional-order controllers for industry applications and fractional proportional-integral-derivative controllers, we refer to, e.g., Refs. [9] and [10], respectively.

The coupled hydrodynamical and mechanical response of hydraulic cylinders, especially with lower acceleration rates, can be challenging when it comes to accurate description by Newton's classical lumped-mass laws of motion. Instead, the viscoelastic behavior of the dashpot type can be seen to be beneficial. The integer-order Jeffreys [11] viscoelastic model takes into account not only the dashpot, which is Newtonian fluid, but also the Kelvin-Voigt viscoelastic structure, which appears to be suitable for mimicking the compressibility of a hydraulic medium. While the classical Jeffreys model is well established and suitable for non-Newtonian fluids, generally allowing for relaxation of stress and strain at different rates, its integer-order dynamics can be suboptimal for different-type hydromechanical couplings. Note that an integer order in the Jeffreys model implies a chain of two integrators, which would restrict the phase response of transfer characteristics by –180 deg at higher frequencies. Involvement of fractional-order dynamics, in contrast, provides additional flexibility in shaping the frequency characteristics of a hydromechanical transducer. At the same time, it allows the dimension of the free model parameters to be kept as low as possible. This is quite the opposite of lumped-parameter rheological modeling approaches, where more sophisticated structures of parallel and serial connections of the viscoelastic elements are involved, in order to approximate the input–output transfer characteristics of interest.

In this brief, we make use of fractional-order viscoelastic behavior to describe the principal dynamics of standard dashpot-type hydraulic cylinders, with a low-parameter and low-order model that can still accurately predict the frequency response characteristics. The proposed solution follows the Jeffreys [11] viscoelastic model (though extending this by an additional stiffness), generalizing to fractional derivatives, and setting physically reasonable constraints on the parameters and differential order. While providing all necessary preliminaries of the fractional-order calculus (see Sec. 2), for the basics on applied hydraulics, we refer to Ref. [12]. The main contribution of this paper is in Sec. 3, and the experimental evaluation and discussion are provided in Sec. 4.

## 2 Fractional Differentiation

In accordance with the seminal literature on fractional-order calculus, e.g., Refs. [1315], we will use a fractional α-order operator on the a and t limits, defined by
$aDtα={dαdtαfor α>01for α=0∫atdt−αfor α<0$
(1)

For the sake of practical relevance, we confine ourselves to real fractional orders, i.e., $α∈ℝ$, often denoted as a noninteger differentiation, correspondingly noninteger integration. For the same reasoning as in practical (engineering) applications, we will consider zero initial time, i.e., a =0.

The classical definition of the Riemann–Liouville α-derivative of a continuous function f(t), with $n−1≤α, where n is an integer, is given by
$0Dtαf(t)=1Γ(n−α)dndtn∫0tf(τ)(t−τ)α−n+1 dτ$
(2)
Here, $Γ(·)$ is the gamma function. We should point out that solving differential equations in terms of Riemann–Liouville derivatives (2) requires the initial conditions
$[0Dtα−kf(t)]t→0=ck for k=1,2,…,n$
(3)
to be known and determined accordingly. That is, one needs the initial values of the $(α−k)$-fractional derivatives of the function f(t). This is particularly evident when considering the Laplace transform $F(s)=L{f(t)}$ of the Riemann–Liouville fractional derivative (cf. Ref. [15]), given by
$L{0Dtαf(t)}=sαF(s)−∑k=0n−1sk[0Dtα−k−1f(t)]t=0$
(4)
The Grünwald–Letnikov definition (cf. Ref. [15])
$0Dtαf(t)=limh→0 h−α∑i=0[th−1](−1)i(αi)f(t−ih)$
(5)
of the fractional-order derivative is valid for any $α∈ℝ$ and is particularly suitable for numerical implementations, since it constitutes the limit of the difference quotient
$Δhαf(t)≈ 0Dtαf(t)$
(6)
with the time-step $h→0$. In (5), the operator $[x]$ means the integer part of x, while i is the index of the discrete time series of t. The binomial coefficients, which are sign-alternating and summarized as
$wi(α)=(−1)i(αi) for i=0,1,2,…$
(7)
can be evaluated recursively (cf. Ref. [15]) by
$w0(α)=1; wi(α)=(1−α+1i)wi−1(α) for i=1,2,3,… .$
(8)
Further, we will also make use of the fact that for zero initial conditions, the Laplace transform is given by
$L0{0Dtαf(t)}=sαF(s)$
(9)
and we will write the Fourier transform as
$F{0Dtαf(t)}=(jω)αF(ω)$
(10)

Note that for $F$, the magnitude and phase responses are conventionally given as in the case of $α∈Z$.

## 3 Modeling

Our starting point is the standard dashpot element used in viscoelasticity for analyzing and representing an ideal linear viscous fluid. A constantly applied stress σ produces a constant strain rate so that
$σ=μddtε$
(11)
where μ is the coefficient of viscosity of a Newtonian fluid. Note that a linear dashpot, as a basic mechanical element, is approximating a simple piston cylinder, with one degree-of-freedom x, where the stress equivalent of the hydraulic pressure is replaced by the applied force τ. Correspondingly, the strain rate is equivalent to the rate of displacement so that, in the Laplace domain, one obtains
$x=μ−1s−1τ$
(12)

Obviously, Eq. (12) constitutes a free integrator dynamic, this way yielding a simple first-order lumped model of the viscous driven motion when inertial effects are neglected. Note that the latter is justified for multiple hydraulic cylinders, deployed as the actuators, in which the viscoelastic forces largely dominate over the inertial.

Next, we elaborate whether the hydraulic cylinder force, induced by the pressure difference, itself undergoes the dynamic transients to be captured by viscoelastic behavior. The stress–strain relation of the standard linear solid model, also denoted as Zener model, is given by
$E[1ϕddt+1]ε(t)=[1φddt+1]σ(t)$
(13)

where E is a suitable elastic modulus. The positive constants $ϕ−1>φ−1$ refer to the retardation and relaxation times (see, e.g., Ref. [11]). It can be noted that $ϕ=φ$ simplifies (13) to a purely elastic material, i.e., $σ=Eε$, while Eq. (13) includes equally the Maxwell and Kelvin-Voigt models of the viscoelastic fluid, correspondingly solid, as limiting cases [3]. The otherwise unequal and nonzero time constants of retardation and relaxation shape the frequency characteristics of a standard linear solid within $ϕ<ω<φ$, while keeping it at different constant levels for $ω→0$ and $ω→∞$ limits. Since the dynamic behavior (13) has zero relative degree, its input–output transfer characteristics can also be considered as entirely unitless. For approaching a viscoelastic fluid in the closed hydraulic circuit (of a powered cylinder) as a quasi-solid medium, the model (13) can be used to describe its dynamic compressibility. Thus, the applied stress σ appears in equivalence to the supplied pressure difference, while the strain ε mimics a stiff kinematic excitation, and therefore force, acting on the piston interface. In this equivalence, the elastic modulus appears as the bulk modulus [12] of the hydraulic medium. The retardation time approaches the time constant of the effective force response to a stepwise change of the differential pressure. Similarly, the relaxation time relates to the time constant of the pressure response to an instantaneous change in the piston stroke and thus force on the piston interface.

### 3.1 Jeffreys Model.

The Jeffreys model [11] comprises, essentially, the Kelvin-Voigt and dashpot elements connected in series (see structural arrangement in Fig. 1). The overall strain (or relative displacement) rate is a superposition of both
$ddtε=ddtε1+ddtε2$
(14)
while the stress (corresponding force) in both elements is the same, meaning
$σ=Kε1+μ1ddtε1=μ2ddtε2$
(15)
Fig. 1
Fig. 1
Close modal
Combining Eqs. (14) and (15), and eliminating $ε1$ and $ε2$, results in the Jeffreys model of the form (see Ref. [11] for details)
$μ2[μ1Kd2dt2+ddt]ε(t)=[μ1+μ2Kddt+1]σ(t)$
(16)

Here, the relaxation and retardation times are given by $(μ1+μ2)K−1$ and $μ1K−1$, correspondingly. Note that, as in the case of the Zener model, elimination (set to zero) of μ1 or μ2 simplifies the Jeffreys model (16) to the Maxwell or Kelvin–Voigt model, respectively.

### 3.2 Fractional-Order Formulation.

Allowing for the fractional-order derivatives and considering the input–output pair of differential pressure (equivalent to hydraulic force) and relative displacement, correspondingly, the fractional-order Jeffreys model, in the Laplace domain, is written as
$μsγ(λ2sα+1)x(s)=(λ1sβ+1)τ(s)$
(17)

Note that a similar modified Jeffreys model has also been proposed, though for viscoelastic fluids only, in Ref. [16]. It can be recognized that for the same time constants $λ1=λ2$ and fractional orders $α=β$, the proposed model (17) is reduced to the simple viscous dashpot, cf. with Eq. (12).

Unlike the original Jeffreys model (cf. Ref. [11]) and the modified one [16], the inequality $λ2>λ1$ of parameters is required for Eq. (17). Note that this implies the lag characteristics of the $x(s)/τ(s)$ transfer function, which means a phase-lowering region and amplitude drop within $λ2−1<ω<λ1−1$. From a structural viewpoint, it would require an additional spring that would lead to a serial connection of the Zener model with a viscous dashpot (see above in Sec. 3), while the dynamics order and relative degree remain unchanged this way. The viscoelastic damping properties of the Zener model (13) argue in favor of the lag characteristics of Eq. (17). Indeed, the force propagation through the hydraulic medium, from the differential pressure source to the effective piston force, is weakened for higher frequencies and additionally lagged for a certain frequency range. Here, we recall that for the Zener model to be dissipative and, therefore, physically reasonable for real materials (and corresponding media), the thermodynamic constraints of the $α=β$ parameters should be additionally satisfied. This has been observed previously in multiple rheological studies and also theoretically proved in Ref. [17].

While the general form Eq. (17) allows for all derivatives to be fractional-order, we further require γ = 1 to avoid violating causality and physical reasoning of the force–displacement transfer characteristics. This is analyzed below in the context of the impulse response and its steady-state (final) value. The above physical constraints on the time constants and differential orders lead us to the overall four-parameter linear fractional-order model (17) of viscoelastic hydraulic cylinders with negligible inertial terms.

### 3.3 Initial Value and Impulse Response.

Next, we need to clarify the initial conditions for fractional differential Eq. (17). The initial condition for the integer-order dashpot, which is the relative displacement at t =0, can be assumed to be zero without loss of generality. The remaining fractional-order dynamics of the Zener model require the single initial condition $0Dtα−1x(t)$ (cf. Ref. [1]), provided the input differential pressure (corresponding actuation force) is a known exogenous value. It has been shown in Ref. [1] that for a physically reasonable (i.e., continuous) loading, or even in the case of a step discontinuity, zero initial conditions for the strain dynamics apply. Nonzero initial conditions will only be valid in the case of a $B δ(t)$ stress impulse (cf. Ref. [1]), resulting in
$[λ2 0Dtα−1x(t)]t→0=B$
(18)
We now have to take a closer look at the impulse response of Eq. (17), in terms of the final value problem from which the γ = 1 constraint is enlightened. Rewriting (17) in the time domain and integrating both sides with respect to $dtγ$ results in
$μλ2 0Dtαx(t)+μx(t)= 0Dt−γ(λ1 0Dtβτ(t)+τ(t))$
(19)
The steady-state of Eq. (19) can be obtained via the final value theorem, since for $α,β>0$, both differential terms in Eq. (19) vanish for $t→∞$. Solving the final value problem for x(t) this way yields
$[x(t)]t→∞=μ−1∫0tτ(t)dtγ$
(20)
For the applied Dirac impulse $δ(t)$, the integral solution (20) can be directly evaluated (see, e.g., Ref. [15]), given as
$0Dtηδ(t)=t−η−1Γ(−η)$
(21)

It can be shown that for $η<0$, the value of Eq. (21) remains, $∀ t>0$, constant and equal to one for $η=−1$ only. Otherwise, it converges to zero for $−1<η<0$ and diverges for $η<−1$. For the final value (20), it means that a constant nonzero impulse response of the relative displacement can only be achieved when γ = 1. Indeed, the free integrator in Eq. (17) requires a constant finite displacement as a result of a force impulse which has finite energy content. An illustrative numerical example of the Eq. (17) model response to the input impulse is shown in Fig. 2 for $γ={0.9, 1, 1.1}$.

Fig. 2
Fig. 2
Close modal

## 4 Experimental Evaluation and Discussion

Experimental evaluation of the above model was made with the data recorded from a standard linear-stroke hydraulic cylinder in a laboratory setting (see Fig. 3). The standard one-side-rod hydraulic cylinder, with a full cross section of 25 mm, is actuated via a 4/3 servovalve with the operational supply pressure set to 100 bar. The relative displacement x(t) is directly measured by an absolute resistor-based linear potentiometer, while the total effective stroke is about 200 mm. Two pressure sensors, with a measurement range up to 250 bar, are installed on the fittings, close to both cylinder chambers. The differential pressure $τ(t)$ is directly obtained from the measurements. All experimental data are generated in an open-loop control manner by designing and feeding the reference control signal for the servovalve for both frequency domain and time domain measurements. The frequency domain data, used below for identification, were collected as single steady-state points for the set of harmonic excitations, and that using averaged sine waves of different frequencies between 0.005 Hz and 1.6 Hz. Further details of the experimental setup can be found in [18].

Fig. 3
Fig. 3
Close modal

### 4.1 Frequency Response Function.

The measured frequency response function $G(jω)=x(jω)/τ(jω)$ is used for the least-squares fit of the fractional-order model (17). Note that both the amplitude response in dB and phase response in deg have been incorporated into the objective function. Given the comparable range of dB and deg units, the amplitude and phase response errors are equally included in the overall objective function of numerical minimization. For the sake of comparison, the least-squares best fit has also been found (from the same experimental data) for the integer-order model with the structure as in Eq. (17). The measured and both-ways modeled frequency response functions are shown opposite to each other in Fig. 4. The fractional-order model coincides with the measurements, while the integer-order one misses the phase response. The identified parameters of the proposed model (17) are $λ1=0.013, λ2=0.047$, $α=β=1.571, μ=171e3$.

Fig. 4
Fig. 4
Close modal

### 4.2 Time Series.

A similar observation can be made when comparing the time series of the identified fractional- and integer-order models with the initial phase of the measured response to the slope-shaped input $τ(t)$. From Fig. 5, one can see that only the fractional-order model captures a lagged transient of the relative displacement. Note that further progress of the measured response (not shown in the figure) diverges from each model in all cases, due to an inherent free integrator error. Yet the fractional-order model reproduces the initial transition of the viscoelastic response to a linearly increasing input force.

Fig. 5
Fig. 5
Close modal

### 4.3 Concluding Remarks.

The fraction-order Jeffreys model manifests as more accurate, the opposite of its integer-order counterpart, when describing the dynamic behavior of a viscoelastic hydraulic cylinder. This becomes particularly visible in the frequency domain for the phase response and, as a logical consequence, within the transient shape in the time domain. The identified four-parameter model (17) takes into account the physical reasoning of force–displacement transfer characteristics and thermodynamic constraints.

## Acknowledgment

This research received funding from the IS-DAAD program under RCN project number 294835. Laboratory measurements performed by Daniela Kapp are also acknowledged.

## Funding Data

• Norges Forskningsråd (Grant No. 294835; Funder ID: 10.13039/501100005416).

## References

1.
Heymans
,
N.
, and
Podlubny
,
I.
,
2006
, “
Physical Interpretation of Initial Conditions for Fractional Differential Equations With Riemann-Liouville Fractional Derivatives
,”
Rheol. Acta
,
45
(
5
), pp.
765
771
.10.1007/s00397-005-0043-5
2.
Mainardi
,
F.
, and
,
G.
,
2011
, “
Creep, Relaxation and Viscosity Properties for Basic Fractional Models in Rheology
,”
Eur. Phys. J. Spec. Top.
,
193
(
1
), pp.
133
160
.10.1140/epjst/e2011-01387-1
3.
Caputo
,
M.
, and
Mainardi
,
F.
,
1971
, “
A New Dissipation Model Based on Memory Mechanism
,”
Pure App. Geop.
,
91
(
1
), pp.
134
147
.10.1007/BF00879562
4.
Meral
,
F.
,
Royston
,
T.
, and
Magin
,
R.
,
2010
, “
Fractional Calculus in Viscoelasticity: An Experimental Study
,”
Commun. Nonlinear Sci. Numer. Simul.
,
15
(
4
), pp.
939
945
.10.1016/j.cnsns.2009.05.004
5.
Ionescu
,
C.
,
Lopes
,
A.
,
Copot
,
D.
,
,
J. T.
, and
Bates
,
J.
,
2017
, “
The Role of Fractional Calculus in Modeling Biological Phenomena: A Review
,”
Commun. Nonlinear Sci. Numer. Simul.
,
51
, pp.
141
159
.10.1016/j.cnsns.2017.04.001
6.
,
J. T.
,
2013
, “
Fractional Order Modelling of Dynamic Backlash
,”
Mechatronics
,
23
(
7
), pp.
741
745
.10.1016/j.mechatronics.2013.01.011
7.
Liu
,
Y.
,
Shan
,
J.
, and
Qi
,
N.
,
2013
, “
Creep Modeling and Identification for Piezoelectric Actuators Based on Fractional-Order System
,”
Mechatronics
,
23
(
7
), pp.
840
847
.10.1016/j.mechatronics.2013.04.008
8.
Freeborn
,
T. J.
,
Maundy
,
B.
, and
Elwakil
,
A. S.
,
2013
, “
Measurement of Supercapacitor Fractional-Order Model Parameters From Voltage-Excited Step Response
,”
IEEE J. Emer. Sel. Top. Circuits Syst.
,
3
(
3
), pp.
367
376
.10.1109/JETCAS.2013.2271433
9.
Monje
,
C. A.
,
Vinagre
,
B. M.
,
Feliu
,
V.
, and
Chen
,
Y.
,
2008
, “
Tuning and Auto-Tuning of Fractional Order Controllers for Industry Applications
,”
Control Eng. Prac.
,
16
(
7
), pp.
798
812
.10.1016/j.conengprac.2007.08.006
10.
Shah
,
P.
, and
Agashe
,
S.
,
2016
, “
Review of Fractional PID Controller
,”
Mechatronics
,
38
, pp.
29
41
.10.1016/j.mechatronics.2016.06.005
11.
Joseph
,
D. D.
,
1990
,
Fluid Dynamics of Viscoelastic Liquids
,
Springer
,
Berlin
.
12.
Jelali
,
M.
, and
Kroll
,
A.
,
2012
,
Hydraulic Servo-Systems: Modelling, Identification and Control
,
Springer
,
London
.
13.
Oldham
,
K.
, and
Spanier
,
J.
,
1974
,
The Fractional Calculus Theory and Applications of Differentiation and Integration to Arbitrary Order
,
Elsevier
,
New York and London
.
14.
Samko
,
S. G.
,
Kilbas
,
A. A.
, and
Marichev
,
O. I.
,
1993
,
Fractional Integrals and Derivatives
,
Gordon and Breach Science
,
Amsterdam, The Netherlands
.
15.
Podlubny
,
I.
,
1998
,
Fractional Differential Equations
,
Elsevier
,
San Diego/London, CA
.
16.
Song
,
D. Y.
, and
Jiang
,
T. Q.
,
1998
, “
Study on the Constitutive Equation With Fractional Derivative for the Viscoelastic Fluids–Modified Jeffreys Model and Its Application
,”
Rheol. Acta
,
37
(
5
), pp.
512
517
.10.1007/s003970050138
17.
Bagley
,
R. L.
, and
Torvik
,
P. J.
,
1986
, “
On the Fractional Calculus Model of Viscoelastic Behavior
,”
J. Rheol.
,
30
(
1
), pp.
133
155
.10.1122/1.549887
18.
Pasolli
,
P.
, and
Ruderman
,
M.
,
2018
, “
Linearized Piecewise Affine in Control and States Hydraulic System: Modeling and Identification
,”
IEEE 44th Annual Conference of the Industrial Electronics Society
, Washington, DC, Oct. 21–23, pp.
4537
4544
.10.1109/IECON.2018.8591572