A comprehensive method for the contact detection of a translational clearance joint and dynamic response after its application in a crank-slider mechanism

Mechanism and Machine Theory 145 (2020) 103717

Contents lists available at ScienceDirect

Mechanism and Machine Theory journal homepage: www.elsevier.com/locate/mechmachtheory

Research paper

A comprehensive method for the contact detection of a translational clearance joint and dynamic response after its application in a crank-slider mechanism Mengbo Qian a,b, Zhe Qin b, Shaoze Yan a,∗, Lin Zhang a a b

State Key Laboratory of Tribology, Department of Mechanical Engineering, Tsinghua University, Beijing 100084, China Department of Mechanical Engineering, Zhejiang A&F University, Hangzhou 311300, China

a r t i c l e

i n f o

Article history: Received 21 October 2019 Revised 7 November 2019 Accepted 15 November 2019

Keywords: Translational clearance joint Nonlinear dynamic Detection Contact modes Experiment

a b s t r a c t The clearance in a translational joint affects the dynamic characteristics and leads to chaotic response of the mechanism. An experiment using a 3D translational joint with clearances was conducted to reveal the characteristics of the relative motion between the slide way and the slider. Then, a 3D model of translational joint with clearance was proposed. This model presents a comprehensive description of the contact modes in which a group of contact force models were employed to reflect different contact-impact phenomena. A method of contact detection, which aims to distinguish the positional relationship between two elements during the collision, was proposed according to the vector relationship of the apices of the contact segment between the slide way and the slider. To assess the nonlinear dynamic behaviour of a mechanism with translational clearance, the detection method was applied to numerical calculation and simulation of a crank-slider mechanism to illustrate its application, and investigate the nonlinear dynamics of this system. By combining theoretical calculation, simulation and experiment, the nonlinear dynamic characteristics of the mechanism with a clearance translational joint are revealed. © 2019 Elsevier Ltd. All rights reserved.

1. Introduction A mechanical system with translational clearance joints is well known as a nonlinear dynamic system. Much experimental and theoretical work has been done to study the issues associated with translational joints. Bauchau et al. [1] established a mechanical model with a translational joint on a fully-flexible micro-motion robot. Dupac et al. [2] analysed the translational joint of the flexible slider mechanism and the key structural parameters of the flexible translational joint were also analysed. Li et al. [3] studied the stiffness characteristics of a flexible parallel plate translational joint. Most studies of its dynamic characteristics are conducted under ideal conditions. Clearance will lead to nonlinear behaviour in the dynamic characteristics of the mechanism, thus affecting its dynamic performance [4–8]. Therefore, scholars have investigated mechanisms with clearances [9–16]. The clearance is necessary, but has its consequences. In the 1980s, Haines [17] and Ravn [18] obtained the error transfer coefficient regression equation for a plane linkage mechanism with clearances by a multiple linear regression method. Bu [19], Ambrósio [20], and Machado [21] found that the dynamic

∗

Corresponding author. E-mail address: [email protected] (S. Yan).

https://doi.org/10.1016/j.mechmachtheory.2019.103717 0094-114X/© 2019 Elsevier Ltd. All rights reserved.

2

M. Qian, Z. Qin and S. Yan et al. / Mechanism and Machine Theory 145 (2020) 103717

response of the system was affected by clearances, resulting in a deviation from expected performance and even serious accidents. Furthermore, experimental research into the dynamic performance of mechanical systems with clearance joints has been undertaken. Khemili designed and built an experimental set-up which consists of a slider–crank mechanism with clearance revolute joint to validate the simulation results [22]. To analyse the influences of clearance, crank speed, and crank length on the acceleration of a clearance rotating pair, Bengisu designed a crank slider mechanism test device with clearance between the connecting rod and rocker, and carried out comparative experiments by measuring the acceleration of the connecting rod and rocker in a horizontal direction [23]. Erkaya designed a crank-slider mechanism with two clearance pairs and installed three acceleration sensors to measure vibration [24]. The clearance effect of the kinematic pair was verified experimentally. At the same time, Flores designed a crank-slider mechanism with clearance between its kinematic pairs to validate the effect of clearance on the dynamic performance thereof. By changing the crank speed and the clearance of the kinematic pair, and based on the amplitude of the slider acceleration, he verified the influence of clearance on the dynamic performance of the mechanism [25]. Funabashi and Ogawa [26] reveal a four-bar mechanism with clearance to establish the clearance model of kinematic pairs under various conditions, and validated a previous theoretical model. Orden [27] tested the rigid and flexible mechanism of crank-slider mechanism with clearances, and measured the acceleration response to verify the theoretical collision model. Flores [28] conducted theoretical and experimental research into a crank-slider mechanism with clearances, and verified the theoretical research results thereby. Khemili [29] investigated the dynamic characteristics of a flexible linkage mechanism with clearances and carried out experimental verification, thus revealing three motion states of the mechanism with clearances: rigid and flexible linkages were simulated and compared. Erkaya [30,31] elucidated the influence of clearance between hinges on the dynamic characteristics of a crank-slider mechanism through experiments. The results show that the existence of clearance increases the vibration, noise amplitude, and joint impact force and moment on the structure: the larger the clearance, the greater the vibration and noise amplitude of the structure [32]. In addition, Muvengei et al. [33] obtained the non-linear dynamic response of a large stretching mechanism with articulated clearance through theoretical and experimental methods. Li et al. [34]designed and manufactured an experimental platform with a clearance mechanism. A piezoelectric acceleration sensor was used to measure the steady-state acceleration response of the rigid and elastic links of a mechanism with clearances and the experimental results validated the correctness of the theoretical modelling. Varedi-Koulaei et al. [35] studied the dynamic characteristics of space-deployable structures with articulated clearances in a microgravity environment by means of experiments. Jia et al. [36] carried out an experimental study on the three-ball pin pair mechanism with clearances and the experimental results are in good agreement with those obtained by theoretical calculation. Soong et al. [37] studied the motion state of the mechanism with clearances through experiments: they further refined the three-state model into four states, which made the relative motion relationship between the pair elements of the mechanism with clearances more precise. Haines [38] analysed the contact deformation characteristics of dry friction rotary clearance hinges without considering the lubrication between hinges through experiments: results show that the local contact deformation of the elements of clearance hinges is nonlinear. Besides, a mechanical system with clearance joints may exhibit chaotic responses under certain conditions, which are unique to nonlinear systems and are similar to random motion and, in the long-term, are unpredictable [39,40]. The analysis and research thereof in clearance joints is conducive to predicting the dynamic behaviour of mechanisms, especially in precision machinery [41,42]. To sum up, there is little research available on the contact force of translational joints considering complex contact modes when the clearance is considered. There is no systematic study of the nonlinear dynamics of the translational joints with clearances. Here, an experiment on a 3D clearance translational joint mechanism was conducted to reveal the characteristics of the relative motion between the slide way and the slider. From the experimental results, it can be demonstrated that relative motion between the slider and the slide way in a translational joint with clearances is complex. Based on evidence from the experiment, a 3D model of a translational joint with clearances was proposed. A comprehensive description of contact modes was developed in this model, in which different formulations of continuous contact forces were employed to describe different contact-impact phenomena. Furthermore, the detection method was applied to numerical calculation and simulation of a crank-slider mechanism to illustrate its application, and investigate the nonlinear dynamics of this system.

