Efficiency comparison of various friction models of a hydraulic cylinder in the framework of multibody system dynamics

Dynamic simulation of mechanical systems can be performed using a multibody system dynamics approach. The approach allows to account systems of other physical nature, such as hydraulic actuators. In such systems, the nonlinearity and numerical stiffness introduced by the friction model of the hydraulic cylinders can be an important aspect to consider in the modeling because it can lead to poor computational efficiency. This paper couples various friction models of a hydraulic cylinder with the equations of motion of a hydraulically actuated multibody system in a monolithic framework. To this end, two static friction models, the Bengisu–Akay model and Brown–McPhee model, and two dynamic friction models, the LuGre model and modified LuGre model, are considered in this work. A hydraulically actuated four-bar mechanism is exemplified as a case study. The four modeling approaches are compared based on the work cycle, friction force, energy balance, and numerical efficiency. It is concluded that the Brown–McPhee approach is numerically the most efficient approach and it is well able to describe usual friction characteristics in dynamic simulation of hydraulically actuated multibody systems.


Introduction
Dynamic simulation of mechanical systems can be performed using multibody system dynamics. In this approach, the equations of motion describe a force equilibrium for the mechanical system under consideration. This approach allows to describe systems of other physical nature, such as hydraulic actuators [32,33,43]. Hydraulic systems are highly nonlinear systems, which are caused, in part, by the friction model of the hydraulic cylinders. The friction model of the hydraulic cylinders can introduce numerical stiffness [38], and as a consequence, time integration of hydraulically actuated systems may become cumbersome, especially, in real-time applications [30,31,50]. Moreover, a physically incorrect friction model of the hydraulic cylinders can lead to inaccurate simulations and control errors.
Friction models can, in general, be categorized into two main groups, namely, static friction models and dynamic friction models [5,37,45]. Static friction models describe only the steady-state relation between the friction force and the relative velocity [5]. The wellknown friction model was developed by Coulomb in 1785, and over the years, various friction models have been proposed that capture Coulomb, stiction, viscous, and Stribeck [29] friction effects, as in [14,27]. The exponential Stribeck friction model described by Bo and Pavelescu [14] is popular, however, it can lead to numerical problems because of the discontinuity of the friction force at zero relative velocity. To overcome this numerical problem, a separate equation can be used to describe friction at zero relative velocity [2][3][4]10]. For example, Bengisu and Akay [10] utilize two separate equations for the adhesion regime and Stribeck effect. However, when this approach is adopted, the use of different equations during the simulation can lead to poor numerical efficiency. To address this problem, Brown and McPhee [16] recently developed a continuous velocity-based friction model suitable for computationally challenging applications and optimization problems. Their model [16] does not consider the time lag [11], micro-slip [2] or pre-slip displacement [11], and dwell-time dependence of friction [23]. In this study, the friction models developed by Bengisu and Akay [10] and Brown and McPhee [16] are referred to as the Bengisu-Akay and Brown-McPhee friction models, respectively.
Dynamic friction models evaluate the friction force based on the actual state of contact and the contact history [5,38]. In general, extra state variables, such as the average bristle deflection, are used with relative velocity to model the evolution of the friction force [37]. Unlike most static friction models, dynamic friction models consider frictional phenomena such as preslip displacement [11] and frictional lag, and dynamic friction models eliminate the discontinuity at zero relative velocity. The Dahl friction model [20] was the first dynamic friction model in which the friction force is considered proportional to the average bristle deflection. It did not, however, consider the static friction. Based on the Dahl friction model, a number of friction models, such as the reset integrator model [24], Bliman-Sorine model [12,13], LuGre model [21], Leuven model [51], and generalized Maxwell slip model [1,36], are described accounting for various frictional phenomena. The LuGre friction model [21] is widely used and it captures the Stribeck and stiction effects, as well as the pre-slip displacement and frictional lag. Subsequently, a modified LuGre friction model [57] was developed that considered the dynamics of the lubricant film formation. In the modified LuGre friction model, the lubricant film dynamics is computed using a first-order lag element in which the time constant is varied between the deceleration, acceleration, and dwell periods [57]. It should be noted that dynamic friction models can introduce numerical stiffness to systems and, consequently, increasing the computational cost [38].
In the literature, experiments have been conducted to measure the friction force of hydraulic cylinders with different reciprocating seals operating under various conditions [42]. Additionally, a number of studies have been carried out that utilize different seal friction models for simulation of hydraulic cylinders. For example, the exponential Stribeck friction model [14] was used in [53,54] and the stick-slip friction model [3,4] in [39]. Furthermore, the dynamic friction behavior of hydraulic cylinders was studied using the LuGre friction model in [41,53] and a modified LuGre friction model [57] in [52,53]. Yanada and Sekikawa [57] found that compared with the LuGre friction model, the modified LuGre friction model [57] can simulate the real friction characteristics of a hydraulic cylinder with a good accuracy. Other work has used simple pressuredependent friction models [15,25] that either had discontinuity at zero relative velocity [25] or were casedependent [15]. Ylinen et al. [58] combined a pressuredependent function with the LuGre friction model to study multibody simulation of a hydraulic cylinder.
In the framework of multibody system dynamics [35,46], a large number of studies have been presented that demonstrate coupling of multibody systems and hydraulic actuators [22,40]. Coupling has been carried out using either monolithic [48,49] or cosimulation approaches [44,47]. There are even studies using a monolithic approach that demonstrate this coupling for real-time applications, such as [6] and [31]. Despite widespread interest in the field, detailed comparative study on the inclusion of various friction models of a hydraulic cylinder with the equations of motion of hydraulically actuated multibody systems has been overlooked in the literature.
The objective of this paper is to couple various friction models of a hydraulic cylinder with the equations of motion of a hydraulically actuated multibody system in a monolithic framework. To this end, the Bengisu-Akay model [10] and the Brown-McPhee model [16], in the category of static friction models, and the LuGre model [21] and the modified LuGre model [57], in the category of dynamic friction models, are considered in this work. In the monolithic coupling scheme, the multibody system is described using a semi-recursive formulation [8,18] and the hydraulic system is described using the lumped fluid theory [56]. In the study, only a fixed step-size integrator is considered because the interest of the authors lies in the area of real-time simulation. The coupled systems are referred to with the name of the friction force model used, that is, the Bengisu-Akay approach, Brown-McPhee approach, LuGre approach, and modified LuGre approach. Although a case study of a hydraulically actuated planar four-bar mechanism is shown, the methods presented in this study are general and applicable to three-dimensional mechanisms as well. Using the numerical example, the four approaches are compared based on the simulation work cycle, friction force, energy balance, and numerical efficiency.

