A hydraulic cylinder model for multibody simulations

Computers and Structures 138 (2014) 62–72

Contents lists available at ScienceDirect

Computers and Structures journal homepage: www.elsevier.com/locate/compstruc

A hydraulic cylinder model for multibody simulations Antti Ylinen a,⇑, Heikki Marjamäki a, Jari Mäkinen b a b

Tampere University of Technology, Department of Mechanical Engineering and Industrial Systems, P.O. Box 589, FIN-33101 Tampere, Finland Tampere University of Technology, Department of Civil Engineering, P.O. Box 600, FIN-33101 Tampere, Finland

a r t i c l e

i n f o

Article history: Received 17 April 2013 Accepted 26 February 2014 Available online 27 March 2014 Keywords: Hydraulic cylinder Coupled problem Multibody dynamics Friction

a b s t r a c t In this paper we present a model for a linear hydraulic actuator for multibody simulations. In addition to translational degrees of freedom also the cylinder chamber pressures are taken as state variables. In static analysis the chamber pressures are embedded thus a purely mechanical simulation can be performed. In dynamic analysis we end up to a coupled problem where the chamber pressures are separate variables. The sealing friction is taken into account in both static and dynamic analysis. For numerical solution we use a monolithic algorithm and the coupling matrices between the hydraulic and mechanic variables are presented. Ó 2014 Elsevier Ltd. All rights reserved.

1. Introduction Hydraulic driven working machines are commonplace in the industry of today. The applications can vary from very robust excavators to very precise manufacturing robots but in outline the hydraulic actuators are quite similar in mechanical sense. Simulation of these applications can derive from need of fatigue assessment or simply to determine the peak stresses due to motion. Engineers are also interested in the positioning accuracy of the systems. In this paper we will propose a new hydraulic cylinder model for dynamic simulations of hydraulic driven multibody systems. The proposed cylinder element couples the hydraulic and mechanical variables in a monolithic way. We introduce the tangent matrices concerning the mechanical system and the hydraulic system but the analytical coupling matrices as well. In addition the cylinder element is implemented with a novel friction model. However, the dynamic friction model is not suitable for static computations. Therefore we derive a static friction model based on the parameters of the dynamic friction for defining the initial state of the coupled system. Hydraulic cylinder produces linear movement for the multibody system, for instance a lifting boom of a personnel carriage. The movement is then controlled by the flow rate lead to the cylinder chambers. However the linear movement for the simulation models is traditionally dealt with constraint equations. The classical formulation for constraint equations can be found in [1] but ⇑ Corresponding author. Tel.: +358 401981631. E-mail addresses: antti.ylinen@tut.fi (A. Ylinen), (H. Marjamäki), jari.m.makinen@tut.fi (J. Mäkinen). http://dx.doi.org/10.1016/j.compstruc.2014.02.006 0045-7949/Ó 2014 Elsevier Ltd. All rights reserved.

heikki.marjamaki@tut.fi

more modern finite element approach is presented in [2]. Solving systems with constraint equations can be done by utilizing for example the Lagrange multiplier method or penalty method. However, as shown in [3] this method can introduce spurious oscillation to the response of the system. More sophisticated way to model the hydraulic cylinder is a length controlled rod element where the unstressed length of the element is given as a function of time [4]. The compressibility of the hydraulic fluid can is then taken into account as a reduced Young’s modulus of the material attached to the element. However, the length change has to be a predefined function estimated from the stationary state of the hydraulic control system in the actual appliance which is simulated. In addition the length controlled rod element has not an embedded friction model. The friction plays an important role, especially when the extension rates are low and there is a possibility for a stick–slip effect. The cylinder element presented in this paper changes length as a function of the inbound flow rate. The flow rate induces pressure to the cylinder chamber and this force extends the cylinder. Noting that the cylinder force is a function of the chamber pressures and cylinder displacements we find that the cylinder element is a coupling element. Friction plays an important role in dynamical analyses by introducing damping to the system. The very basic friction model is the Coulomb friction where the friction force is dependent only on the contact force [5]. However, in systems as the hydraulic cylinder, the friction force is known to be a function of the sliding velocity thus a dynamical friction model is required to capture the velocity dependent properties [6]. Furthermore the stick–slip phenomenon where the sliding is not continuous is an important factor to take into account. In stick–slip phenomenon the friction alternates

63

A. Ylinen et al. / Computers and Structures 138 (2014) 62–72

between the static friction and sliding friction. The stick–slip phenomenon and the reduction in friction force from static to sliding can be modeled using approach in [7]. Through the stick–slip effect the extension of the cylinder is not continuous and it can induce high stress peaks to the system. The dynamical friction model however is not suitable for initial state computation because of the velocity dependency of the model. Moreover, the dynamical friction model introduces a new variable to the system. Therefore we derive a friction model for the initial state computations based on the parameters of the dynamical friction model. In this paper the numerical treatment of the coupled system is treated in a monolithic ways, as it was first introduced in [8]. For time integration multi rate integration is exploited where the time stepping scheme for hydraulic cylinder is different from the one of the mechanical system. Information is however changed during each time step. Also in [9] the system was solved monolithically and different integration schemes were compared however. In this paper we derive a hydraulic cylinder model for multibody simulations and the paper is organized as follows. First we discuss the equations of motion for the coupled system. From linearized equations of motion for the coupled system we present the tangential matrices for the system and in Section 3 we present these matrices. In Section 3 we also present the cylinder element for initial state computations. In Section 4 we discuss the time integration shortly and in the final section we present three different numerical examples. In the first example we deal with the initial state of the cylinder element with constant load. The second example follows from the first example where we present a simple dynamic simulation. In the final example we study a lifting boom and draw comparison between different linear actuator modeling methods.

where we identify the tangent matrices for the mechanical system as well as for the hydraulic system. The subscript m refers to the mechanical system and c to the cylinder variables. The notation is as follows: The first subscript refers to the state equation to be differentiated and the latter to the variable with respect the differentiation is taken to. The mass matrix is the derivative of the kinetic energy

M¼

@2T _ q_ @ q@

The total mass of the system changes due to the changes in cylinder chamber volumes thus the center of the mass also varies. The damping matrix of the mechanical system is