2. Experimental study of a mechanism with translational joint The experimental set-up with clearance translational joint is shown in Fig. 1. This set-up comprises a real translational joint with clearances which can be used to reveal the characteristics of the relative motion between the slide way and the slider. The structure of the translational joint is illustrated in Fig. 2, in which the slider way is fixed to the frame while the slider is constrained to move within the boundary thereof. In this experiment, two acceleration sensors, are respectively, fixed to the top and the front surface of the slider and make it move constantly. The driving device is an AC motor with 1.5 kW rated power, and the speed can be controlled and the inclination angle of the slide way can also be adjusted. The parameters used in the experiment are listed in Table 1. The clearance size is defined by the difference in width and height between the slider and the slide way. We can change the internal diameter of the slider to obtain a different clearance. The lifting handle can be used to adjust the deflection angle and the speed of the slider can be controlled by the speed regulating switch. It should be highlighted that the joint clearance size in actual mechanisms will not generally exceed 0.5 mm, therefore, when studying the influence of the clearance sizes, several key dimensions were selected in the range of 0–0.5 mm.

M. Qian, Z. Qin and S. Yan et al. / Mechanism and Machine Theory 145 (2020) 103717

3

Fig. 1. Experimental setup with clearance translational joint.

Fig. 2. The diagrammatic sketch. Table 1 Parameters used in the four cases. Num.

Clearance (mm)

Operating speed (cm/s)

Deflection angle (°)

1 2 3 4

0.25 0.15 0.25 0.25

6 6 9.3 6

0 0 0 30

The vibration of the mechanism is mainly due to the clearance, the operating speed and the deflection angle. In this section, the vibration of the mechanism is indirectly obtained by analysing the variation of the acceleration with the size of the clearance, the operating speed and the deflection angle, which can verify the complexity of nonlinear motion of translational clearance joint. Experimental results are shown in Figs. 3–5 which were obtained after steady state motion of the slider has been reached, where the representative result refers to general result, which is got from general mechanism with 0.25 mm clearance size, 6 cm/s operation speed and 0° deflection angle, it can be used for comparative analysis with results in other cases. Fig. 3 shows the time-domain signal of the acceleration: due to the presence of the clearance, the slider is accompanied by vibration during movement, so the acceleration of the slider fluctuates slightly around zero because collisions between the two elements are affected thereby. Fig. 4 shows a cumulative area diagram in different cases: the acceleration fluctuation periods in different cases differ from each other. The acceleration amplitude under the clearance case, the operating speed case, and the deflection angle case is greater from the perspective of the representative result. Fig. 5 shows the frequency-domain signals of acceleration: differences between the signals on Channels 1 and 2

4

M. Qian, Z. Qin and S. Yan et al. / Mechanism and Machine Theory 145 (2020) 103717

Fig. 3. Influences of different factors on acceleration. (a) representative result; (b) clearance; (c) operating speed; (d) deflection angle.

Fig. 4. Cumulative area diagram (acceleration) in different cases.

support the conclusions obtained from time-domain analysis. The signal on Channel 1 comes from the top face of the slider. The signal of channel 2 is from the front face of the slider. From the experimental results, it can be demonstrated that relative motion between the slider and the slide way in a translational joint with clearances is complex because misalignment may occur between the two elements and the slider can move outside the slide way boundary. Besides, the contact-impact phenomena are more complicated than in the idealised model, thus, a more precise model of a translational joint with clearances is required. 3. Translational joint model with clearance In this section, a model for 3D translational joints with clearances is proposed. According to the geometric description of a translational joint with clearances, kinematic variables associated with the translational joint were determined. Based on the results, different contact force models were used to describe the contact-impact phenomena, and the forces between the two joint elements could be calculated.

M. Qian, Z. Qin and S. Yan et al. / Mechanism and Machine Theory 145 (2020) 103717

5

Fig. 5. Effects of different factors on acceleration. (a) representative result; (b) clearance; (c) operating speed; (d) deflection angle. ∗ dB represents the reference value and g represents the acceleration due to gravity.

Fig. 6. Positional relationship between the slider and the slide way.

Fig. 7. A translational joint with clearance.

3.1. Composition of translational joint with clearance Translational joints of rectangular cross-section can be simplified into two elements: a slide way and a slider. The slide way is simplified as a solid cuboid, and the slider is simplified as a hollow cuboid. The translational joint can be expressed as the structure shown in Figs. 6 and 7. Generally, in translational joint mechanisms, the slide way is connected to other components in the mechanical structure. Considering that the slider does not move to both ends of the slide in the actual working process, therefore, this paper will not discuss the contact between the slide and the slider along the moving axis.

6

M. Qian, Z. Qin and S. Yan et al. / Mechanism and Machine Theory 145 (2020) 103717

Fig. 8. Contact-segment of the translational joint.