Multibody system dynamics
In this study, the dynamics of a constrained mechanical system is described using a semi-recursive multibody formulation. In this method, the equations of motion of a system are formulated in the relative joint coordinates using the dynamics of an open-loop system. Closed-loop systems are modeled by employing the cut-joint constraints using the penalty-based augmented Lagrangian method [8,17,18]. The multibody system is hydraulically actuated and the dynamics of the hydraulic system is described using the lumped fluid theory [56]. It should be noted that the hydraulic cylinder-piston is not considered as separate bodies, and therefore, their masses are assumed to be neglectable. Nevertheless, the internal dynamics of the hydraulic system, including the friction model of a hydraulic cylinder, are computed and the resultant force is fed into the multibody equations of motion in a monolithic approach.

Semi-recursive multibody formulation
Consider an open-loop system, which has N b bodies. The absolute velocity, Z j , and acceleration,Ż j , of the bodies can be described using Cartesian velocities and accelerations [28,34] as Here,ṡ j ands j are, respectively, the velocities and accelerations of a point on body j ∈ [1, N b ], which are instantaneously coincident with the origin of the inertial reference frame, and ω j anḋ ω j are their angular counterparts. The kinematics for an open-loop system, shown in Fig. 1a, can be calculated, for example, from the base to the leaves using classical kinematic relations [19]. In Fig. 1, the relative joint coordinates, z j , represent the position of the bodies. The absolute velocities and accelerations can be recursively expressed with respect to the previous bodies as [34] whereż j andz j are, respectively, the relative joint velocities and accelerations, and b j and d j are vectors that depend on the joint-type connecting bodies j − 1 and j [19]. As the system may branch, the indexes j −1 and j may not be successive.
, into a set of relative joint velocities,ż = ż 1 ,ż 2 , . . . ,ż N b T , as [34] where T ∈ R 6N b ×6N b is the constant path matrix that represents the topology of the open-loop system and R d ∈ R 6N b ×N b is a diagonal matrix whose elements are vectors b j , arranged in an ascending order. In Eq. (4), the termṘż can be expressed in terms of vectors d j using Eq. (2) [19].
In an open-loop system, the virtual power of the inertia and external forces can be written as [34] where an asterisk (*) denotes the virtual velocities, and the mass matrixM j and force vectorQ j can be expressed as and where m j and J j are, respectively, the mass and inertia tensor of body j, I 3 is a (3 × 3) identity matrix, the vector g j is the position of the center of mass of body j, a tilde (∼) denotes the skew-symmetric matrix of a vector, and τ j and f j are, respectively, the vector of external moments with respect to the center of mass and external forces applied on body j. By substituting Eqs.
wherez ∈ R N b is the vector of the relative joint accelerations, andM ∈ R 6N b ×6N b andQ ∈ R 6N b are, respectively, a diagonal matrix and a column vector, whose elements areM j andQ j . By following the penalty-based augmented Lagrangian method [17,18], a set of m loop-closure constraints, = 0, are incorporated in Eq. (8) and the equations of motion for the closed-loop system can be written as [18,49] where α is a penalty factor that can be set to the same value for all constraints, z is the Jacobian matrix of (z) = 0, and λ is the vector of Lagrange multipliers that are iterated at each time-step (k) as where h is the iteration step. The value of λ (0) k is the final value of λ k−1 from the previous time-step. For simplicity, the constraints are assumed to be holonomic. An example of a closed-loop system is shown in Fig. 1b.
In this study, an implicit single-step trapezoidal rule as in [7] is used. In this method, the relative joint velocities,ż, and accelerations,z, are corrected using massdamping-stiffness-orthogonal projections as [9,17] whereż andz are, respectively, the relative joint velocities and accelerations obtained from the Newton-Raphson iteration, t is the time-step, W = M + t 2 C + t 2 4 K, with C and K representing, respectively, the contributions of damping and stiffness in the system, and t is the partial derivative of the constraints with respect to time t. In this study, the lumped fluid theory [56] is used to describe the hydraulic pressures in a hydraulic circuit. In this approach, a hydraulic circuit is partitioned into discrete volumes where the pressures are assumed to be equally distributed. Thus, the effect of acoustic waves is assumed to be negligible. The modeling concept is illustrated in Fig. 2, where one incoming flow, Q k1 , and two outgoing flows, Q k2 and Q k3 , are associated with a volume, V k . The pressure, p k , within a hydraulic volume, V k , can be expressed as a function of the effective bulk modulus, B e k , volume, V k , and incoming and outgoing flows using the following differential equatioṅ where Q ks is the sum of incoming and outgoing flows related to volume V k , and n f is the total number of volume flows going in or out of the volume V k . The effective bulk modulus, B e k , associated with the volume V k can be expressed as where B oil is the bulk modulus of oil, B s is the bulk modulus of the volume V s that forms the volume V k , and n c is the total number of sub-volumes V s .