Cmm ¼ 

Kmm ¼ 

t2

_ dLðq; qÞdt ¼

t1

Z

½dT  dVdt ¼ 0

ð1Þ

t1

€ ¼ gðq; q; _ z; tÞ Mq _ z_ ¼ f cyl ðz; q; qÞ

ð2Þ

where M is the mass matrix and the vector g is the sum of external, internal and complementary inertial forces. Vector z collects the state variables of the hydraulic cylinder. Coupling between the mechanical system and the hydraulic cylinder is visible in the right hand side terms of Eq. (2) where we note the crosswise dependency of the systems. For the solution of the dynamical system we linearize the state equations yielding a linear set of equations with monolithic coupling



M 0 0

0



€Þ @g @ðMq þ @q @q

ð6Þ

Kcc ¼ 

@f cyl @z

ð7Þ

The off-diagonal matrices in the linearized equations of motion in (3) are the coupling matrices. They rise from the crosswise derivation of the state equations thus yielding in non-square matrices. The coupling matrices are then written as

@f cyl @ q_

Kcm ¼ 

@f cyl @q

Kmc ¼ 

@g @z

ð8Þ

_ z; tÞ  Mq € and s ¼ zþ _ The right hand side terms r ¼ gðq; q; _ f cyl ðz; q; qÞ in Eq. (3) are the residuals of the mechanical system and hydraulic cylinder equilibriums respectively. We also note that the coupling damping matrix Ccm is an actual damping matrix yielding from the friction model. 3. Cylinder element

t2

where T is the kinetic energy of the system and V is the potential energy. The potential energy consists of the potential of the internal forces and external forces. Vector q is the generalized coordinates describing the system. From the Hamilton’s principle we can write the equations of motion for the mechanical system [10] together with the state equation of the hydraulic cylinder