3.2. Mathematical model of a translational joint Marking the slide way and the slider as L and K, respectively, in two coordinate systems, Ol − el1 el2 el3 and Ok − ek1 ek2 ek3 , whose origins are at the centroids of L and K, as shown in Fig. 7, the coordinate axes are perpendicular to the planes of L and K. The contact-segment can be extracted for study when L and K come into contact, so the contact mode of the two elements can be distinguished upon detection of the positional relationships between the apices, therefore, for ease of analysis, the contact-segment of the slider is simplified into two rectangular steel plates joined by a rigid link, and the apices are marked as point P1 to point P8 , as shown in Fig. 8. In this figure, the solid line represents the slide way and the dotted line represents the slider. 3.2.1. Detection of contact modes In a spatial mechanism, translational joints with clearance not only have different contact modes, but also different attitudes under each contact mode. Zhang et al. proposed a feasible method for determining the position of the slider [43]. For various contact modes of the translational joint, the contact force between the slider and the slide way is not constant. To distinguish the contact modes from the possible relative positions, it is necessary to quantitatively detect each contact mode on the translational joint on the basis of the qualitative division so that the contact force calculation can be performed. The function of SYSARY is a function subroutine, which is a special subroutine provided by ADAMS. It provides system state values such as displacement and acceleration to the user subroutine. The function of SYSARY is called to get the state value of the system when the model is under simulation. The position and velocity vectors of point P1 –P8 in contact with the slider can be identified and are provided to the user sub-routine. Take point P1 as an example: Suppose that, in inertial coordinate system O − i1 i2 i3 , the position vectors of the coordinate system of the slide way and −−→ −−→ −−→ slider in the inertial system are ul and uk , respectively, Ol Ok is d, Ol P1 is r1 , Ok P1 and is R1 . In coordinate system Ol − el1 el2 el3 , −−→ Ol P1 is r1l . Al and Ak are the coordinate transformation matrices of the relative inertial system of the connected coordinate system of the slide way and the slider, respectively. The position vector R1 can be given by

R1 = −d + r1 = ul − uk + Al r1l

(1)

Regardless of the relative position of the slide and slider in the direction of ek3 , the contact situation between the two elements at point P is as follows: Suppose the position vector of P1 is R1k , and the location in the Ok − ek1 ek2 ek3 coordinate system can be expressed as



R1k = ATk R1 = ATk ul − uk + Al r1l



(2)

The translation transformation matrix and the rotation matrix can be expressed as:

⎡

1 ⎢0 Trans(x, y, z ) = ⎣ 0 0

⎡

cos θ ⎢ 0 Rot(y, θ ) = ⎣ − sin θ 0

0 1 0 0 0 1 0 0

0 0 1 0 sin θ 0 cos θ 0

⎤ ⎡ x 1 y ⎥ ⎢0 ; Rot(x, θ ) = ⎣ z ⎦ 0 1

⎤

0

⎡

0 cos θ sin θ 0

0 cos θ 0⎥ ⎢ sin θ ; Rot(z, θ ) = ⎣ 0⎦ 0 1 0

0 − sin θ cos θ 0

− sin θ cos θ 0 0

0 0 1 0

⎤

0 0⎥ ; 0⎦ 1

⎤

0 0⎥ , 0⎦ 1

The composite transformation matrix Al and Ak can be expressed as:

Al = Rot(z, θl )Rot(y, θl )Rot(x, θl )Trans(xl , yl , zl )

(3)

Ak = Rot(z, θk )Rot(y, θk )Rot(x, θk )Trans(xk , yk , zk )

(4)

M. Qian, Z. Qin and S. Yan et al. / Mechanism and Machine Theory 145 (2020) 103717

7

Table 2 Detection of contact situations. Points

Relational expressions

P1

R1k R1k R2k R2k R3k R3k R4k R4k R5k R5k R6k R6k R7k R7k R8k R8k

P2 P3 P4 P5 P6 P7 P8

· · · · · · · · · · · · · · · ·

ek1 ek2 ek1 ek2 ek1 ek2 ek1 ek2 ek1 ek2 ek1 ek2 ek1 ek2 ek1 ek2

> a /2 > b /2

Contact surfaces P1 P1 P2 P1 P2 P4 P1 P4 P1 P1 P2 P1 P2 P4 P1 P4

< −a /2 > b /2 < −a /2 < −b /2 > a /2 < −b /2 > a /2 > b /2 < −a /2 > b /2 < −a /2 < −b /2 > a /2 < −b /2

− P4 − P2 − P3 − P2 − P3 − P3 − P4 − P3 − P4 − P2 − P3 − P2 − P3 − P3 − P4 − P3

− P8 − P6 − P7 − P6 − P7 − P7 − P8 − P7 − P8 − P6 − P7 − P6 − P7 − P7 − P8 − P7

− P5 − P5 − P6 − P5 − P6 − P8 − P5 − P8 − P5 − P5 − P6 − P5 − P6 − P8 − P5 − P8

Table 3 Surface contact.

Order

Position

1

Illustration

Order

In contact with one surface

2

Position

Illustration

In contact with two surfaces

∗ 1. In the table above, black solid lines represent the inner border of the slider, the black dashed line represents the rectangular plate 1, and the grey dashed lines represent the rectangular plate 2. ∗ 2. "•" and "•" represent the contact points of plate 1 and plate 2, respectively.

When R1k · ek1 > a /2, P1 is in contact with the top face P1 − P4 − P8 − P5 of the slide way; when R1k · ek2 > b /2, P1 represents in contact with the front face of the slide way. The relative velocity vector of point P1 and the slider is δ˙ = R˙ 1k . When the two elements are in contact, normal contact force acts in the direction of ek1 or ek2 . We can reach following conclusion in the same way and detection of the contact situation is as listed in Table 2. The contact mode between the two elements can be determined after identifying the contact of eight possible contact points through the aforementioned steps. Only in this way can we determine the kind of contact type the two elements are under surface contact, line contact, or point contact. 3.2.2. Detection of surface contact The apices of the contact-segment of the slide way are marked P1 to P8 , therefore, the four outer surfaces of L can be marked as follows: face 1:P1 − P4 − P8 − P5 , face 2:P2 − P3 − P7 − P6 , face 3: P1 − P2 − P6 − P5 , face 4: P4 − P3 − P7 − P5 , the four inner surfaces of K are marked as face 1 :P1 − P4 − P8 − P5 , face 2 :P2 − P3 − P7 − P6 , face 3 : P1 − P2 − P6 − P5 , face 4 :P4 − P3 − P7 − P8 . According to the method provided above, we can determine if the four apices in each surface are in contact with K separately. When the four apices on one surface are in contact with the inner surface of the slider, it can be determined that L and K are in surface contact, and the contact surface is determined by the four apices. Possible positions between surfaces L and K (when in contact) are listed in Table 3. Due to the existence of clearance, there will be friction-wear phenomena on the contact surfaces of the two elements during the movement of the mechanism, which is caused by frictional force and will result in high temperature. Regardless of wear and tear, lubrication and other factors, the interaction force between the two elements can be decomposed into two parts, normal contact force FN and tangential contact force Ft (namely the friction), respectively. The normal contact force should be calculated according to three different contact modes of two elements, namely the point contact mode, line contact mode, and surface contact mode. The tangential friction force Ft can be calculated as a Coulomb friction after obtaining the normal contact force FN . the tangential friction force is expressed as