Throttle and directional control valves
In this study, a semi-empirical modeling approach [26] is used to describe the valves in the system. In this approach, the volume flow rate, Q t , through a simple throttle valve can be expressed as where C v t is the semi-empirical flow rate coefficient of the throttle valve, p is the pressure difference over the valve, and sgn( p) is the signum function which determines the sign of p. C v t can be calculated as where ρ is the oil density, A t is the area of the throttle valve, and C d is the flow discharge coefficient. Similarly, the volume flow rate, Q d , through a directional control valve can be expressed as where C v d is the semi-empirical flow rate constant of the valve, and U is the relative spool position. The value of C v d can be obtained from the manufacturers' catalogues. The volume flow is assumed to be laminar for pressure differences of less than 2 bar, and Eqs. (15) and (17) are modified such that the volume flow and pressure difference follow a linear relation. The spool position can be described using the following differential equatioṅ where U ref is the reference voltage signal for the reference spool position, and τ is the time constant. The value of τ can be obtained from the Bode-diagram of the valve, which describes the dynamics of the valve spool.

Double-acting hydraulic cylinder
The motion of a double-acting hydraulic cylinder, shown in Fig. 3, causes volume flows that can be written as where Q in and Q out are, respectively, the incoming and outgoing volume flow rates of the cylinder,ẋ is the piston velocity, and A 1 and A 2 are, respectively, the piston and piston-rod areas of the cylinder. The force, F cyl , produced by a cylinder can be written in terms of where p 1 and p 2 are, respectively, the pressures on the piston and piston-rod sides of the cylinder, and F μ is the total friction force caused by sealing.