(

ð5Þ

where we note that the latter matrix rises from the inertial forces. For the cylinder element we obtain the hydraulic tangent as

2. Equations of motion

Z

@g @ q_

The tangential stiffness is then defined as

Ccm ¼ 

In this section we present shortly the equations of motion for the coupled hydromechanical simulation model as well as its linearization. The time integration procedure is then discussed in Section 4. The Hamilton’s principle when only conservative external forces are applied to the system states

ð4Þ

  € Cmm Dq þ Ccm D€z

0 I



  Kmm Dq_ þ Kcm Dz_

Kmc Kcc



   Dq r ¼  Dz s ð3Þ

In this section we present the hydraulic cylinder element and the state Eq. (2) in detail. First we derive the hydraulic cylinder with the dynamic friction model and treat with the coupling between the cylinder and the mechanical system. In the final section we present the cylinder model for initial state computations. The hydraulic cylinder element is a rod like element with two nodes and three translational degrees of freedom in each node, see Fig. 1. In addition, the cylinder element introduces extra variables: The chamber pressures and a friction induced variable discussed later. Cylinder piston position is xc ¼ Ln  L0 where L0 is the initial length of the cylinder and Ln is the current length. The force that the cylinder produces can be written in terms of the chamber pressures and the friction force as

F c ¼ pA AA  pB AB  ffr

ð9Þ

where pA and pB are the pressures of the plus and minus chamber, respectively. The corresponding areas of the cylinder piston are AA and AB . Friction force is denoted ffr . The unit vector in direction of the axis of the cylinder element is given by using the current nodal coordinates as

nc ¼

xB  xA xB  xA ¼ kxB  xA k Ln

Now the internal force vector can be expressed as

ð10Þ

64

A. Ylinen et al. / Computers and Structures 138 (2014) 62–72

where kpr is the pressure coefficient. For the velocity dependent term we the use model introduced in [6]. The friction force is given

fl ¼ k0 z þ k1 z_ þ F v x_ c

ð15Þ

where coefficients k0 ; k1 and F v are stiffness coefficient, damping coefficient and viscous friction, respectively, and z is the bristle deformation, see Fig. 2. The first time derivative of the bristle deformation is defined as [6]

dz jx_ c j ¼ x_ c  z dt gðx_ c Þ

ð16Þ

The function gðx_ c Þ in the denominator is

gðx_ c Þ ¼

Fig. 1. Hydraulic cylinder element with main dimensions and variables.

 f int ¼ F c

nc nc

 ¼

  F c Cx Fc ¼  Ax Ln Cx Ln

ð11Þ

where the vector x ¼ ðxA ; yA ; zA ; xB ; yB ; zB ÞT and matrices A and C are

2

3 1 0 0 1 0 0 6 7 C ¼ 4 0 1 0 0 1 0 5 0 0 1 0 0 1

ð12Þ

 2 1 F Cou þ ðF st  F Cou Þeðx_ c =v Str Þ k0

ð17Þ

where F Cou and F st are the Coulomb the and static friction, respectively. The Stribeck velocity v Str attributes from the Stribeck effect where the friction force is function of the sliding velocity. In this paper we give the static and the Coulomb friction forces as constant values instead of computing them from the contact normal force. Using the approach described in [7] the tangential force can be computed from the contact normal force and it explains also the stick– slip phenomenon. The reduction from static friction to Coulombic friction is through high frequency vibration, in direction of contact normal, when relative movement occurs. Since the contact normal force is hard to define in cylinder applications, the model proposed in [6] seems to be more suitable for present simulations. Moreover, the tangential quantities can be measured. This model is also suitable for the cylinder application where contact can never break and it has been a starting point for this model. In more general contact modeling the approach in [7] is more suitable. 3.2. Cylinder variables

A ¼ CT C The p current ffiffiffiffiffiffiffiffiffiffiffiffi length for the element can be written also as Ln ¼ xT A x which is used in subsequent derivations. 3.1. Friction model We write the friction force as a product of two functions where one depends on the chamber pressures and the other on the velocity of the piston x_ c and the bristle deformation z, see Fig. 2. The bristle can be interpreted as the piston sealing ring deformation. The chosen friction model takes the most important factors such as Stribeck effect, hysteresis and stick–slip phenomenon into account [5,11]. The total friction force is a product of two functions

ffr ¼ fpr ðpA ; pB Þfl ðz; x_ c Þ

ð13Þ

Friction force grows as a function of the chamber pressures thus the linear dimensionless scaling factor is

fpr ðpA ; pB Þ ¼ 1 þ kpr ðpA þ pB Þ

ð14Þ

The state equation for the cylinder element can be written as

_ z_ ¼ f cyl ðz; q; qÞ

ð18Þ

where the vector q and q_ are the nodal displacements and velocities of the mechanical system respectively. The state equation in component form is

3 Q A  x_ c AA B 7 2 3 6 V þx A c A 7 6 A p_ A 7 6 _ 6 _ 7 6 Q B þ xc AB 7 4 pB 5 ¼ 6 B7 6 V B  xc A B 7 7 6 z_ 5 4 jx_ j x_ c  _c z gðxc Þ 2

ð19Þ

where B is the bulk modulus of the hydraulic fluid and V A and V B are initial fluid volumes in cylinder chambers. The bulk modulus is a function of the fluid pressure but in practical engineering applications it can be assumed as constant [8] thus no differentiation with respect to the bulk modulus is needed in following sections. Flow rates in and out of the cylinder are Q A and Q B and the positive flows are shown in Fig. 1. The coupling with the mechanical system is through the piston displacement xc and velocity x_ c . 3.3. Tangent matrices In this section we derive the tangential and coupling matrices for the dynamical simulation model. For the tangential matrices concerning the mechanical system we take the first variation of the internal force vector in Eq. (11) with respect to the mechanical variables

df int ¼  Fig. 2. Piston sealing ring deformation z.

  Fc 1 1  AxdF c Adx  F c Axd Ln Ln Ln

ð20Þ

65

A. Ylinen et al. / Computers and Structures 138 (2014) 62–72

Straightforwardly we compute the variation of the inverse of the current cylinder length as

  1 ¼ d Ln

dLn L2n

! ¼

1 L3n

xT Adx

dz

ð23Þ

where the function Hðz; x_ c Þ is

z g 0 ðx_ c Þ sgnðx_ c Þ þ 2 jx_ c jz gðx_ c Þ g ðx_ c Þ

2 2x_ c g 0 ðx_ c Þ ¼  2 ðF  F Þeðx_ c =v Str Þ v Str k0 st Cou

1 T x Ax_ Ln

dx_ c ¼

ð24Þ

And now the stiffness matrices and damping matrix for the cylinder element can be written

! fpr cl Fc Fc A þ 3 AxxT A  AxB dx Ln Ln Ln ! fpr cl _ 1 þ k2 þ k3 Þdx þ Cmm dx_ þ  2 AxxT A dxðk Ln 

ð26Þ

From Eqs. (25) and (26) we note that the damping matrix Cmm and the geometric stiffness matrices k1 and k2 are symmetric whereas the material stiffness k3 is unsymmetric. The unsymmetricity rises from the variation of the cylinder piston velocity x_ c in Eq. (25). The magnitude of the unsymmetric part, k3 , can be estimated with an expression fpr cl v =Ln where v is a typical velocity for the cylinder. For the symmetric parts, k1 and k2 , we find the estimated magnitude to be F c =Ln and the quotient of these two is fpr cl v =F c  1. With the values from the numerical example 3 we find the quotient to be between 0.01% and 0.1% thus the asymmetry can be neglected. When the typical velocity of the cylinder is higher like in helicopter blades, the unsymmetric part improves convergence. The Jacobian matrix of the cylinder element can be computed as a variation of the cylinder state equation (18) to the cylinder variables yielding

0 0

0

The row numbers in Kmc correspond to the nodal displacement vector dx and the columns to the variations of the hydraulic cylinder variables ½ dpA dpB dz T . The second coupling matrix Kcm rises from the variation of the hydraulic state equation with respect to the mechanical variables:

6 df cyl ¼ 6 4

_

xc AA A B  ðVQ Aþx A Þ2 A A

c A

Q B þx_ c AB ðV B xc AB Þ2

3

2

7 6 6 AB B 7 5dxc þ 4

AA  V A þx B c AA

3

7 B 7 5dx_ c ¼ bcm dxc þ ccm dx_ c _ Hðz; xc Þ

AB V B xc AB

ð30Þ Substituting the variations of the cylinder piston displacement and velocity we obtain the coupling matrices as

df cyl ¼

1 T 1 1 1 _ T A dx þ xT Adx_ ¼ Bdx þ xT Adx_ x_ A  3 xT Axx Ln Ln Ln Ln

2

ð29Þ

  1 1 ccm B þ bcm xT A dx þ ccm xT Adx_ Ln Ln ¼ Kcm dx þ Ccm dx_

!

¼ Kmm dx þ Cmm dx_

1 Axbmc dz ¼ Kmc dz Ln

0

ð25Þ

df int ¼

df int ¼ 

2

Variation of the cylinder piston velocity can then be derived by writing the piston velocity as

x_ c ¼

ð28Þ The coupling matrix can be written by taking account the last term in Eq. (20) and substituting into Eq. (28)

noting that the pressure dependent variation dfpr with respect to the mechanical variables is zero. Variation of velocity term of the friction force is given as

Hðz; x_ c Þ ¼ 1 

h   i dpA _ 6 7 dF c ¼ AA  kpr fl AB  kpr fl fpr k0  k1 gðjxx_ccjÞ 4 dpB 5 ¼ bmc dz

ð22Þ

d f l ¼ ðF v þ k1 Hðz; x_ c ÞÞdx_ c ¼ cl dx_ c

3

ð21Þ

The variation of the cylinder force in Eq. (9) is

d F c ¼ fl d f pr  fpr d f l

2

The rows of these coupling matrices correspond to the cylinder variables and the columns to the nodal displacements and velocities. For the dynamical simulation we derive also the mass matrix of the cylinder element. For this we divide the mass of the cylinder into four parts being the cylinder lining G1 , cylinder rod G2 and hydraulic fluids in plus and minus chambers, G3 and G4 respectively, as Fig. 3 illustrates. Each of these parts are then considered as point masses. To derive the mass matrix we identify the locations for these masses and using the dimensions given in the Fig. 3 we can write the dimensionless positions measured from the bottom of the cylinder

3

6 0 7 df cyl ¼ 4 0 0 5dz ¼ Kcc dz _ 0 0  gðjxx_ccjÞ

ð27Þ

The interaction or coupling matrices are derived as crosswise variations of the internal force and cylinder state equation leading to rectangular matrices. First we derive the upper coupling matrix Kmc where we take the variation of the cylinder force with respect to the cylinder variables

ð31Þ

Fig. 3. Point masses of the cylinder element.

66

A. Ylinen et al. / Computers and Structures 138 (2014) 62–72

g p Lp

Lap þ xc þ gr Lr ; n2 ¼ Ln Ln 2Lap þ xc Lr þ Lap þ xc n3 ¼ ; n4 ¼ 2Ln 2Ln n1 ¼

dm3 ¼ qAA dxc ¼ ð32Þ dm4 ¼ 

where the numbering corresponds to the mass parts in order of cylinder lining, cylinder rod, fluid in plus chamber and finally in minus chamber. The dimensions are shown in Fig. 3 and the parameters gp and gr are to define the distance of center of gravity of the lining and rod. For interpolation we use linear shape functions and they are collected to a matrix as

N i1 ¼ 1  ni ; Ni2 ¼ ni 2 Ni 0 0 Ni2 6 1 i Ni ¼ 6 0 4 0 N1 0 0 0 N i1 0

0 Ni2 0

0

¼

ð33Þ

7 0 7 5 Ni2

4 4 1X 1X _ T Ni q_ mi u_ Ti u_ i ¼ mi ðNi qÞ 2 i¼1 2 i¼1 4 1 TX 1 mi NTi Ni q_ ¼ q_ T Mq_ q_ 2 i¼1 2

ð34Þ

where the mass mi is constant for the cylinder lining and rod. However, the mass of the hydraulic fluid is

m3 ¼ qðV A þ xc AA Þ m4 ¼ qðV B  xc AB Þ

ð35Þ

Because the mass matrix is a function of the cylinder piston displacement it induces gyroscopic stiffness to the system. This stiffness matrix is then derived by taking the first variation of the inertial forces 4 X € df iner ¼ d mi NTi Ni q

ð36Þ

i¼3

  i¼1  € ¼ d NTi Ni q € where q € remains constant. After noting that d NTi Ni q some algebra we obtain the variation as

!   ci di T € ¼ ð2ni A þ DÞq € dni ¼ ð2ni A þ DÞq € dLn d Ni Ni q þ Ln L2n ! 1 ci di € xT Adx ¼ þ 2 ð2ni A þ DÞq Ln Ln Ln

ð37Þ

where the coefficients ci and di are shown in Table 1. The coefficients ci and di are derivatives of Eq. (32) to xc and Ln respectively. The auxiliary matrix D is

D¼

2I

I

I

0

xT Adx ð39Þ

xT Adx

Substituting these equations to Eq. (36) we obtain the gyroscopic stiffness matrix. We define the external load due to the gravitation as

fG ¼

4 X

NTi Gi

ð40Þ

i¼1

Gi ¼ mi g x

gy

gz

T

¼ mi g

ð41Þ

Now we can derive the stiffness matrices due to the gravitational load by taking the first variation of the external load vector

df G ¼

" # 4 X Ni1 gdmi i¼3

Ni2 gdmi

þ

" # 4 X mi Ni1 gdNi1 i¼1

mi Ni2 gdNi2

ð42Þ

where the mass variations are given in (39). The latter part is then varied as follows because dNi1 ¼ dni and dN i2 ¼ dni " # " # !" # 4 4 4 X X X mi Ni1 gdNi1 mi Ni1 g mi Ni1 g T 1 c i di ¼ dni ¼ x Adx þ i i L Ln L2n mi Ni2 g mi Ni2 g i¼1 mi N 2 gdN 2 i¼1 i¼1 n ð43Þ

Substituting these to (42) we get the external load stiffness matrix. In traditional engineering applications such as in lifting booms or excavators the gyroscopic tangent and the stiffness of external loads can be neglected since their influence is small in comparison with the geometric and material stiffness. However in computational cases where the angular velocity of the cylinder element is significant the gyroscopic tangent matrix improves the convergence of the solution algorithm. Such applications are for instance the servo cylinders of the helicopter blade adjustment cylinders. 3.4. Initial state

i¼1 4 4   X X € dmi þ € ¼ mi NTi Ni q mi d NTi Ni q



Ln

Ln

The vector Gi can be expressed using the vector of gravitational accelerations and mass as

3

where the index i refers to the mass particles. We derive the mass matrix from the kinetic energy of the point masses and using the linear interpolation for the velocities of the point masses u_ i

T ¼

qAB

qAA

 ð38Þ

The mass variations in Eq. (36) are zero for the cylinder lining and rod but for the hydraulic fluid volumes we get Table 1 The coefficients of the ni derivatives. i

ci

di

1

0

g p L p

2

1

ðLap þ xc þ gp Lp Þ

3 4

1/2 1/2

ðLap þ 1=2xc Þ ðLap þ xc þ Lp Þ=2

For the dynamical simulation we need to compute the initial state of the system. Because the chosen friction model is purely dynamical friction model we also need a static model for the friction. In this section we derive a cylinder element for initial state computations as well as the static friction model. We also define a procedure for changing the static friction to dynamic friction. First we start by deriving the static cylinder model without friction by taking the definition of the bulk modulus for the hydraulic fluid as

Dp DV ¼ B V0

ð44Þ

where V 0 is the initial volume and DV and Dp are change in volume and pressure respectively. The bulk modulus is assumed as constant. From this we can write the pressure change in cylinder chambers

Dxc A A B V A þ Dxc AA Dxc AB DpB ¼ B V B  Dxc A B DpA ¼ 

ð45Þ

where the change in cylinder piston displacement is Dxc ¼ Ln  L0 . When the initial state or the unstressed length is known we can compute the absolute pressures in cylinder chambers using Eq. (45) thus we can write the cylinder force of the frictionless cylinder as

A. Ylinen et al. / Computers and Structures 138 (2014) 62–72

F c ¼ DpA AA  DpB AB

ð46Þ

Using same methods as presented, see Eq. (20), we can write the tangential stiffness matrix as a first variation of the internal force vector df int ¼ 

! ! Fc Fc B A2A V A A2B V B T Axx dx A þ 3 AxxT A þ 2 þ A Ln Ln Ln ðV A þ AA Dxc Þ2 ðV B  AB Dxc Þ2

¼ ðkg1 þ kg2 þ kmat Þdx ¼ Kdx ð47Þ

Unlike in the dynamical cylinder model, in the static model the material stiffness matrix is also symmetric. When adding the friction model to the static cylinder we take it into account by subtracting the friction force from Eq. (46) as in Eq. (9). The friction force introduced in Section 3.1 is purely dynamical so we derive a new model for the initial state computation. The static friction is derived as a function of the piston displacement fls ¼ fls ðDxc Þ. For the total friction we use same pressure dependency as for the dynamical model in Eq. (13) giving the static friction as

ffrs ¼ fprs ðDpA ; DpB Þfls ðDxc Þ

ð48Þ

The schematic friction curve is shown in Fig. 4 where we have marked the points A–C. The point A corresponds to the initial state with zero friction and no displacements and we define the slope to be the stiffness coefficient k0 from the dynamical friction. At point B we define the static friction with given displacement and at point C the Coulomb friction. For the friction curve we have fitted 5th order polynomial through points A–C to get almost linear response with small displacements. However the 5th order polynomial behaves badly on the softening part of the curve from B to C thus we have replaced it with Hermite polynomial. From point C on we have constant Coulomb friction. The friction curve is C1 differentiable in its domain. The parameters are taken from the dynamic friction model to get smooth transition from the initial state computation to the dynamical simulation. The switching between the models is done when the Newton–Raphson algorithm has converged to a solution. While the dynamic friction model introduces a new variable, piston sealing ring deformation z, we need a special method for computing the initial state for this variable. The friction force is a function of the piston displacement xc thus a stiffness matrix rises from the variation.

dffrs

¼

s fpr

! @fprs 1 dxc ¼ csl xT Adx þ fl Ln @xc @xc @fls

s

ð49Þ

where fpr is given in Eq. (14). And finally the material stiffness matrix of the friction is

1 s c AxxT A Ln l

ð50Þ

1 Friction force fµ/N

kl ¼

C

flini ¼ fl ðz; x_ c Þ

where flini is the displacement term of the static friction model and the right hand side corresponds to the dynamic friction model. The velocity dependent part in Eq. (15) can be simplified because cylinder velocities are zero. From these conditions the initial bristle deformation is

zini ¼

flini k0

0.2 A 0.1

ð52Þ

Using the methods discussed in this section it is possible to compute initial state for a system with a hydraulic cylinder element. In static computations we note that in comparison with the dynamical model the pressures are secondary variables because they can be expressed directly with the cylinder piston displacements. Therefore they can be eliminated from the equations unlike in the dynamical simulation model where the pressures are actual state variables. 4. Time integration In this section we discuss the time integration. The Newmark scheme is a widely used method for solving second order ordinary differential equations but it is not specifically tuned for stiff differential equations. For stiff differential equations one needs typically some energy decaying or conserving algorithms as shown in [12,13]. As shown in [8] the hydraulic system is integrated using the Crank Nicolson scheme and mechanical system with the Newmark method. The multi rate integrator concept has also been discussed in [14,9] and it improves the performance of the time stepping. However, due to simplicity the single step Crank Nicolson scheme, i.e. the trapezoidal rule is implemented for the whole coupled system [15]. First we define the predictors for the new time step which are taken from the converged solution of the previous time step. We write the zero acceleration predictors as

€ 0nþ1 ¼ 0 q 1 €n q_ 0nþ1 ¼ q_ n þ hq 2 1 2 €n q0nþ1 ¼ qn þ hq_ n þ h q 4 z0nþ1 ¼ zn   1 1 z_ 0nþ1 ¼ z0nþ1  zn  hz_ n 2 ch

1 _k Dq q_ kþ1 nþ1 ¼ qnþ1 þ ch   1 1 _ € € kþ1 q_ kþ1 q nþ1 ¼ nþ1  qn  hqn 2 ch

0.4

0.05

ð51Þ

ð53Þ

k qkþ1 nþ1 ¼ qnþ1 þ Dq

0.6

0

The geometric stiffness of the friction force is taken into account in the first two terms in Eq. (47) where the friction force is computed into the cylinder force F c . The initial bristle deformation for the dynamical simulation can be computed as follows: the static friction force and the chamber pressures are known and therefore we can write

where the parameter c ¼ 1=2 for the Crank Nicolson scheme. The corrections for the predictors can then be given as

B

0.8

0

67

0.15

0.2

0.25

Piston displacement xc /mm

Fig. 4. Schematic presentation of the friction curve.

0.3

k zkþ1 nþ1 ¼ znþ1 þ Dz   1 1 _ zkþ1 z_ kþ1 nþ1 ¼ nþ1  zn  hzn 2 ch

ð54Þ

68

A. Ylinen et al. / Computers and Structures 138 (2014) 62–72

where the corrections are computed from the linearized equation of motion



ST

Uc

Lc

Hc



   Dq r ¼  Dz s

ð55Þ

The components of the system matrix can then be given

1 M 2 ðchÞ Lc ¼ chKcm þ Ccm

ST ¼ Kmm þ Uc ¼ Kmc

1

ch

Cmm þ

ð56Þ

Hc ¼ chKcc þ I and the right hand side for the mechanical system is as introduced in Section 2. For the hydraulic cylinder we use as residual [8]

_ s ¼ zkþ1 nþ1 þ zn þ chðzn þ f cyl Þ

ð57Þ

where the cylinder state f cyl is computed as shown in (18) and at iteration step k. 5. Numerical examples In this section we present three different numerical examples to show the performance of the proposed cylinder element. In our first example we study the convergence of the computation in the initial state calculations. Results are presented for the both frictionless as well as for frictional element. In the second example we consider lifting of a point mass initially resting at a static equilibrium state. A constant flow rate is prescribed to the plus chamber of the cylinder for a period of one second, after which the flow rate is immediately stopped. The computation is again carried out in frictionless and frictional case. In our final example we are seeking the peak normal stresses of a lifting boom due to bending. Three different simulation models are used. First Lagrange multiplier method and the length controlled rod element are used to produce the lifting movement. These stress results are then compared to the results obtained by using the proposed cylinder element in frictionless case. Finally the frictionless cylinder element results are compared to the results brought by the frictional element with two different friction coefficients. 5.1. Static analysis of a cylinder In the first example a point load is attached to the rod eye of a vertical cylinder as shown in Fig. 5. The convergence of the computation is studied in a frictionless case as well as in frictional case using three different friction coefficients. The results are computed in both cases cylinder up and upside down. The number of Newton–Raphson iterations and the chamber pressure for the frictionless cylinder and additionally the bristle deflection z and friction force for the frictional cylinders are gathered. In static analysis only the chamber areas and comparable fluid volumes, the oil compressibility and the friction parameters are needed. Data of the cylinder is given in Table 2. In our test case the external load has been chosen to induce chamber pressure of 20 MPa when no friction occurs. The load is given as a point load and gravity loads are not included. The static friction force is then determined as a fraction of the external force. In this example the Coulomb friction is chosen to coincide to the static friction. Thus there is no softening part in the friction force, see Fig. 4. It should be pointed out that at the first iteration step both cylinder chambers carry the load. At the second step only the chamber with positive pressure provides stiffness to the cylinder. This procedure is implemented to figure the direction of the load.

Fig. 5. Load cases for the cylinder statics. On the left cylinder up and on the right cylinder down.

Table 2 Cylinder data for static analysis. Quantity

Symbol

Value

Unit

Chamber area Chamber area Chamber volume Chamber volume Cylinder initial length Oil compressibility Stiffness coefficient

AA AB VA VB L0 B k0

0.0038 0.0023 0.00084 0.0013 1.118 1493 29.86

m2 m2 m3 m3 m MPa MPa

Convergence of the Newton’s iteration is attained when the norm of the residual vector is millionth part of the norm of the external load. This is justified in sense of engineering applications. In this example we have chosen the Coulomb friction equal to the static friction and the stiffness coefficient k0 is relatively high. Thus friction grows from the static region to the sliding friction with small displacement changes. The results of the example are presented in Table 3 where on the upper part we show results for the case 1 (cylinder up) and on the bottom part for the case 2 (cylinder down), see Fig. 5. First to notice is the good convergence of the element within friction coefficient range 0 6 l 6 0:25. It takes only 4 steps to converge to the given tolerance. An extra iteration round is required in first case with friction coefficient l ¼ 0:5 because the friction is not fully developed and iterations are done in the rising part of the

Table 3 Cylinder statics result information.

l

Fl

z (mm)

Nr. steps

pA (MPa)

Cylinder up 0.00 0.10 0.25 0.50

0.0 7696.9 19242.3 38484.5

0.0 0.258 0.644 1.149

4 4 4 5

20 18 15 11.85

l

Fl

z (mm)

Nr. steps

pB (MPa)

0.00 0.10 0.25 0.50

0.0 4516.0 11290.1 22580.2

0.0 0.151 0.378 0.56

4 4 4 5

20 18 15 10

Cylinder down

69

A. Ylinen et al. / Computers and Structures 138 (2014) 62–72

friction curve. In the second case an extra iteration round is needed due to poor guess of the initial stiffness. We can say that the static cylinder element gives realistic results and the convergence of the Newton’s iteration is quadratic. The effects of the sealing friction can be easily taken into account. The static cylinder element is an effective way to compute the initial state of the system and to obtain the chamber pressures as well as the bristle deflection. Surely the initial state could be computed using dynamic analysis with damping and long settling time but running an initial static analysis is computationally much more effective.

5.2. Lifting of a point mass In the second example we consider lifting of a point mass, see Fig. 6, where a mass element is attached to the rod eye. Static analysis is first performed to obtain the initial state of the system. The mass is chosen to induce 20 MPa chamber pressure. A constant flow rate is then prescribed to the cylinder plus chamber for a one second period. After that the flow is immediately stopped. In this numerical example the gravitational load of the cylinder is included. The cylinder has same dimensions as in example 1. Additional data needed in the dynamic analysis is shown in Table 4. In addition we have changed the friction parameters slightly. In this example the static friction is scaled to the lifted mass. At time t ¼ 0 s when the initial state has been computed and the system is at equilibrium we give constant flow rate to the system for a period of one second. Then the inward flow rate is reset to zero. For the outward flow Q B we use the orifice flow model by Ellman et al. [16]. Time integration is carried out using constant time step of 0.001 s for the coupled system. Each time step converges on average with 3 iterations. In our computations we have not used any artificial damping like the Rayleigh damping. In this way we can compare the frictionless and frictional cylinder more precisely and observe clearly the dampening due to friction. For both models we note that the outbound flow rate is computed using the dissipative turbulent orifice model thus some damping is introduced for both of the models. The cylinder plus chamber pressures are presented in Fig. 7 where we have pressure curves for three different models: a

Table 4 Initial values for lifting of a point mass in addition to the values given in Table 2. Quantity

Symbol

Value

Unit

Gravitational acceleration Flow rate Hydraulic fluid density Steel density Cylinder lining mass Cylinder lining mass coeff.

g QA

9.81 10 870 7850 60 0.5

m/s2 l/min kg/m3 kg/m3 kg kg

v Str

40 0.5 0.7 F st 149.3 28.284 0 1

kpr

109

qf qs gp

Cylinder rod mass Cylinder lining mass coeff. Coulombic friction Stiffness coefficient Damping coefficient Viscous friction Stribeck velocity Pressure coefficient

gr F Cou k0 k1 Fv

N MPa kN/s m N/m m/s N/Pa

frictionless cylinder and two frictional cylinders with different friction coefficients. The dotted line presents the chamber pressure for the frictionless cylinder. The pressure in cylinder chamber is a bit higher than 20 MPa due to the gravity load and the back pressure in minus chamber. When the flow rate is lead to the cylinder the pressure in plus chamber rises roughly to 26 MPa and starts to oscillate around the initial pressure. This pressure oscillation settles slowly due to dissipative outward flow from the minus chamber. When the incoming flow rate is suddenly set to zero the pressure begins to oscillate around the stationary value. When friction is added to the cylinder the results are somewhat different. In its initial state the plus chamber pressure is lower due to the sealing friction. When a flow rate is lead to the cylinder the peak pressures are higher than in frictionless cylinder because the friction opposes the movement. When the friction coefficient is increased the pressure oscillation settles faster but the stationary value is of course higher than with the frictionless model. We also notice the phase shift between the models. This phase shift occurs because the friction slows the response of the system. In starting of the movement we can see from Fig. 8 that the friction dampens the behavior of the systems therefore it is somewhat slower. The initial displacement is higher with the frictionless cylinder than with the frictional cylinder as we observed in example 1. The chamber pressure is also lower because the friction carries part of the load. With high static friction the cylinder moves first due to the elastic deformation of the bristle and when the plus chamber pressure is high enough the static friction changes to sliding friction.

Chamber pressure [bar]

300 250 200 150 100

0

Fig. 6. Initial configuration for lifting of a point mass.

No friction Friction 0.1 Friction 0.25

50

0

0.2

0.4

0.6 0.8 Time [s]

1

1.2

1.4

Fig. 7. Plus chamber pressures of the lifting example.

70

A. Ylinen et al. / Computers and Structures 138 (2014) 62–72

Mass vertical displacement [mm]

10 8 6 4 2 0 No friction Friction 0.1 Friction 0.25

−2 −4

0

0.05

0.1

0.15 Time [s]

0.2

0.25

0.3

Fig. 8. Starting of the cylinder movement.

In stopping phase of the system we notice completely different response when the frictionless cylinder is compared with the frictional cylinder, see Fig. 9. When the frictionless system starts to oscillate around the stationary state, the frictional systems stop faster and there is practically no oscillation. We observed same kind of response in the chamber pressures, see Fig. 7. When the flow rate is reset to zero, the mass continues to rise because of the inertia. Thus in the frictionless model the kinetic energy conserves and the mass starts to oscillate with constant angular velocity around the static equilibrium. With frictional cylinder the response is different because the friction sticks and there is practically no oscillation. The only vibration is through the elastic deformation of the bristle. From this example we can see, that the frictionless and frictional dynamical cylinder model can be applied to simulations effectively. The convergence of the Newton’s iteration is quadratic and the results shown are sensible. In the engineering sense the friction offers damping and the cylinder movement stops in a realistic way. However, the frictionless model requires some sort of external damping model to settle the pressure oscillation and thus giving realistic results. 5.3. Raising boom In our final example we study the lifting boom configuration shown in Fig. 10. Movement of the boom is achieved with three different methods. In the first case we give a constraint equation for the distance between the nodes D and B as a function of time and use the Lagrange multiplier method to solve the system response. In the second case a length controlled rod element is used where the unstressed length follows the same linear function as in previous case. In the final case we use the hydraulic cylinder element. In this last case the incoming flow Q A is adjusted so that

Mass vertical displacement [mm]

42

41

the same movement is achieved as in the first cases assuming incompressible fluid. Finally the frictionless cylinder element results are compared to the results brought by the frictional elements with two different friction coefficients. The boom system without the hydraulic cylinder is modeled using the geometrically exact 3D total Lagrangian Reissner beam elements introduced by Mäkinen [17]. Similar beam elements are presented in [10,18]. The main dimensions are as shown in the Fig. 10 and additional data is shown in Table 5. To draw comparison between these three different methods for producing movement to the system we study the normal stress history in the cross section at the top edge of the node B, see Fig. 10. We are interested in the maximum values of the stress history. The distance between the nodes D and B is determined by assuming that the whole lifting movement of the boom is done in time period of 5 s in which time the boom rises from horizontal plane to the angle of 60°. The constant time step for the system is set to 0.001 s and each time step converges on average with 3 iterations. For the Lagrange multiplier method we construct a constraint equation where we define the distance of the cylinder attachment nodes as a linear function of time. The length controlled rod element is presented in [4]. Instead of the engineering strain in this example the element is based on the Green strain and the unstressed length of the element can be controlled as function of time. The same linear function as we used in the Lagrange multiplier method is used for the unstressed length of the element.

5.3.1. Results In Fig. 11 we see the normal stresses at the top edge of the cross section at node B compared with three different models: the constraint equation (CE), length controlled rod element (RVL) and frictionless cylinder element (NoFric). As expected the model with constraint equations solved by Lagrange multipliers gives the highest natural frequencies and normal stresses because the connection between the nodes B and D is assumed to be rigid. It also introduces some spurious oscillation to the system. The length controlled rod element evens out slightly the behavior of the system in comparison to the model with constraint equations. Although the first time step where the extension rate of the rod is not differentiable there is not a remarkable influence to the stress history. In the case of the hydraulic cylinder model we see that the behavior is smoother. In addition a slight reduction at the lowest natural frequency compared to the RVL model can be observed. In Fig. 12 the stress histories computed with three different friction coefficients are compared. When we add friction to the system the highest peak stress lowers and it also occurs a bit later which is known to be true from the natural frequencies of damped systems [19]. The observed peak

40

39

38

No friction Friction 0.1 Friction 0.25 1

1.1

1.2 1.3 Time [s]

1.4

1.5

Fig. 9. Mass displacement when the flow rate is set to zero.

Fig. 10. The initial configuration of the lifting boom example where the dash-dot – line represents the actuator and the lengths are given in meters.

A. Ylinen et al. / Computers and Structures 138 (2014) 62–72

400

a sophisticated friction model implemented which takes the key elements of the sliding friction between the cylinder lining and piston into account. The proposed element has several benefits compared to the traditional methods for modeling the linear hydraulic actuator. Firstly there are no constraint equations involved with the system thus no excessive variables are needed. In addition it is known, that the constraint equations can induce spurious oscillation to the system [3]. When we compared the hydraulic cylinder element with the length controlled rod element we observed that the cylinder element softens the system in starting and stopping of the movement. Moreover, the length change of the cylinder element is determined by the flow rates and the chamber pressures unlike in the length controlled rod element where we have to give explicit function for the length change. On the downside the cylinder element makes the differential equations stiffer and effective solution algorithms have to be found [20]. In the future research we shall take the flexural stiffness of the cylinder into account thus we can model the buckling of the hydraulic cylinder. The cylinder element is then modeled using beam elements and the relative sliding can be modeled using local constraints as presented in [21]. The sliding constraint can also be embedded to the element thus no Lagrange multipliers are needed [22]. Natural follow up is to extend the modeling to the hydraulic control system to compute accurate flow rates to the hydraulic cylinder instead of giving the flow as predefined function [23]. The hydraulic control system can be assembled as a network model consisting of transmission line models and therefore the flow rate to the cylinder depends on the state of the hydraulic control system [24]. Developing a new time stepping method is a valid research question because the coupled system is stiff and damped. Time stepping can be improved by using multi rate integration or by choosing an integrator specifically designed for stiff differential equations.

350

Acknowledgments

300

We would like to thank the Academy of Finland and the Finnish Cultural Foundation for the financial support of this work.

Table 5 Additional given values for raising boom example. Quantity

Value

Unit

Young’s modulus GAs Beam height Beam width Wall thickness Beam cross-sectional area Beam second moment of area

200 EA 0.12 0.08 0.005 0.0019

GPa N m m m m2 m4

3:756  106 14.92 400 18.94 0.3

Beam linear density Point mass at the boom end Rod adjusted Young’s modulus Calculation time tmax

Stress at point B [MPa]

kg/m kg GPa s

CE RVL NoFric

400

350

300

250

200

0

0.05

0.1

0.15 Time [s]

0.2

0.25

0.3

Fig. 11. Normal stresses at point B in the lift boom.

Stress at point B [MPa]

71

250

200

NoFric Friction 0.1 Friction 0.25 0

0.05

0.1

0.15 Time [s]

0.2

0.25

0.3

Appendix A. Supplementary data Supplementary data associated with this article can be found, in the online version, at http://dx.doi.org/10.1016/j.compstruc.2014. 02.006.

Fig. 12. Normal stresses at point B using different friction coefficients.

References stresses in simulation model with the frictional cylinder are smaller in comparison with the frictionless cylinder. In addition the peaks are not as steep. The comparison between the frictional cylinder and the constraint equation shows this behavior more clearly. This is explained partly by the damping factor of the friction model. The other part is through the turbulent orifice model used in the outbound flow rate because turbulent flow is energy consuming. 6. Conclusions In this paper we have presented a new hydraulic cylinder model for flexible multibody simulations. In the numerical simulations we have used a monolithic coupling to treat with the interaction between the hydraulic fluid in the cylinder chambers and the mechanical system. Thus no constraints are required. In addition,

[1] Haug EJ. Computer aided kinematics and dynamics of mechanical systems. Massachusetts: Allyn and Bacon; 1989. [2] Bauchau OA. Flexible multibody dynamics. Springer; 2001. [3] Marjamäki H, Mäkinen J. The raising movement of a hydro-mechanical lift frame, Rethymno, Greece; 13–16 June 2007. [4] Marjamäki H, Mäkinen J. Modelling a telescopic boom – the 3D case: part II. Comput Struct 2006;84:2001–15. [5] Wriggers P. Computational contact mechanics. 2nd ed. Berlin, Heidelberg: Springer; 2010. [6] Canudas de Wit C, Olsson H, Åström KJ, Lischinsky P. A new model for control of systems with friction. IEEE Trans Automat Control 1995;40:419–25. [7] Ibrahimbegovic A, Wilson EL. Unified computational model for static and dynamic frictional contact analysis. Int J Numer Methods Eng 1992;34:233–47. [8] Cardona A, Géradin M. Modeling of a hydraulic actuator in flexible machine dynamics simulation. Mech Mach Theory 1990;25. [9] Naya MA, Cuadrado J, Dopico D, Lugris U. An efficient unified method for the combined simulation of multibody and hydraulic dynamics: comparison with simplified and co-integration approaches. Arch Mech Eng 2011;58:223–43. [10] Géradin M, Cardona A. Flexible multibody dynamics: a finite element approach. Chichester: J. Wiley & Sons; 2001. [11] Pfeiffer F. Mechanical system dynamics. Berlin, Heidelberg: Springer; 2008.

72

A. Ylinen et al. / Computers and Structures 138 (2014) 62–72

[12] Ibrahimbegovic A, Mamouri S. Nonlinear dynamics of flexible beams in planar motion: formulation and time-stepping scheme for stiff problems. Comput Struct 1999;79:1–22. [13] Ibrahimbegovic A, Mamouri S. Energy conserving/decaying implicit timestepping scheme for nonlinear dynamics of three-dimensional beams undergoing finite rotations. Comput Methods Appl Mech Eng 2002;191:4241–58. [14] Bathe KJ. Finite element procedures. New Jersey: Prentice-Hall; 1996. 1039pp. [15] Hughes TJR. The finite element method, linear static and dynamic finite element analysis. Prentice-Hall; 1987. [16] Ellman A, Piché R. A two regime orifice flow formula for numerical simulation. ASME J Dyn Syst Meas Control 1999;121(4):721–4. [17] Mäkinen J. A total Lagrangian Reissner’s geometrically exact beam element without singularities. Int J Numer Methods Eng 2007;70:1009–48. [18] Ibrahimbegovic´ A, Frey F, Kozar I. Computational aspects of vector-like parametrization of three-dimensional finite rotations. Int J Numer Methods Eng 1995;38:3653–73.

[19] Thomson WT. Theory of vibration with applications. 4th ed. New Jersey: Prentice-Hall; 1993. [20] Hairer E, Wanner G. Solving ordinary differential equations II, stiff and differential-algebraic problems. Springer series in computational mathematics, vol. 14. Berlin: Springer; 1991. [21] Ibrahimbegovic A, Mamouri S. On rigid components and joint constraints in nonlinear dynamics of flexible multibody system employing 3D geometrically exact beam model. Comput Methods Appl Mech Eng 2000; 188:805–31. [22] Marjamäki H, Mäkinen J. Total Lagrangian beam element with C1 -continuous slide-spring. Comput Struct 2009;87:534–42. [23] Mäkinen J, Piché R, Ellman A. Fluid transmission line modeling using a variational method. ASME J Dyn Syst Meas Control 2000;122:153–62. [24] Pfeiffer F, Borchsenius F. New hydraulic system modelling. J Vib Control 2004;10:1493–515.