Ft = μFN where μ is the coefficient of sliding friction.

(5)

8

M. Qian, Z. Qin and S. Yan et al. / Mechanism and Machine Theory 145 (2020) 103717 Table 4 Line contact.

Order

Position

Illustration

Order

Position

Illustration

1

In contact with one line

2

In contact with one line

3

In contact with two lines

4

In contact with two lines

5

In contact with two lines

6

In contact with two lines

7

In contact with three lines

8

In contact with four lines

9

In contact with four lines

10

In contact with one line and one plane

Considering the stick-slip effect on the contact surface and to avoid the sudden change of friction caused by the change of the speed direction, the friction coefficient is modified to [44,45]

⎧ ⎨μd ,   |V | μ = μS sin π2 Vs , ⎩ μs +μd + 1 (μ − μ ) cos π |V |−Vs , S d 2 2 Vd −Vs

|V | > Vd |V | < VS Vs ≤ |V | ≤ Vd

(6)

where Vs and Vd are the stick-slip switch velocity and static-sliding friction switch velocity, respectively, and μs and μd are the static and sliding friction coefficients, respectively. For the translational joint with clearance, the contact area is that rectangular surface when L and K are in surface contact. According to the mechanics of materials [46], the normal contact force FN (δ  )between L and K is expressed as

FN (δ ) =

EAδ  l

(7)

where E is the elastic modulus of the material, A is the contact area, δ is the distance through which the material is compressed, and l denotes the thickness of the material in the principal direction of compression. 3.2.3. Detection of line contact There are 12 boundary lines in the contact-segment (P8 P4 , P5 P8 , P5 P1 , P1 P4 , P5 P6 , P8 P7 , P4 P3 , P1 P2 , P6 P7 , P6 P2 , P2 P3 , and P3 P7 ). When two apices on any boundary line are in contact with one of the surfaces of the slider, L and K are in line contact. Possible positions between L and K (when in line contact) are listed in Table 4. For a translational joint with clearance, the contact region can be regarded as a narrow rectangle when L and K are in line contact. The contact mode is a cylinder-plane contact. Setting the contact area to be a rectangle with length l and width 2a, the distribution of the contact force F presents the shape of a semi-elliptical cylinder [47], as given by:

F (x ) =

2P

π a2 l



a2 − x2

(8)

where x is the radial distance of the contact point from the centre of the contact area, and P is the pressure per unit length.

M. Qian, Z. Qin and S. Yan et al. / Mechanism and Machine Theory 145 (2020) 103717

9

Puddock derived an expression for the penetration depth of a cylinder under pressure δ as approximately equal to [48]



δ=

P 1 − In π E∗l

R=

r2 r1 r2 − r1

 PR  π E∗l3

(9) (10)

where r1 and r2 are the equivalent contact radius of L and K, respectively. To compute the pressure P, the function of “Lambert W” is introduced, which is defined as: If y · ey = x, y = lambertw(x ), where x ∈ [− 1e , +∞ ). 

e

In

− −



PR

= e[1−

π E∗ l3

δπ E ∗ l P

δπ E ∗ l P

· e−

δπ E ∗ l

δπ E ∗ l P

P

]

=−

(11)

δR el 2



= l ambertw −

(12)

δR



el 2

(13)

So the pressure P is

P=−

δπ E ∗ l 

lambertw − eδlR2



(14)

When x = 0, the radial contact force can be expressed as

2P FN = = π al



PE ∗ π Rl

(15)

where R is the equivalent radius, and E∗ is the equivalent elastic modulus.

E∗ =

E 1 − v2

 a=

4P R

π E∗l

(16)

(17)

3.2.4. Detection of point contact After the above identification, if there are other apices in contact with K, it means that L and K are in point contact. Possible positions between L and K are shown in Table 5 when they are in point contact. For the translational joint with clearance, the contact point is one point on the edge of the contact-segment when L and K are in point contact. We apply a nonlinear spring damping model as is often used elsewhere and the normal contact force can thus be expressed as [49]:

˙ FN = K δ e + C δδ

(18)

where K is the Hertz contact stiffness, C is the damping factor, δ is the penetration depth at the contact points, and δ˙ is the normal relative velocity of the contact points. For metals, the Hertz contact force index e is generally between 1.3 and 1.5. In the nonlinear spring damping model, the damping coefficient C is generally 0.1 to 1% of the equivalent contact stiffness K [50]. It should highlight that Eq. (18) is the Lankarani–Nikravesh contact-force model which accounts for not only the elastic nature of the colliding bodies but also the energy dissipation during impact. However, the coupling relationship between penetration depth and contact stiffness is neglected in L-N contact force model, which is only suitable for large clearance and small load case. At the same time, a large recovery coefficient will lead to failure of the model, since the scope of the damping coefficient is limited by the recovery factor. So there are some problems when starting and ending the impact. 4. Case study To further study the nonlinear dynamic behaviour of mechanism with translational clearance, the 3D model and the detection method was applied to numerical calculation and simulation of a crank-slider mechanism to illustrate the application of the model presented, and investigate the nonlinear dynamics of this system. The model and method can reflect different contact-impact phenomena of the clearance translational joint more accurately. Therefore, the simulated results are closer to the actual values.

10

M. Qian, Z. Qin and S. Yan et al. / Mechanism and Machine Theory 145 (2020) 103717

Table 5 Point contact.

Order

Position

Illustration

Order

Position

Illustration

1

In contact with one point

2

In contact with two points

3

In contact with two points

4

In contact with three points

5

In contact with four points

6

In contact with one point and one line

7

In contact with two points and two lines

8

In contact with two points and one line

Fig. 9. Slider-crank mechanism with clearance translational joint.

Table 6 Geometric and inertial properties of the slider-crank mechanism. Rod 1 2 3 4

Length (mm) 50 142.2 40 250

Width (mm) 5 14.2 20 20

Depth (mm) 2.5 7.1 30 30

Density (kg/m3 ) 7.8 7.8 7.8 7.8

× × × ×

3

10 103 103 103

Elastic modulus (N/m2 ) 2.07 2.07 2.07 2.07

× × × ×

11

10 1011 1011 1011