Friction force models
Friction force is defined as a force that opposes the relative motion between two surfaces in contact. In the literature, several friction force models have been proposed. These models capture a combination of the Coulomb, stiction, viscous, and Stribeck effects. However, most of the friction models in the static friction model category may suffer from numerical instability because of the discontinuity of the friction force at zero relative tangential velocity. To overcome this numerical inconvenience, many alternative continuous friction models have been developed, for example, as shown in Fig. 4.
In the category of static friction models, the Bengisu-Akay and Brown-McPhee friction models [10,16] are considered in this work. The friction models in the dynamic friction category consider the pre-slip displacement [11]. It is for this reason that these models eliminate the discontinuity of the friction force at zero relative tangential velocity. In general, they consider the evolution of friction force based on the actual state of contact and the contact history [38]. In a dynamic model, an extra state variable, such as the average bristle deflection [21], is used, which represents the behavior of surface asperities during the contact [38]. During the sticking phase, the bristles behave like springs. In the category of dynamic friction models, the LuGre and modified LuGre friction models [21,57] are considered in this work.

Bengisu-Akay friction model
An approximation of the friction force can be presented by the Coulomb friction law, which established that the magnitude of the friction force depends on the relative tangential velocity between the surfaces in contact. Furthermore, it has been established that the static friction is greater than the kinetic or Coulomb friction and it decreases continuously with the increment of relative tangential velocity. This effect is known as the Stribeck effect [29], as shown in Fig. 4. By considering the numerical stability at zero relative tangential velocity, the friction force can be modeled using two equations, one for the slope and the other for the Stribeck effect, as [10,37] where F c is the Coulomb friction, F s is the static friction, F μ is the total friction force, v s is the Stribeck velocity, v t is the relative tangential velocity between the two surfaces in contact, and ξ is a positive parameter that represents the negative slope of the sliding regime. Note that F c = μ c |F n | and F s = μ s |F n |, where μ c and μ s are the coefficients of Coulomb and static frictions, respectively, and F n is the normal contact force. It should be noted that the viscous friction is not considered in Eq. (21). However, for a hydraulic cylinder where lubricants are present, as in this study, the viscous friction cannot be ignored. Therefore, a corrected version of the Bengisu-Akay friction model is presented in this study by incorporating the viscous friction as where the function f (v t ) describes the viscous friction.
In lubricated contacts, viscous friction is generally considered to be linear viscous friction as [21,38] f Fig. 4 Representation of the classical Stribeck friction model (left [14]), and two alternative friction models (middle [37] and right [10,16]) to overcome discontinuity at zero relative tangential velocity where σ 2 is the coefficient of viscous friction. The effect of viscous friction is crucial, especially, in the presence of fluid lubricant or at high relative velocity [37]. A representation of this model is shown on the right side of Fig. 4. It should be noted that if the slope at v t = 0 in Eq. (22) is large, small time-steps are required to properly capture the friction effects at low velocity. This may burden the computational efficiency.

Brown-McPhee friction model
A suitable alternate method to the classical Stribeck friction model that avoids the discontinuity at zero relative tangential velocity has been presented by Brown and McPhee [16]. This model incorporates the Coulomb, stiction, and viscous friction, and it is valid for both positive and negative relative tangential velocity. This velocity-dependent continuous differential model can be written as [16] In Eq. (24), the Coulomb friction term approaches F c as |v t | increases from zero to v s . Similarly, the static friction term, which represents the stiction and Stribeck effect, has its maximum value as (F s − F c ) at |v t | = v s and it then decreases to zero for larger relative tangential velocity, where the Coulomb and viscous friction is dominant. Here, the viscous friction term is linearly scaled with velocity and physically corresponds to thin-film viscous friction [16]. A representation of this model is shown on the right side of Fig. 4. This

LuGre friction model
In LuGre friction model, the friction is considered as an interaction of the bristles between two surfaces, as shown in Fig. 5, and the model utilizes the average bristle deflection. The LuGre (Lund-Grenoble) model [21] is an extension of the Dahl model [20] that captures the Stribeck and stiction effects. When a tangential force is applied, the bristles start to deflect like a spring, giving rise to the sticking phase. The moment the tangential force is sufficiently large, some of the bristles deflect so much that they start to slip [21].
The average bristle deflection is quantified by introducing an internal state variable, q, such that the model can be expressed as [21,37] with where q is the average bristle deflection,q is the rate of average bristle deflection, σ 0 is the stiffness of the elastic bristles, the function g(v t ) describes the Stribeck effect [21,53], and n is the exponent that affects the slope of the Stribeck curve. The total friction force, F μ , generated from the bending of the bristles can be written as where σ 1 is the damping of the elastic bristles, which can be a constant or a function of v t .

Modified LuGre friction model
In this model, the LuGre friction model (Sect. 3.3) is modified by incorporating the dynamics of the lubricant film formation [57]. The lubricant film dynamics is computed using a first-order lag element in which the time constant is varied between the deceleration, acceleration, and dwell periods [57]. Thus, in this model, the friction is obtained as a combination of the bristle model (Sect. 3.3) and the lubricant film dynamics. Similar to the LuGre friction model, the average bristle deflection is quantified by an internal state variable, q, such that [57] with where h is the dimensionless lubricant film thickness parameter, and the function g(v t , h) describes the Stribeck effect. The lubricant film dynamics can be written as [52,57] with where τ hp , τ hn , τ h0 are, respectively, the time constants for the acceleration, deceleration, and dwell periods, h ss is the dimensionless steady-state lubricant film thickness parameter, K f is a proportional constant, and v b is the velocity at which the magnitude of the steady-state friction force becomes almost minimum [52,57]. In Eq. (31), |h| ≤ |h ss | represents the acceleration period and |h| > |h ss | represents the deceleration period. Here, it is assumed [57] that the lubricant film thickness increases with velocity only in the negative resistance regime, that is, |v t | ≤ |v b |, otherwise, it has a constant maximum value, as shown in Eq. (32). The total friction force, F μ , generated from the bending of the bristles and lubricant film dynamics can be written as