Poisson’s ratio 0.29 0.29 0.29 0.29

M. Qian, Z. Qin and S. Yan et al. / Mechanism and Machine Theory 145 (2020) 103717

11

Table 7 Parameters. Notation

Parameter

Notation

Parameter

F IM JM CM1 CM2 CM3 CM4 CM5 CM6

Friction flag Marker point of L Marker point of K Marker for point C1 Marker for point C2 Marker for point C3 Marker for point C4 Marker for point C5 Marker for point C6

CM7 CM8 Ak Bk Lk E1 E2 μ1 μ2

Marker for point C7 Marker for point C8 Length a of the rectangular hole Width of the rectangular hole Height of the rectangular hole Elastic modulus of the material of the slider Elastic modulus of the material of the slide way Poisson’s ratio of the material of the slider Poisson’s ratio of the material of the slide way

Table 8 Parameters of the inner hole of the slider. Cases

Length (mm)

Width (mm)

Depth (mm)

speed (rad/s)

1 2 3 4

40 40 40 40

2.016 2.02 2.06 2.1

3.016 3.02 3.06 3.1

π π π π

4.1. Modelling of a slider-crank mechanism Fig. 9 shows the configuration of the slider-crank mechanism. The geometric and physical parameters of the mechanism are listed in Table 6. In this mechanism, translational joint between the slider and the slide way is considered as a 3D clearance joint, while other joints are assumed to be ideal. In the dynamic simulation, the direction of action of gravity is taken as the negative Y-direction and the motion of the system is defined in a vertical plane. 4.2. Numerical calculation and simulation The contact force model of a clearance translational joint was established based on the contact detection of the two elements, which is written as a user subroutine in FORTRAN language. In the case, the contact force was solved by calling a sub-routine according to the system state values obtained by the function subroutine “SYSARY” in ADAMS, and the results are returned to the translational joint model to complete the dynamic simulation analysis. 4.2.1. Specification of the parameters User parameters should be input to the subroutine first so as to get the simulation result. All input parameters are listed in Table 7. Contact forces between the two elements are calculated using the subroutine based on user-input parameters. 4.2.2. Solution procedure R The contact force is obtained by using software including Fortran, MATLAB , and ADAMS. The Lambert W function R in MATLAB is called to calculate the contact force model (Fig. 10). We can obtain the contact force between the two elements of the translational joint by calling the subroutine according to the flow chart. 4.2.3. Contact force under different contact modes Results shown in Fig. 11 are obtained from the mechanism with a clearance size of 0.15 mm at an operating speed of π rad/s. As shown in Fig. 11, the contact force, under different modes, is different at any time, as represented by taller columns on the histogram. This is because the magnitude of the contact force is related to the position of the point of contact between the two elements. We can also find that the contact force undergoes irregular changes, consistent with the variation of contact force as studied previously. It can be demonstrated that the contact force between the slider and the slide way in a translational joint with clearances is complex, because misalignment occurs between the two elements, the contact-impact phenomena are more complicated than in the idealised model. 4.3. Effect of clearance on kinematic characteristics To understand how clearance affects the dynamics of the slider-crank mechanism, as presented in Fig. 12, we focus on the velocity and acceleration of the slider in the z-direction and the angular velocity of the connecting rod. The parameters of the slider used in different cases are listed in Table 8 and the simulation results are shown in Fig. 12. It should be highlighted that the clearance is determined by the parameters of the slider: we obtain different clearances by changing the size of the inner hole. As shown in Fig. 12 and Table 8, the slider shows low kinematic characteristics corresponding to a relatively small clearance under gravity and connecting rod forces. Upon incorporating clearance, in addition to higher

12

M. Qian, Z. Qin and S. Yan et al. / Mechanism and Machine Theory 145 (2020) 103717

Fig. 10. The flow chart for the calculation.

velocity in the Z-direction, the acceleration is also found to gradually increase with increasing clearance, which can be verified with the experimental results shown in Fig. 5. We find similar regularities by comparing the angular velocity curves shown in Fig. 12(c): the bigger the clearance, the greater the angular velocity. It can be shown that clearances can lead to the vibration of the slider and misalignment, which changes the kinematic state of the system. 4.4. Nonlinear phenomena pertinent to the slider position The position of the slider in the vertical direction is constant in an ideal mechanism, the coordinate position of the slider in the Z-direction will fluctuate due to the clearance. Fig. 13 shows the nonlinear behaviour of the slider position, which refers to the nonlinear change of the slider position caused by the clearance. As presented in Fig. 13(a) and (b), the positions in the y and z-direction are irregular with a size of 0.1 and 0.5 mm at an operating speed of π rad/s. In addition to strong vibration bands of the slider at 0.5 rad/s and 1 rad/s, indicative of the presence of clearance and this is further evinced by the vibration peak located at 0.832 and 1.323 rad/s. The greater the clearance, the greater the peak value of the slider position in the Y-direction. Fig. 13(c) and (d) reflects the relationship between the positions in the y and z-directions at operating speeds of 1.22 and 2π rad/s under a clearance of 0.1 mm: the higher the speed, the higher the amplitude of slider motion. It is reasonable to conclude that the larger the clearance and the higher the speed, the more obvious the

M. Qian, Z. Qin and S. Yan et al. / Mechanism and Machine Theory 145 (2020) 103717

Fig. 11. Contact force under different contact modes.

Fig. 12. Effect of clearance size on kinematic characteristics. (a) velocity; (b) acceleration; (c) Angular velocity.

13

14

M. Qian, Z. Qin and S. Yan et al. / Mechanism and Machine Theory 145 (2020) 103717

Fig. 13. Nonlinear behaviour of the slider position (a) c = 0.1 mm, π rad/s; (b) c = 0.5 mm, π rad/s (c) c = 0.1 mm, ω = 1.22π rad/s; (d) c = 0.1 mm, ω = 2π rad/s.

fluctuations and the higher the slider. It can be shown that the slider position in the Y-direction is related to the contact force between the two elements. The clearance and the speed can lead to changes in microscopic positional relationships between the two elements, so the position where the contact force changes abruptly varies due to different clearances and speeds, which can lead to the more obvious fluctuations.