Coupling the dynamics of the multibody and hydraulic systems
In this section, the multibody formulation in Sect. 2.1 is extended to incorporate the dynamics of the hydraulic system, introduced in Sect. 2.2, with the friction force models of a hydraulic cylinder, introduced in Sect. 3, in a monolithic approach [32,49]. The force vector, Q , in Eq. (9), is incremented with the pressure variation equations and friction force model such that the combined system of equations can be written as where M and Q are the respective accumulated mass matrix and accumulated external force vector as shown in Eq. (9), z ∈ R N b is the vector of the relative joint coordinates,ṗ is the time derivative of the vector of hydraulic pressures, p, and u 0 (p, z,ż) are the pressure variation equations. Here, it is assumed that the dependency of Q with respect to z,ż, p, and F μ , and the dependency of function u 0 with respect to z,ż, and p are known. In static friction models, such as the Bengisu-Akay and Brown-McPhee models, the inclusion of friction force, F μ , in Eq. (35) is a straightforward task, where F μ is dependent on z andż. However, in dynamic friction models, such as the LuGre and modified LuGre models, the inclusion of F μ in Eq. (35) utilizes extra state variables to describe the friction, such aṡ where u 1 (q, z,ż) and u 2 (q, z,ż, h) are the equations of average bristle deflection variation in the LuGre and modified LuGre friction models, respectively, and u 3 (h, z,ż) are the equations of lubricant film thickness variation. The state variables in Eqs. (36) and (37) can be integrated over time in two ways. In the first approach, which is employed in this study, the frictional state variables can be added to the state variables of the coupled equations of motion in Eq. (35) and integrated as shown in this section. In the second approach, the frictional state variables can be integrated locally during resolution of the friction problem [38]. The use of dynamic friction models can introduce numerical stiffness to the system, which may decrease the numerical efficiency of the simulation. In this study, coupled systems are integrated [55] using an implicit single-step trapezoidal rule, as in [7]. It has already been established in [7,40] that the trapezoidal rule performs satisfactorily for coupled multibody and hydraulic systems by considering position and pressure as the primary variables instead of accelerations and pressure derivatives. Similarly, average bristle deflection for the LuGre friction model and average bristle deflection and lubricant film thickness for the modified LuGre friction model are considered as primary variables when applying the trapezoidal rule. The implementation details of the trapezoidal rule can be followed from Appendix A. By following Appendix A, the residual vector, f (x) , can be written as

A case study of a hydraulically actuated four-bar mechanism
In this study, a hydraulically actuated four-bar mechanism, as shown in Fig. 6, is used to illustrate the effect of friction models of a hydraulic cylinder on the four-bar mechanism. The four-bar mechanism is modeled using a semi-recursive formulation, as explained in Sect. 2.1.
As can be seen in Fig. 6, three relative joint coordinates, z 1 , z 2 and z 3 , are used in modeling of the planar structure, and two loop-closure constraints are used for the cut-joint (revolute joint) at point E. The mechanism has one degree of freedom. Note that z 1 , z 2 and z 3 represent, respectively, the positions of the crank, coupler, and rocker in the inertial reference frame, XY Z. The four-bar mechanism, as shown in Fig. 6, is hydraulically actuated using a hydraulic circuit that consists of a pump, a tank, a directional control valve, a throttle valve, a double-acting hydraulic cylinder, and connecting hoses. For simplicity, leakage in the hydraulic components is neglected. Following the lumped fluid theory, the hydraulic circuit is partitioned into discrete volumes (shown in Fig. 6), V 1 , V 2 , and V 3 , whose respective pressures p 1 , p 2 , and p 3 are computed aṡ where Q d1 and Q 3d are the volume flow rates from the directional control valve (Eq. 17), Q t is the volume flow rate from the throttle valve (Eq. 15), B e1 , B e2 , and B e3 are the effective bulk modulus of the respective control volumes (Eq. 14), A 2 and A 3 are the areas of the piston and piston-rod sides, respectively, andṡ is the actuator velocity. For the presented hydraulic circuit, the control volumes can be written as where V h 1 , V h 2 , and V h 3 are the volumes of the respective hoses, and l 2 and l 3 are the chamber lengths of the piston and piston-rod sides, respectively, which can be calculated as l 2 = (l 2 ) t=0 + | s | t=0 − | s | and where | s | is the actuator length of the hydraulic cylinder. The values of s andṡ are, respectively, calculated as whereṙ F is the velocity vector of point F calculated using classical kinematic relations, as in [19,30]. The force, F c , produced by the hydraulic cylinder (Eq. 20) is expressed as where s X , s Y , and s Z are the components of vector s in the inertial reference frame, XY Z. From the static configuration at t = 0, (F c ) t=0 = 681.23 N, ( p 1 ) t=0 = ( p 2 ) t=0 , and where the friction force is assumed as zero and ( p 3 ) t=0 as 3.5 MPa. The system is hydraulically actuated for 10 s using the reference voltage signal, U ref , as where t is the simulation run time. The parameters of the hydraulic circuit are shown in Table 1.
In this study, the coupled system is integrated by incorporating the friction models of the hydraulic cylinder, as shown in Sect. 4. The parameters of the friction models are shown in Table 2, where the values are selected based on typical values found in the literature. Static friction models, that is, the Bengisu-Akay    x = z T , p T , q T T (LuGre approach) In the trapezoidal integration scheme, the error tolerance for the position is 1 × 10 −7 rad; for the average bristle deflection and lubricant film thickness, it is 1 × 10 −7 m; and for the pressure 1 × 10 −2 Pa. The penalty factor, α, in Eq. (9) is considered as 1 × 10 11 in this study. The voltage, which represents the relative spool position, is integrated using the trapezoidal rule, and the error tolerance is considered to be 1 × 10 −7 V. All the approaches in this study are implemented in the Matlab environment.

Results and discussion
This section presents simulation results of the hydraulically actuated four-bar mechanism, shown in Fig. 6, and demonstrates the effect of friction models of the hydraulic cylinder on the four-bar mechanism. The simulation frames of the four-bar structure at different instants of time, representing the positions of the bodies, are shown in Fig. 7. Here, the four coupled system approaches, namely, the Bengisu-Akay approach, Brown-McPhee approach, LuGre approach, and modified LuGre approach, are compared based on the simulation work cycle, friction force, energy balance, and numerical efficiency. In this study, all the simulations are run with a time-step of 1 ms.

The work cycle representation
In the presented case study, the hydraulic cylinder tilts the four-bar structure outward between 1 and 2.5 s and inward between 5 and 8 s. In the subsequent plots, the regions between the opening and closing of the directional control valve are highlighted in purple. For demonstration, only the crank angle is presented in Fig. 8a. However, the coupler and rocker angles follow a similar pattern. Similarly, only pressure p 3 is presented in Fig. 8b; pressures p 1 and p 2 follow a similar pattern. Even though the differences in the relative joint coordinates and pressures may be small, they differ based on the friction model utilized in the model- ing. Such small differences can be an important aspect to consider in applications where high precision and/or low vibration is a pre-requisite. Thus, the friction model can significantly affect the control system performance.
It should be noted that the semi-recursive multibody formulation utilizes the full set of relative joint coordinates in the equations of motion.

Friction force
A comparison of the friction force models of the hydraulic cylinder is shown in Fig. 9. The friction forces in all four approaches are continuous at zero relative piston velocity, thus, avoiding numerical instabilities. Figure 9b compares Fig. 9b, the latter two approaches can incorporate pre-slip displacement and frictional lag. The pre-slip displacement implies a small motion of the bristles in the elastic range, when the applied force is less than the break-away force [38]. It can help to determine the static friction as well as the Stribeck effect in the steady state. Frictional lag is a delayed behavior of the friction force with respect to velocity such that the friction force is higher for acceleration than during deceleration.

Energy balance
The energy balance in the approaches is analyzed by comparing the potential energy, kinetic energy, and work done by the hydraulic cylinder. The energy comparison is shown in Fig. 10, where the energy balance in all four studied approaches showed good agreement with respect to each other. With the Bengisu-Akay approach, the peak energy drift is 0.45 J, which is 0.37% of the maximum cylinder work, 121.02 J, whereas with the Brown-McPhee approach, the peak energy drift is 0.45 J, which is 0.37% of the maximum cylinder work, 120.99 J. With the LuGre approach, the peak energy drift is 0.44 J, which is 0.36% of the maximum cylinder work, 120.94 J, and in the modified LuGre approach, the peak energy drift is 0.44 J, which is 0.36% of the maximum cylinder work, 120.97 J. It should be noted that in all four approaches, the peak energy drift occurred at the opening of the directional control valve at around 1 s. An analogy can be drawn between the hydraulic system and a stiff spring supporting the fourbar structure.