4.5. Effect of different factors on contact force As presented in Fig. 14(a) and (b), the contact forces on the slider are irregular at clearances of 0.1 and 0.5 mm at an operating speed of π rad/s. In addition to strong contact force bands of the slider at 0.3 and 0.8 s, indicative of the presence of clearance, this is further evinced by the peak contact forces located at 0.4 and 0.75 s. The greater the clearances, the greater the magnitude and amplitude of the contact force fluctuations. The contact forces show good symmetry with sliding time. The main reason for this is that vertical fluctuation of the slider is symmetrical, and the magnitude of the contact force is consistent with the position of the centre of the slider. In addition to the clearance, the speed affects the dynamic behaviour of the slider position. In this section, the crank speed is increased to 1.22π rad/s and 2π rad/s, respectively, and the clearance remains 0.1 mm. Fig. 15(c) and (d) reflects the relationship between the contact force and speed: the higher speed, the greater the magnitude and amplitude of fluctuation of the contact force. Such dynamic behaviour is primarily due to the contact force which may be caused by collisions in the system, so the contact force upon collision changes abruptly due to different clearances and speeds. This fluctuation of the contact force demonstrates that the clearances and speed can lead to chaotic motion of the mechanism via collision interactions, which can lead to dynamic behaviour as manifest in the slider position.

M. Qian, Z. Qin and S. Yan et al. / Mechanism and Machine Theory 145 (2020) 103717

15

Fig. 14. Effect of different factors on contact force (a) c = 0.1 mm, π rad/s; (b) c = 0.5 mm, π rad/s; (c) c = 0.1 mm, ω = 1.22π rad/s; (d) c = 0.1 mm, ω = 2π rad/s.

Fig. 15. The experimental and numerical results in clearance case. (a) experimental results; (b) numerical results.

4.6. Comparison of experimental results and numerical results To compare the experimental and numerical results, the acceleration of the slider in a case with clearance as found from experiment and simulation is shown in Fig. 15. Since the contact mode of the two elements switches between point contact, line contact, and surface contact, the direction of the acceleration changes cyclically as a result of vibration of the mechanism caused by clearance, the amplitude and frequency of the vibration of the mechanism are different under different clearances.

16

M. Qian, Z. Qin and S. Yan et al. / Mechanism and Machine Theory 145 (2020) 103717

Z

S

r P0

P1

P2 P3

P(x)

P4

x

S O Poincare section

r

Y

X Fig. 16. Schematic diagram of Poincaré map.

Fig. 17. Poincare maps for different radial clearance sizes. (a) Rc =0.08 mm; (b) Rc =0.1 mm; (c) Rc =0.3 mm; (d) Rc =0.5 mm.

From the experimental results shown in Fig. 15(a), the amplitude and frequency of the acceleration vary under different clearances, and the acceleration curve is serrated. On the whole, the acceleration increases with bigger clearance, except the 0.22 mm clearance case. Due to the limitations of the experimental conditions, it is impossible to identify all the factors affecting the results and find the deviation in the motion characteristics of the mechanism caused by clearance to sufficient accuracy, therefore, the contact force model of the clearance translational joint was established and numerical simulation was conducted. It can be verified in Fig. 15(b) that the clearance causes the mechanism to vibrate at a higher frequency, and the bigger the clearance, the greater the magnitude of the acceleration.

M. Qian, Z. Qin and S. Yan et al. / Mechanism and Machine Theory 145 (2020) 103717

17

Fig. 18. Flow chart of small-data sets method.

Fig. 19. Estimation of maximum Lyapunov exponent. (a) linear fitting region; (b) original result.

4.7. Detection of chaotic phenomena A mechanical system with translational clearance joints is a nonlinear dynamic system that may exhibit chaotic responses under certain conditions [51]. Chaos detectionmethods are mainly divided into qualitative analysis and quantitative analysis. In this paper, Poincaré mapping method is used to quantitatively analyse the chaotic characteristics of the mechanism while the Lyapunov exponential method is applied for qualitative analysis. 4.7.1. Poincaré map The system state variables of continuous systems change with time, and the periodic trajectories in the state space intersect with a lower dimension in a series of points. This is the Poincaré map and its space is named Poincaré section.This method is characterized by the transformation of the trajectory of the power system into the intersection of the trajectory and a cross section, so that the evolution of the system morphology over time can be observed on the Poincaré section. Fig. 16 is the schematic diagram of this method. The trajectory of a three-dimensional space intersects the cross-section S (z is constant) on a revolution direction (z<0) at discrete points P0 , P1 , P2 ,…, which form a Poincaré map.

Pk+1 = T (Pk ) = T (T (Pk−1 ) ) = T 2 (Pk−1 ) = · · ·

(18’)

Fig. 17 is a Poincaré map showing the displacement and velocity components of the slider in X direction. The dynamic responses were obtained for 150 full crank rotations, and the displacement and velocity at 2 s intervals during each cycle were selected for analysis. In ideal mechanical systems, the solutions will be the same in each period while they are different in the mechanism with clearance due to the development of chaotic behaviour. As can be seen from the figure, the map of the mechanism is disorganized, while mapping points are distributed in the figure, indicating that the mechanism has no periodic solution and has certain fractal dimension characteristics, which illustrates that the mechanism is in a chaotic state. At the same time, the distribution of mapping points vary with the change of clearance value, which indicates

18

M. Qian, Z. Qin and S. Yan et al. / Mechanism and Machine Theory 145 (2020) 103717

that the complexity of the dynamic response of the system also changes. Observably, the larger the clearance value, the wider distribution of the mapping points. 4.7.2. Lyapunov exponent Lyapunov exponent can be used to evaluate the sensitivity of dynamic system to initial values [52,53]. A positive exponent means that the system has chaotic characteristics. This is a quantitative detection method for chaotic phenomena. It is cumbersome and unnecessary to obtain all the Lyapunov exponents of the system. Most applications only need to solve for maximum Lyapunov exponent (Lym ). In this paper, Lym is estimated by means of small-data sets method. Fig. 18 is the flow chart of the method. In this paper, the time series of the displacement component of the slider along the X direction is selected for chaos analysis. The time delay (tau), time window(tw), and embedding dimension(m) can be obtained by C–C method [54], the average period(P) can be obtained by FFT method [55].

tau = 8, tw = 80, m =

tw + 1 = 11, P = 0.2147 tau

(19)

These data are substituted into small data volume algorithm, and the calculated result is shown in Fig. 19, where i represents time series and dv represents dispersion degree. Select the linear region shown in the dotted bordered rectangle and fit straight line in this region by least squares method. The slope of the fitted line is the maximum Lyapunov exponent [56]. Obviously, Lym > 0, it is further verified that the translational joint with clearance has chaotic characteristics. 5. Conclusion An experiment on a 3-D translational joint with clearances was conducted to reveal the characteristics of the relative motion between the slide way and the slider. Based on evidence from the experiment, a comprehensive description of contact modes was developed in this model, in which different formulations of continuous contact forces were used to describe different contact-impact phenomena. Furthermore, a numerical simulation was performed and the nonlinear dynamic behaviour of a crank-slider mechanism was investigated using this method. Through the numerical calculation and simulation, we analysed the influences of clearance on position, kinematics, and contact force of a translational joint mechanism: it can be concluded that a larger clearance and higher operating speed may lead to greater contact force and more fluctuations of the slider. At the same time, larger clearances may cause increases in mechanical characteristics of the mechanism, such as speed, acceleration, and angular velocity of the slider. We also verified the presence of chaotic behaviour in the translational joint with clearances by use of the Poincaré map and the Lyapunov exponent. It provides a method for the quantitative analysis of nonlinear dynamic characteristics of the mechanism with a clearance translational joint. Declaration of Competing Interest The authors declare that they have no conflict of interest. Acknowledgments This work is supported by National Science Foundation of China under Contract nos. 11872033, 51875531, and 51505432, and Natural Science Foundation of Beijing under Contract no. 3132017. Informed consent Informed consent was obtained from all individual participants included in the study. References [1] O.A. Bauchau, J. Rodriguez, Modeling of joints with clearance in flexible multibody systems, Int. J. Solids Struct. 39 (1) (2002) 41–63. [2] M. Dupac, D.G. Beale, Dynamic analysis of a flexible linkage mechanism with cracks and clearance, Mech. Mach. Theory 45 (12) (2010) 1909–1923. [3] Y.Y. Li, G.P. Chen, D.Y. Sun, Y. Gao, K. Wang, Dynamic analysis and optimization design of a planar slider-crank mechanism with flexible components and two clearance joints, Mech. Mach. Theory. 99 (2016) 37–57. [4] S.X. Wang, Y.H. Wang, B.Y. He, Dynamic modeling of flexible multibody systems with parameter uncertainty, Chaos Solitons Fractals 36 (3) (2008) 605–611. [5] K.L. Ting, K.L. Hsu, J. Wang, Clearance-induced position uncertainty of planar linkages and parallel manipulators, J. Mech. Rob. 9 (6) (2017) 061001. [6] X. Zheng, F. Zhang, Q. Wang, Modeling and simulation of planar multibody systems with revolute clearance joints considering stiction based on an LCP method, Mech. Mach. Theory 130 (2018) 184–202. [7] W.W.K. Xiang, S.Z. Yan, J.N. Wu, X.G. Robert, Complexity evaluation of nonlinear dynamic behavior of mechanisms with clearance joints by using the fractal method, Proc. Inst. Mech. Eng. C: J. Mech. Eng. Sci. 228 (18) (2014) 3482–3495. [8] C.S. Liu, K. Zhang, L. Yang, The compliance contact model of cylindrical joints with clearances, Acta Mech. Sin. 21 (5) (2005) 451–458. [9] J. Zhang, Q. Wang, Modeling and simulation of a frictional translational joint with a flexible slider and clearance, Multibody Syst. Dyn. 38 (4) (2016) 367–389. [10] L.X. Yu, C.S. Liu, Dynamic simulation and kinetic description of revolute joint with spatial clearance, Acta Sci. Natural. Univ. Pekinesis. 41 (5) (2005) 679–687.

M. Qian, Z. Qin and S. Yan et al. / Mechanism and Machine Theory 145 (2020) 103717

19