Numerical efficiency
The numerical efficiencies of all four approaches are compared in Figs. 11 and 12, and the average and maximum iterations and total integration time are shown in Table 3. The maximum number of iterations occurred during the opening and closing of the directional control valve. Otherwise, the iterations are mostly one iteration in the Bengisu-Akay and Brown-McPhee approaches and two iterations in the LuGre and modified LuGre approaches. During closing of the directional control valve, the LuGre and modified LuGre approaches required more iterations than the Bengisu-Akay and Brown-McPhee approaches (Fig. 11). The integration time of the LuGre and modified LuGre approaches is greater than that of the Bengisu-Akay and Brown-McPhee approaches, which is in agreement with the literature [37,38]. The poor numerical efficiency of the LuGre and modified LuGre approaches can be attributed to the numerical stiffness introduced by use of extra state variable(s), such as average bristle deflection in the LuGre approach and  Brown-McPhee approach is relatively the most efficient approach and is a potential candidate for real-time applications [16]. Real-time capabilities can be achieved by implementing the approach in a lower language, such as C++ or Fortran, and by limiting the number of iterations to a lower value.

Conclusion
In simulation of coupled multibody and hydraulic systems, an important aspect in the modeling is the nonlinearity and numerical stiffness introduced by the friction model of a hydraulic cylinder. The objective of this paper was to couple various friction models of a hydraulic cylinder with the equations of motion of a hydraulically actuated multibody system in a monolithic framework. The study considered two static friction models: the Bengisu-Akay model and the Brown-McPhee model, and two dynamic friction models: the LuGre model and the modified LuGre model. Within the monolithic framework, the multibody system was described using a semi-recursive formulation and the hydraulic actuator was described using the lumped fluid theory. A hydraulically actuated four-bar mechanism was examined as a case study and the four approaches were compared based on the work cycle, friction force, energy balance, and numerical efficiency. The relative joint coordinates and hydraulic pressures were affected by the friction force model selected. This is a significant aspect that should be considered in modeling, especially, in applications where high precision and/or low vibration is a pre-requisite. The Bengisu-Akay and Brown-McPhee approaches showed traditional friction-velocity characteristics, namely the Coulomb, stiction, Stribeck, and viscous friction effects. The LuGre and modified LuGre approaches demonstrated frictional phenomena such as pre-slip displacement and frictional lag. All four approaches studied showed acceptable results for energy balance. The modified LuGre approach had the poorest numerical efficiency, which was attributed to the solution of additional state variables: average bristle deflection and fluid film thickness. The Bengisu-Akay approach was more efficient than the LuGre approach but less efficient than the Brown-McPhee approach because of the change in friction equation during the simulation. It can be concluded that the Brown-McPhee approach appears to be numerically the most efficient approach and it was well able to describe usual friction characteristics in dynamic simulation of hydraulically actuated multibody systems.
In future work, the Brown-McPhee approach could be implemented in a lower language, such as C++ or Fortran, to assess real-time capabilities further. The methods presented in this study may be extended to study of state and force estimators. Studies may also be directed toward coupling the dynamic friction models with other multibody formulations. Future studies can investigate the comparative responses of multibody systems based on the high and low magnitudes of the stiffness/damping of dynamic friction models of a hydraulic cylinder and based on throttling in pressure lines. This can lead to a possible source for a stick-slip behavior between the hydraulic cylinder and piston. Furthermore, static friction models that introduce discontinuity at zero relative velocity can be considered in future studies to demonstrate the numerical instability caused by them.
Service Business from Physics Based Digital Twins-DigiBuzz) and in part by the Academy of Finland under Grant #316106.

Conflicts of interest
The authors declare that they have no conflict of interest.
Availability of data and material Not applicable.
Code availability Not applicable.
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/.

∂f (x) ∂x
where f (x) is the residual vector and ∂f(x) ∂x is the approximated tangent matrix.