[11] K.L. Ting, K.L. Hsu, Z. Yu, J. Wang, Clearance-induced output position uncertainty of planar linkages with revolute and prismatic joints, Mech. Mach. Theory 111 (2017) 66–75. [12] X.D. Zheng, J. Li, Q. Wang, Q.M. Liao, A methodology for modeling and simulating frictional translational clearance joint in multibody systems including a flexible slider part, Mech. Mach. Theory 142 (2019) 103603. [13] Q. Tian, P. Flores, H.M. Lankarani, A comprehensive survey of the analytical, numerical and experimental methodologies for dynamics of multibody mechanical systems with clearance or imperfect joints, Mech. Mach. Theory 122 (2018) 1–57. [14] F. Marques, F. Isaac, N. Dourado, P. Flores, An enhanced formulation to model spatial revolute joints with radial and axial clearances, Mech. Mach. Theory 116 (2017) 123–144. [15] L. Skrinjar, J. Slavicˇ , M. Boltežar, A validated model for a pin-slot clearance joint, Nonlinear Dyn. 88 (1) (2017) 131–143. [16] X.Y. Geng, M. Li, Y.F. Liu, W. Zheng, Z.J. Zhao, Non-probabilistic kinematic reliability analysis of planar mechanisms with non-uniform revolute clearance joints, Mech. Mach. Theory 140 (2019) 413–433. [17] R.S. Haines, Survey: 2-dimensional motion and impact at revolute joints, Mech. Mach. Theory 15 (5) (1980) 361–370. [18] P. Ravn, A continuous analysis method for planar multibody systems with joint clearance, Multibody Syst. Dyn. 2 (1) (1998) 1–24. [19] W.H. Bu, Z.Y. Liu, J.R. Tan, S. Gao, Detachment avoidance of joint elements of a robotic manipulator with clearances based on trajectory planning, Mech. Mach. Theory 45 (6) (2010) 925–940. [20] P. Flores, J. Ambrósio, J.C.P. Claro, H.M. Lankarani, Spatial revolute joints with clearances for dynamic analysis of multi-body systems, Proc. Inst. Mech. Eng. K: J. Multibody Dyn. 220 (4) (2006) 257–271. [21] M. Machado, P. Flores, J.C.P. Claro, J. Ambrósio, M. Silva, A. Completo, H.M. Lankarani, Development of a planar multibody model of the human knee joint, Nonlinear Dyn. 60 (3) (2010) 459–478. [22] I. Khemili, L. Romdhane, Dynamic analysis of a flexible slider-crank mechanism with clearance, Eur. J. Mech. A Solids 27 (5) (2008) 882–898. [23] M.T. Bengisu, T. Hidayetoglu, A. Akay, A theoretical and experimental investigation of contact loss in the clearances of a four-bar mechanism, Strahlentherapie 108 (2) (1986) 237–244. [24] S. Erkaya, I. Uzmay, Experimental investigation of joint clearance effects on the dynamics of a slider-crank mechanism, Multibody Sys. Dyn. 24 (1) (2010) 81–102. [25] P. Flores, C.S. Koshy, H.M. Lankarani, J. Ambrósio, C.P. Claro, Numerical and experimental investigation on multibody systems with revolute clearance joints, Nonlinear Dyn. 65 (4) (2011) 383–398. [26] H. Funabashi, K. Ogawa, M. Horie, H. Lida, A dynamic analysis of the plane crank and rocker mechanisms with clearances, Bull JSME 23 (177) (1980) 446–452. [27] J.C.G. Orden, Analysis of joint clearances in multibody systems, Multibody Sys. Dyn. 13 (4) (2005) 401–420. [28] P. Flores, J. Ambrósio, Revolute joints with clearance in multibody systems, Comput. Struct. 82 (17–19) (2004) 1359–1369. [29] I. Khemili, L. Romdhane, Dynamic analysis of a flexible slider-crank mechanism with clearance, Eur. J. Mech. A/Solids 27 (2008) 882–898. [30] S. Erkaya, I. Uzmay, Investigation on effect of joint clearance on dynamics of four-bar mechanism, Nonlinear Dyn. 58 (1–2) (2009) 179. [31] S. Erkaya, Experimental investigation of flexible connection and clearance joint effects on the vibration responses of mechanisms, Mech. Mach. Theory 121 (2018) 515–529. [32] J. Liu, B. Li, Theoretical and experimental identification of clearance nonlinearities for a continuum structure, J. Comput. Nonlinear Dyn. 11 (4) (2016) 041019. [33] O. Muvengei, J. Kihiu, B. Kua, Dynamic analysis of planar rigid-body mechanical systems with two-clearance revolute joints, Nonlinear Dyn. 73 (1–2) (2013) 259–273. [34] B. Li, R.Y. Chen, Z.J. He, Dynamic model and experimental research on pitch mechanism of gun with clearances, Adv. Mater. Res. 228 (2011) 125–129. [35] S.M. Varedi-Koulaei, H.M. Daniali, M. Farajtabar, Reducing the undesirable effects of joints clearance on the behavior of the planar 3-RRR parallel manipulators, Nonlinear Dyn. 86 (2) (2016) 1007–1022. [36] X. Jia, D. Jin, L. Ji, Investigation on the dynamic performance of the tripod-ball sliding joint with clearance in a crank–slider mechanism. Part 1: theoretical and experimental results, J. Sound Vib 252 (5) (2002) 919–933. [37] K. Soong, B.S. Thompson, A theoretical and experimental investigation of the dynamic response of a slider-crank mechanism with radial clearance in the gudgeon-pin joint, J. Mech. Des. 112 (2) (1990) 183–189. [38] R.S. Haines, A theory of contact loss at revolute joints with clearance, J. Mech. Eng. Sci. 22 (3) (1980) 129–136. [39] T.Y. Li, J.A. Yorke, Period three implies chaos, Am. Math. Mon. 82 (10) (1975) 985–992. [40] M.J. Feigenbaum, Quantitative universality for a class of nonlinear transformations, J. Stat. Phys. 19 (1) (1978) 25–52. [41] A. Wolf, J.B. Swift, H.L. Swinney, A.V. John, Determining Lyapunov exponents from a time series, Physica D 16 (3) (1985) 285–317. [42] P. Melby, J. Kaidel, N. Weber, Adaptation to the edge of chaos in the self-adjusting logistic map, Phys. Rev. Lett. 84 (26) (20 0 0) 5991. [43] Y.M. Zhang, X.K. Tang, Determination of slide postures of prismatic pair with clearance in spatial mechanisms, J. Mech. Sci. Technol. 17 (6) (1998) 889–892. [44] S.Z. Yan, P.F. Guo, Kinematic accuracy analysis of flexible mechanisms with uncertain link lengths and joint clearances., Proc. Inst. Mech. Eng. C –J. Mech. 225 (2011) 1973–1983. [45] F. Marques, P. Flores, J.C.P. Claro, H.M. Lankarani, Modeling and analysis of friction including rolling effects in multibody dynamics: a review, Multibody Syst. Dyn. 45 (2) (2019) 223–244. [46] K.L. Johnson, Contact Mechanics, Cambridge Univ. Press, 1987. [47] Q. Tian, C. Liu, M. Machado, P. Flores, A new model for dry and lubricated cylindrical joints with clearance in spatial flexible multibody systems, Nonlinear Dyn. 64 (1–2) (2011) 25–47. [48] M.J. Puttock, E.G. Thwaite, Elastic Compression of Spheres and Cylinders At Point and Line Contact, Commonwealth Scientific and Industrial Research Organization, Melbourne, Australia, 1969. [49] S. Yigit, A.G. Ulsoy, R.A. Scott, Spring-dashpot models for the dynamics of a radially rotating beam with impact, J. Sound Vib. 142 (3) (1990) 515–525. [50] X.P. Wang, G. Liu, S.X. Ma, Dynamic characteristics of mechanisms with revolute clearance joints, J. Vib. Shock 35 (7) (2016) 110–115. [51] X.Y. Geng, M. Li, Y.F. Liu, W. Zheng, Z.J. Zhao, Non-probabilistic kinematic reliability analysis of planar mechanisms with non-uniform revolute clearance joints, Mech. Mach. Theory 140 (2019) 413–433. [52] W.W.K. Xiang, S.Z. Yan, J.N. Wu, Dynamic analysis of planar mechanical systems considering stick-slip and Stribeck effect in revolute clearance joints, Nonlinear Dyn. 95 (1) (2019) 321–341. [53] L.D. Seneviratne, S.W.E. Earles, Chaotic behaviour exhibited during contact loss in a clearance joint of a four-bar mechanism, Mech. Mach. Theory 27 (3) (1992) 307–321. [54] A.Y. Loskutov, Dynamical chaos: systems of classical mechanics, Phys. Usp 50 (9) (2007) 939. [55] F.C. Moon, Experiments on chaotic motions of a forced nonlinear oscillator: strange attractors, J. Appl. Mech. 47 (3) (1980) 638–644. [56] M.T. Rosenstein, J.J. Collins, C.J.D. Luca, A practical method for calculating largest Lyapunov exponents from small data sets, Physica D 65 (1–2) (1993) 117–134.