# Numerical and Experimental Investigations on Dynamic Response of Hydraulic Cylinder with 3D Spatial Joints considering Radial and Axial Clearances

Authors: Ke Chen; Guojun Zhang; Haishuo Wang; Rui Wu; Hongmei Zheng; Shunhua Chen
Journal: Shock and Vibration
 (2019)

---

## Abstract

Radial clearance, particularly the axial clearance in the 3D joint of a mechanism owing to the assemblage, manufacturing tolerances, wear, and other conditions, has become a research focus in the field of multibody dynamics in recent years. In this study, a hydraulic cylinder model with 3D clearance joints was constructed by combining various potential contact scenarios. The novelty of this study is that potential contact points between the bearing wall and journal were calculated when the bearing wall circle was projected to an ellipse owing to misalignment of axes. Moreover, the simulation model considered the effective bulk modulus of the hydraulic oil and applied the Lagrange multiplier method. Subsequently, an experiment was conducted to verify the simulation results. The simulation and experimental results indicated that the dynamic responses of the hydraulic cylinder with 3D clearance joints can be classified as free, rebound, slide, and contact. The effects of input force, frequency, and clearance size on the dynamic behavior of the hydraulic cylinder were also investigated. Increasing the input force and clearance size will degrade the hydraulic cylinder dynamic response; however, the input force frequency can reduce the deterioration of the dynamic response. This study aids in providing improved understanding of the hydraulic cylinder with 3D clearances in the theoretical field and for practical engineering applications.


---

## Body

## 1. Introduction

The asymmetric hydraulic cylinder, which is used as an activator in construction equipment, agricultural machinery, special vehicles, and other heavy loading machinery, has attracted substantial research interest in recent years. Chen et al. analyzed the static and dynamic performance of a hydraulic cylinder with adaptive variable clearance around the piston, which can be applied in the high-frequency response and high-speed fields [1]. Ylinen et al. presented a linear hydraulic actuator model for multibody simulations, considering the state variables and the sealing friction, and found that the cylinder would soften the system dynamics in the starting and stopping states [2]. Moreover, Nizhegorodov et al. studied the frequency spectrum of the vibration exposure in the hydraulic cylinder [3]. Hydraulic cylinders experience impulsive shock and various other loading conditions in engineering machinery, such as breakers, drifters, or drills. The impulsive responses in mechanical systems with multiple clearances have been studied when the excitation and system load changed abruptly [4]. Ficarella et al. analyzed the impact energy of the breaker to minimize the required power to the hydraulic system using numerical and experimental data [5, 6]. Li et al. studied a hydraulic rock drill and investigated the overlapped reversing valve to propose a fault diagnosis method [7, 8].Owing to manufacturing and assembling errors, the actual joints obtain radial clearance. The phenomenon of clearance has received the attention of researchers since the 1970s [9]. Erkaya and Uzmay substituted clearance joints with a virtual massless link connecting the journal and bearing centers, which have been widely applied to estimate the kinematic error [10]. Furthermore, a series of scholars determined that the clearance can dramatically change the dynamic behavior of a mechanical system, and the interaction force between joint parts exhibits peaks [11, 12]. The point contact is used extensively and commonly referred to as the penalty approach [13], which allows the relative penetration of bodies to be applied to evaluate the reaction contact force. The equations for the dynamic analysis of constrained multibody systems are based on the Lagrange multiplier and Baumgarte stabilization method [14]. Several contact force models depending on different assumptions have been used to describe revolute joint models with clearance. In the nonlinear Hertz model [15], the contact forces between contacting bodies form an ellipse in the contact area, which is evaluated depending on the penetration depth. Moreover, a damping factor was introduced to remedy the disadvantage of energy dissipation during the impact process. The models of Lankarani and Nikravesh [16] and Hunt and Crossley [17] are the most commonly applied approaches for dealing with point contact situations. Dietl et al. analyzed the Hertz pressure distribution and replaced the nonlinear exponent of 1.5 in the Lankarani–Nikravesh model with 1.08, and the results were experimentally validated [18]. Koshy et al. used an experiment to illustrate that an appropriate contact force model with proper dissipative damping has a significant effect on the dynamic response of a mechanical system [19]. Wang et al. excluded the Lankarani–Nikravesh model’s restitution coefficient and replaced it with a damping coefficient [20]. The contact scenarios in spatial revolute joints with clearance appear as line contacts between cylindrical contact surfaces, and the shape of the line contact is rectangular. Different models have been developed for the relationships between the force and penetration. The model presented by Johnson [21] was based on the Hertz theory and subsequently accepted in the ESDU Tribology Series [22], which provided a comprehensive review covering the contact force and penetration expressions for engineering applications. Meanwhile, other contact force models are available to describe the contact between the two cylinders, including the Goldsmith [23], Radzimovsky [24], and Dubowsky–Freudenstein [25] models. Pereira et al. compared the precisions of the above line contact models and concluded that the Johnson and Radzimovsky models are suitable for describing the contact involving the cylinder length in most practical applications [26]. However, the method for solving the line contact force model using the Lambert function is excessively complicated [27]. Thus, it is necessary to determine a suitable method for solving the line contact model rapidly and precisely. The literature mentioned above mainly focused on clearance joint models. However, actual mechanisms have multiple joints, and thus, multiple joints with clearance make the dynamics of practical mechanisms more complex. It is therefore necessary that much attention should be paid to mechanisms with more than one clearance joint. Lu et al. [28, 29] and Wei et al. [30] analyzed the dynamics of a vehicle shimmy system model with multiple clearance joints. Moreover, several factors have been considered for decreasing the vibration dynamics owing to clearances.However, the abovementioned existing studies have only focused on the radial clearance, omitting the axial clearances existing in the joints to simplify 3D revolute joint models as planar models. The 3D revolute clearance joint models are more complex than planar ones, and additional degrees of freedom are incorporated into systems with spatial revolute joints [31]. Marques et al. applied a spatial slider-crank mechanism to study the effects of dry spatial revolute joints with radial and axial clearances [32]. Marques et al. accounted for all possible contact scenarios in the clearance joints and demonstrated that the clearance in the mechanical joints would have a significant impact on the dynamic response of the mechanical system [33]. Yang et al. illustrated that the multiple spatial revolute clearance joint would influence the dynamic behaviors of the open-loop mechanism to a greater extent than that of planar mechanisms, and the joints would interact with one another to cause fiercer vibrations, such as frequency, amplitude, and penetration [34]. Cavalieri and Cardona replaced the journal with two connecting spheres to construct a model of 3D revolute joints with clearance [35]. Ambrósio and Pombo proposed a unified formulation using the same kinematic information to combine perfect and clearance/bushing joints. This method enables similar modeling of two joint models and easy permutation of multibody systems [36]. Isaac et al. proposed a 3D finite element analysis approach to deal with nonideal spatial revolute joints [37].Several experiments have been conducted to verify the effects of clearance joints on the multibody dynamics. The surface wear on the axle and bushing of the revolute joint were investigated by Tasora et al. [38]. Flores et al. presented an experimental verification of a model of the slider-crank mechanism with revolute clearance joints [39]. Erkaya and Uzmay used three accelerometers and two microcameras to study the vibrations of a mechanism with two clearance joints [40]. Haroun and Megahed compared the results of a wear experiment and simulation, which proved the accuracy and efficiency of the simulation method [41]. Lai et al. validated the numerical wear prediction and experimental results of a planar mechanism with a revolute joint at a low velocity [42]. Brito et al. studied the operating effectiveness of an oil bearing with a twin groove or single groove [43]. The sensors for measuring the relative shaft position inside the bearing were arranged in the orthogonal directions [44–46], while vertical position arrangement was used to measure the displacements in the directions of the X- and Y-axes. Tian et al. [47] surveyed the models of planar and spatial multibody mechanical systems with clearance joints and conducted experiments on the mechanical systems with clearance joints. Thereafter, they discussed the phenomena commonly associated with clearance joint models, such as wear, nonsmooth behavior, optimization, chaos, and uncertainty [48]. Erkaya analyzed a flexible connection to minimize the clearance-induced dynamic deterioration by computational and experimental investigations [49, 50]. This is because the flexible connections have a clear suspension on the dynamics of a mechanism. In another study, Erkaya attempted to address the clearance problem in joints and improve the accuracy of robotic manipulators [51]. Yan et al. proposed a model for 3D revolute joints with radial and axial clearances and verified the model experimentally [27]. In the study, a comprehensive contact model was used to describe the contact scenarios, such as the point, line, and face contacts. The simulation results provided a realistic description of the 3D motion. However, several errors existed owing to the simplification of the potential contact points. The journal axle was misaligned with the bearing and eccentricity emerged. Therefore, the projection of the bearing wall circle on the lateral view was an ellipse, for which the angle error should be studied.Nevertheless, the effects of the revolute 3D joint with axial and radial clearances on the dynamic properties of hydraulic cylinders in practical engineering applications have rarely been reported. Under actual working conditions, the clearances existing at the joints interact with the hydraulic cylinders, which will deteriorate the dynamic characteristics of the cylinder. Therefore, understanding the interaction effect between hydraulic cylinders with 3D clearance joints has scientific and industrial significance. In this study, we focused on a hydraulic cylinder with a 3D joint with axial and radial clearances considering the contact force model and hydraulic cylinder model in Section2. The novelty of this study is that potential contact points between the bearing wall and journal were calculated when the bearing wall circle was projected to an ellipse. Moreover, we simplified the line contact force by using the Newton iteration and not the Lambert function as described in a previous study [27]. Numerical and experimental examples of the hydraulic cylinder model with the 3D clearance joint and its parameters are described in Section 3. Subsequently, the dynamic responses of the simulation and experiment are compared. The effects of different inputting force frequencies and clearance sizes were studied, and the results are presented in Section 4. Finally, several conclusions of the study are provided in Section 5.


## 2. Dynamic Model of Hydraulic Cylinder System with Clearance

The top view of the hydraulic cylinder, including a cylinder barrel, piston, rod, forward camber, and return camber, is presented in Figure1. On the bilateral side of the hydraulic cylinder, there is one spatial revolute joint, each of which has a journal, bearing, and flange.Figure 1
Model of the hydraulic cylinder with two clearances on bilateral sides.


### 2.1. Model of 3D Revolute Joints with Radial and Axial Clearances

In this section, the geometric model for the 3D revolute joints with radial and axial clearances is presented. The 3D revolute joints contain a journal and bearing, as illustrated in Figure2. Two flanges exist on the bilateral side of the journal, which limit the range of the bearing axial motion.Figure 2
Structure of 3D revolute joint showing clearances.The bearing wall is on the inner side of the bearing, and the cylindrical surface of the journal faces the corresponding bearing wall. The radius of the bearing wall is greater than that of the journal surface in the radial direction. Therefore, a radial clearancecr is generated. In the axial direction, the bearing length is smaller than the distance between the two flanges connected to the journal. Thus, an axial clearance ca is created.The coupling interactions of the radial and axial clearances create complex contact scenarios, as illustrated in Figure3. These scenarios can be classified as four different types of journal motion inside the bearing [33, 52]:(1)
Free flight motion, where no contact occurs between the two elements, as illustrated in Figure3(a)(2)
The journal surface contacts the bearing wall at a line, as indicated in Figure3(e)(3)
The contact area is between the journal flange and bearing base, as illustrated in Figure3(f)(4)
The journal contacts the bearing wall at certain points: one point in Figures3(b) and 3(c), two points in Figures 3(d), 3(j), and 3(k), three points in Figures 3(i) and 3(l), and four points in Figure 3(h)Figure 3
Potential contact scenarios in 3D revolute joint with clearances.
(a)
(b)
(c)
(d)
(e)
(f)
(g)
(h)
(i)
(j)
(k)
(l)Figure4 illustrates the angle error θerr between points A and B when the eccentricity of the outer ellipse is 0.6, which is obvious and should not be neglected. Thus, it is necessary to study the eccentricity effect of the ellipse.Figure 4
Schematic overview of angle error when the bearing ellipse contacts the journal circle.For the description of the 3D joint model with clearances, several coordinate systems are indicated on parts of the joints in Figure5. Firstly, the global coordinate system Og-xyz is fixed on the ground and provides a basic reference system for the entire model, as illustrated in Figure 5. Secondly, the local coordinate system Oi-xiyizi is positioned at the mass center of the part connected to the bearing, and the matrix Agi describes the transformation from the global coordinate system to the local coordinate system of the bearing. Finally, the local coordinate system Oj-xjyjzj is located at the mass center of the part connected to the journal. The transformation matrix from the ground to the journal is defined by Agj. Here, Agi and Agj are Euler angle transformation matrices, in which α, β, and γ are Euler angle parameters. In this study, the translation sequence was 313. Note that there are two points Oa and Ob that can easily be used to describe the position of the journal and bearing. Point Oa is the center of one ring of the bearing wall and Ob is the journal center. Meanwhile, there exist local coordinate systems on points Oa and Ob. The axes of these local coordinate systems Oa and Ob are parallel to the coordinate systems Oi-xiyizi and Oj-xjyjzj, respectively.Figure 5
Coordinate systems in the 3D revolute joint with clearances.Because the journal is mounted on the slider of the vibration bench, only theY-direction movement exhibits freedom. Thus, the local coordinate system Oj-xjyjzj provides a simplified description of the 3D joint model. The transformation matrix Aij represents the transformation from the local coordinate system Oj-xjyjzj to Oi-xiyizi, and can be expressed as(1)Aij=AgjAig.The pointsQ in the global coordinate system Og-xyz can be described as follows:(2)rgQ=rgOi+AigriOiQ,where rgOi is the vector from the origin point Og to the circle center point Oi in the global coordinate system, which is subscripted as g, while riOiQ is the vector connecting the circle center Oi to any point Q in the journal local coordinate system Oj-xjyjzj. The point Q can be described as(3)rjQ=AgjrgOjOi+AgjAigriOiQ,where the vector rgOjOi is expressed as(4)rgOjOi=rgOi−rgOj.In the lateral view, the circle of the bearing is tilted and the shape is changed into an ellipse. The major radiusa and minor radius b of the ellipse can be calculated as follows:(5)a=Ri,b=Ri⋅cosβij,where Ri is the radius of the inner bearing wall and the angle βij is the Euler angle in the translation Aij.In (5), the pointQ on the ellipse boundary can be illustrated in the coordinate system Oj-xjyjzj:(6)xjQ=xjOa+a⋅cosθ+αij,yjQ=yjOa+a⋅cosθ⋅sinβij,zjQ=zjOa+b⋅sinθ+αij.The pointOa is the center of one ring of the bearing wall, which is depicted in the same manner as Q. The initial angle αij is the Euler angle in the translation Aji. Moreover, xjOa, yjOa, and zjOa are the ellipse centers in the journal original system.The distance fromQ to the journal center can be expressed as(7)d=xjQ−xjOb2+zjQ−zjOb2,where Ob is the journal center and xjOb, yjOb, and zjOb are the coordinate values of the journal center.The penetration depthδr between the journal and bearing in Figure 6(a) can be determined as follows:(8)δr=Rj−d,where d is the distance between the center Ob and Qmin. When the penetration depth δr is shorter than zero, a potential collision may occur. The contact conditions can be summarized as follows:δr>0: free movement between the journal and bearingδr≤0: contact-impact between the journal and bearingFigure 6
Schematic of 3D joint geometry inQmin point calculation process.
(a)
(b)Q
0 is located at the intersection point of the line and bearing wall. However, the point Qmin is closest to the center Ob of the journal. The following part demonstrates the calculation process for determining point Qmin to cause d to be minimal, for which the gradient descent method is applied. The target function of the gradient descent method is defined as(9)fθ=a⋅cosθ+αji−xjOb2+b⋅sinθ+αji−zjOb2.In this function, the coordinate valuesxjOb and yjOb of the journal center Ob are constant. Then, the derivative equation is calculated as follows:(10)f′θ=2a⋅cosθ+αji−xjOb⋅asinθ−2b⋅sinθ+αji−zjOb⋅bcosθ.The flowchart for the gradient descent method is presented in Figure7. The angle θ is set as the following initial value:(11)θinit=arctanzjOaObxjOaOb.Figure 7
Flowchart of the gradient descent method.This is becauseQmin is close to the line connecting Oa to Ob, which can enable the calculation to converge rapidly. In the gradient descent method, the maximum iteration step is 7 and the learning rate α = 0.0005 for determining θmin. The point Qmin coordinate can be calculated as follows:(12)xjQmin=xjOa+a⋅cosθmin+αij,yjQmin=yjOa+a⋅cosθmin⋅sinβij,zjQmin=zjOa+b⋅sinθmin+αij.The contact pointP in the journal can be described as(13)xjP=xjOb+Rj⋅cosθmin+αij,yjP=yjOb+Rj⋅cosθmin⋅sinβij,zjP=zjOb+Rj⋅sinθmin+αij⋅cosβji.Another contact pointQa exists between the flange and bearing base in the axial direction, as shown in Figure 8. The potential point is at the acme of the ellipse and can be expressed as(14)xjQa=xjOa+a⋅cosπ+αji,yjQa=yjOa+a⋅cosθ⋅sinβji,zjQa=zjOa+b⋅sinπ+αji.Figure 8
Schematic view of 3D joint in axial contact.The penetration depthδa is calculated by(15)δa=yjQa−yjflange.Moreover, the axial contact scenarios can be expressed as follows:δa>0: axial free movement between the journal and bearingδa≤0: axial contact-impact between the journal and bearingThe radial relative velocityvr of the contact points is required to calculate the contact and friction forces in the local coordinate system Oj-xjyjzj, which is described as(16)vr=r˙jP−r˙jQmin,where r˙jP and r˙jQmin are the velocities of contact points P and Qmin from previous positions at a simulation step before, respectively, in the local coordinate system Oj-xjyjzj. The velocity vr can be projected as vrT onto a plane that is tangential to the journal through the collision point Qmin, as illustrated in Figure 9. The relative normal velocity vrN is projected onto the normal direction of the above plane.Figure 9
Relative velocity in 3D joint with clearances.Alternatively, during the occurrence of axial contact, the axial relative velocityva can be expressed as(17)va=r˙jPa−r˙jQa,where r˙jPa and r˙jQa are the velocities of the axial impact points Pa and Qa at a simulation step, respectively, in the coordinate system Oj-xjyjzj. The velocity va is projected as vaT onto the plane of the flange and vaN in the normal orientation of the flange plane.Here,δ˙k is the velocity of the relative penetration δk and can be calculated as(18)δ˙k=vkN,k=r,a.The direction of the relative normal penetration is calculated as(19)ekN=rjPk−rjQkδk,k=r,a.The magnitude of the relative tangential velocityvT can be obtained as follows:(20)vkT=vkT,k=r,a.The direction of the relative tangential velocity is described by(21)ekT=vkTvkT,k=r,a,where ekN and ekT are all in the local coordinate system Oj-xjyjzj.The contact force model can be decomposed into the normal and tangential contact forcesFN and FT, respectively. Sections 2.1.1 to 2.1.3 demonstrate the calculation of the normal contact forces FN, while Section 2.1.4 presents the solving process for FT.


#### 2.1.1. Point Contact Model

The point contactFN is evaluated using the Lankarani–Nikravesh contact force model:(22)FN=Kδrn1+31−ce24δ˙rδ˙r−,where ce is the restitution coefficient, the exponent n is usually 1.5 for a collision between metal bodies, and δ and δ˙− are the initial velocities when the impact occurs. The stiffness K can be evaluated as follows:(23)K=4π3E∗RiRjRi−Rj1/2.The value ofE∗ can be determined by(24)E∗=1−vi2Ei+1−vj2Ej,where νi, Ei and νj, Ej indicate Poisson’s ratio and Young’s modulus of each collision body.


#### 2.1.2. Line Contact Model

The line contact area between the journal and bearing wall is a rectangular place, which appears between the cylindrical contact surfaces. The distributed loadP per unit length is expressed as P=F/l. In the Johnson model, the relative penetration δl is defined as(25)δl=PπE∗ln4πE∗ΔRP−1,where the parameter ΔR=Ri−Rj represents the difference between the cylinder radii. To guarantee that, equation (25) satisfies the physical requirements, and the following function must be greater than 1:(26)ln4πE∗ΔRP−1>1.Thus,P must follow the following condition:(27)P≤4πE∗ΔRe2,where e represents natural exponential. In the domain (0, 4πE∗ΔR/e2), the Johnson function linking P to δl increases monotonically. The calculations involved for the line force model using the Lambert function are excessively complicated [27]. Therefore, P has only one value in this domain. We apply the Newton iteration to determine P at a certain penetration δl. The target function is defined as(28)fP=PπE∗ln4πE∗ΔRP−1−δl.The derivative is obtained as follows:(29)f′P=1πE∗ln4πE∗ΔRP−1−1πE∗.The iteration algorithm can be applied iteratively to calculate(30)Pn+1=Pn−fPnf′Pn.The iteration error is set to less than 0.1 N, and the maximum iteration step is 7. Then, the final value of the loadP is obtained. The line contact force is expressed as [37](31)FN=Pl+Dδ˙r,where D is the damping coefficient.


#### 2.1.3. Area Contact Model

When the journal and bearing axles are parallel to one another and the bearing penetrates the journal flange, area contact will occur between the journal flange and bearing base. The contact area is a ring and can be represented by(32)A=πRbase2−πRi2,where Rbase is the radius of the bearing base illustrated in Figure 3.The contact force is denoted by [44](33)FN=πE∗δaA.The normal impact forceFN is obtained as follows:(34)FN=FN⋅eN,where FN represents the point, line, and area contact forces mentioned previously and eN is the unit vector in equation (19).


#### 2.1.4. Friction Model

When relative sliding occurs, the tangential friction forceFT affects the journal and bearing to prevent sliding, which can be described as [42](35)FT=0,vT≤v0,−cfvT−v0v1−v0FN,v0<vT<v1,−cfFN,vT≥v1,where cf is the friction coefficient, vT is the relative tangential velocity, and v0 and v1 are the given tolerances for the tangential velocity, which are selected as v0=0.001mm/s and v1=0.01mm/s in this case.The tangential friction forceFT is obtained as follows:(36)FT=FT⋅eT,where eT is the unit vector in equation (24).As the contact forces between the journal and bearing have been calculated in different scenarios, it is necessary to combine the interaction forces into the dynamic system. Thus, the forces acting on the journal can be expressed as(37)Fj=AjgFT+FN,Mj=Fj⋅rgQM,where rgQMk is the vector from the potential collision impact point Q to the mass center M, M1 and M2 are the mass centers of the cylinder barrel and rod, as illustrated in Figures 3 and 10, respectively, and M3 is the combination of the slider and moving support:(38)Fi=−Fj,Mi=Fi⋅rgPM,where rgPM is the vector connecting the potential collision impact point P on the bearing to the mass center M.Figure 10
Hydraulic cylinder and forces model.


### 2.2. Hydraulic Cylinder Model

The rod moves in the cylinder barrel, causing the volumes of the forward and return cambers to change. When leakage in the hydraulic cylinder is neglected, the effective bulk modulusβe of the hydraulic oil causes pressure fluctuation with the change in volume [53]. The pressure p1 of the forward camber can be calculated as(39)p1=p10−βeA1zrodV10+A1zrod,where A1 is the area of the piston in the forward camber, V10 is the initial volume of the forward camber, p10 is the initial pressure of the forward camber, and zrod is the relative displacement of the piston from its initial position.The pressurep2 of the return camber can be determined in the same manner:(40)p2=p20+βeA2zrodV20+A2zrod,where A2 is the return tank flow, V20 is the initial volume of the return camber, and p20 is the initial pressure of the return camber.The force on the piston and rod is described by Newton’s equation [54]:(41)p1A1−p2A2−f=m2z¨rod−m2gsinθ+FR,where z¨rod is the piston acceleration; m2 is the rod mass; θ is the angle between the axis and the horizon, as shown in Figure 10; and FR is the decomposition of the impact force on the rod journal along the axial direction of the hydraulic cylinder. The impact force is evaluated as follows:(42)FK=AkgFkez,k=L,R,where AKg represents the transformation from the local coordinate system of the journal to the global coordinate system. ez is the Z-axle unit vector in the local coordinate system and is equal to [0, 0, 1]T. Moreover, Fk is the interaction force in the left or right joints when contact-impact occurs. The value f in equation (41) denotes the linear friction with stiction between the rod and cylinder barrel, as described in the following equation [55]:(43)f=FS⋅sgnvT,vT≤vS,FS−vT−vSvC−vS⋅FS−FCsgnvT,vS<vT≤vC,FC·sgnvT,vT>vC,where FS is the static friction and FC is the dynamic friction. In this study, FS = 385 N and FC = 340 N. Furthermore, vT is the rod velocity. The tolerances for the velocity in equation (43) were selected as vS = 0.01 mm/s and vC = 0.05 mm/s.


### 2.3. Dynamic Constraint Equations

The kinematical constraint equation in a multibody system can be expressed as(44)Φqr,t=0,where qr indicates the position and orientation vectors for describing the bodies in global coordinates and t represents the time. The motivation equations for a constrained multibody system in Cartesian coordinates can be represented as follows:(45)Mq¨+ΦqTλ=g,where M is the system mass matrix and the components of M denote the masses and moments of inertia of the bodies in the system. Moreover, ΦqT represents the transformation matrix, λ is the vector of Lagrange multipliers, and g represents the generalized forces, which include the external driving torques or forces, centrifugal force, Coriolis force, and impact forces of the clearance joints:(46)Φqq¨=−Φqq˙qq˙−2Φqtq˙−Φtt=γ.Equation (45) is combined with equation (44) and rewritten as the matrix of differential algebraic equations that describes the motion of a multibody system:(47)MΦqTΦq0q¨rλ=gγ.


## 2.1. Model of 3D Revolute Joints with Radial and Axial Clearances

In this section, the geometric model for the 3D revolute joints with radial and axial clearances is presented. The 3D revolute joints contain a journal and bearing, as illustrated in Figure2. Two flanges exist on the bilateral side of the journal, which limit the range of the bearing axial motion.Figure 2
Structure of 3D revolute joint showing clearances.The bearing wall is on the inner side of the bearing, and the cylindrical surface of the journal faces the corresponding bearing wall. The radius of the bearing wall is greater than that of the journal surface in the radial direction. Therefore, a radial clearancecr is generated. In the axial direction, the bearing length is smaller than the distance between the two flanges connected to the journal. Thus, an axial clearance ca is created.The coupling interactions of the radial and axial clearances create complex contact scenarios, as illustrated in Figure3. These scenarios can be classified as four different types of journal motion inside the bearing [33, 52]:(1)
Free flight motion, where no contact occurs between the two elements, as illustrated in Figure3(a)(2)
The journal surface contacts the bearing wall at a line, as indicated in Figure3(e)(3)
The contact area is between the journal flange and bearing base, as illustrated in Figure3(f)(4)
The journal contacts the bearing wall at certain points: one point in Figures3(b) and 3(c), two points in Figures 3(d), 3(j), and 3(k), three points in Figures 3(i) and 3(l), and four points in Figure 3(h)Figure 3
Potential contact scenarios in 3D revolute joint with clearances.
(a)
(b)
(c)
(d)
(e)
(f)
(g)
(h)
(i)
(j)
(k)
(l)Figure4 illustrates the angle error θerr between points A and B when the eccentricity of the outer ellipse is 0.6, which is obvious and should not be neglected. Thus, it is necessary to study the eccentricity effect of the ellipse.Figure 4
Schematic overview of angle error when the bearing ellipse contacts the journal circle.For the description of the 3D joint model with clearances, several coordinate systems are indicated on parts of the joints in Figure5. Firstly, the global coordinate system Og-xyz is fixed on the ground and provides a basic reference system for the entire model, as illustrated in Figure 5. Secondly, the local coordinate system Oi-xiyizi is positioned at the mass center of the part connected to the bearing, and the matrix Agi describes the transformation from the global coordinate system to the local coordinate system of the bearing. Finally, the local coordinate system Oj-xjyjzj is located at the mass center of the part connected to the journal. The transformation matrix from the ground to the journal is defined by Agj. Here, Agi and Agj are Euler angle transformation matrices, in which α, β, and γ are Euler angle parameters. In this study, the translation sequence was 313. Note that there are two points Oa and Ob that can easily be used to describe the position of the journal and bearing. Point Oa is the center of one ring of the bearing wall and Ob is the journal center. Meanwhile, there exist local coordinate systems on points Oa and Ob. The axes of these local coordinate systems Oa and Ob are parallel to the coordinate systems Oi-xiyizi and Oj-xjyjzj, respectively.Figure 5
Coordinate systems in the 3D revolute joint with clearances.Because the journal is mounted on the slider of the vibration bench, only theY-direction movement exhibits freedom. Thus, the local coordinate system Oj-xjyjzj provides a simplified description of the 3D joint model. The transformation matrix Aij represents the transformation from the local coordinate system Oj-xjyjzj to Oi-xiyizi, and can be expressed as(1)Aij=AgjAig.The pointsQ in the global coordinate system Og-xyz can be described as follows:(2)rgQ=rgOi+AigriOiQ,where rgOi is the vector from the origin point Og to the circle center point Oi in the global coordinate system, which is subscripted as g, while riOiQ is the vector connecting the circle center Oi to any point Q in the journal local coordinate system Oj-xjyjzj. The point Q can be described as(3)rjQ=AgjrgOjOi+AgjAigriOiQ,where the vector rgOjOi is expressed as(4)rgOjOi=rgOi−rgOj.In the lateral view, the circle of the bearing is tilted and the shape is changed into an ellipse. The major radiusa and minor radius b of the ellipse can be calculated as follows:(5)a=Ri,b=Ri⋅cosβij,where Ri is the radius of the inner bearing wall and the angle βij is the Euler angle in the translation Aij.In (5), the pointQ on the ellipse boundary can be illustrated in the coordinate system Oj-xjyjzj:(6)xjQ=xjOa+a⋅cosθ+αij,yjQ=yjOa+a⋅cosθ⋅sinβij,zjQ=zjOa+b⋅sinθ+αij.The pointOa is the center of one ring of the bearing wall, which is depicted in the same manner as Q. The initial angle αij is the Euler angle in the translation Aji. Moreover, xjOa, yjOa, and zjOa are the ellipse centers in the journal original system.The distance fromQ to the journal center can be expressed as(7)d=xjQ−xjOb2+zjQ−zjOb2,where Ob is the journal center and xjOb, yjOb, and zjOb are the coordinate values of the journal center.The penetration depthδr between the journal and bearing in Figure 6(a) can be determined as follows:(8)δr=Rj−d,where d is the distance between the center Ob and Qmin. When the penetration depth δr is shorter than zero, a potential collision may occur. The contact conditions can be summarized as follows:δr>0: free movement between the journal and bearingδr≤0: contact-impact between the journal and bearingFigure 6
Schematic of 3D joint geometry inQmin point calculation process.
(a)
(b)Q
0 is located at the intersection point of the line and bearing wall. However, the point Qmin is closest to the center Ob of the journal. The following part demonstrates the calculation process for determining point Qmin to cause d to be minimal, for which the gradient descent method is applied. The target function of the gradient descent method is defined as(9)fθ=a⋅cosθ+αji−xjOb2+b⋅sinθ+αji−zjOb2.In this function, the coordinate valuesxjOb and yjOb of the journal center Ob are constant. Then, the derivative equation is calculated as follows:(10)f′θ=2a⋅cosθ+αji−xjOb⋅asinθ−2b⋅sinθ+αji−zjOb⋅bcosθ.The flowchart for the gradient descent method is presented in Figure7. The angle θ is set as the following initial value:(11)θinit=arctanzjOaObxjOaOb.Figure 7
Flowchart of the gradient descent method.This is becauseQmin is close to the line connecting Oa to Ob, which can enable the calculation to converge rapidly. In the gradient descent method, the maximum iteration step is 7 and the learning rate α = 0.0005 for determining θmin. The point Qmin coordinate can be calculated as follows:(12)xjQmin=xjOa+a⋅cosθmin+αij,yjQmin=yjOa+a⋅cosθmin⋅sinβij,zjQmin=zjOa+b⋅sinθmin+αij.The contact pointP in the journal can be described as(13)xjP=xjOb+Rj⋅cosθmin+αij,yjP=yjOb+Rj⋅cosθmin⋅sinβij,zjP=zjOb+Rj⋅sinθmin+αij⋅cosβji.Another contact pointQa exists between the flange and bearing base in the axial direction, as shown in Figure 8. The potential point is at the acme of the ellipse and can be expressed as(14)xjQa=xjOa+a⋅cosπ+αji,yjQa=yjOa+a⋅cosθ⋅sinβji,zjQa=zjOa+b⋅sinπ+αji.Figure 8
Schematic view of 3D joint in axial contact.The penetration depthδa is calculated by(15)δa=yjQa−yjflange.Moreover, the axial contact scenarios can be expressed as follows:δa>0: axial free movement between the journal and bearingδa≤0: axial contact-impact between the journal and bearingThe radial relative velocityvr of the contact points is required to calculate the contact and friction forces in the local coordinate system Oj-xjyjzj, which is described as(16)vr=r˙jP−r˙jQmin,where r˙jP and r˙jQmin are the velocities of contact points P and Qmin from previous positions at a simulation step before, respectively, in the local coordinate system Oj-xjyjzj. The velocity vr can be projected as vrT onto a plane that is tangential to the journal through the collision point Qmin, as illustrated in Figure 9. The relative normal velocity vrN is projected onto the normal direction of the above plane.Figure 9
Relative velocity in 3D joint with clearances.Alternatively, during the occurrence of axial contact, the axial relative velocityva can be expressed as(17)va=r˙jPa−r˙jQa,where r˙jPa and r˙jQa are the velocities of the axial impact points Pa and Qa at a simulation step, respectively, in the coordinate system Oj-xjyjzj. The velocity va is projected as vaT onto the plane of the flange and vaN in the normal orientation of the flange plane.Here,δ˙k is the velocity of the relative penetration δk and can be calculated as(18)δ˙k=vkN,k=r,a.The direction of the relative normal penetration is calculated as(19)ekN=rjPk−rjQkδk,k=r,a.The magnitude of the relative tangential velocityvT can be obtained as follows:(20)vkT=vkT,k=r,a.The direction of the relative tangential velocity is described by(21)ekT=vkTvkT,k=r,a,where ekN and ekT are all in the local coordinate system Oj-xjyjzj.The contact force model can be decomposed into the normal and tangential contact forcesFN and FT, respectively. Sections 2.1.1 to 2.1.3 demonstrate the calculation of the normal contact forces FN, while Section 2.1.4 presents the solving process for FT.


### 2.1.1. Point Contact Model

The point contactFN is evaluated using the Lankarani–Nikravesh contact force model:(22)FN=Kδrn1+31−ce24δ˙rδ˙r−,where ce is the restitution coefficient, the exponent n is usually 1.5 for a collision between metal bodies, and δ and δ˙− are the initial velocities when the impact occurs. The stiffness K can be evaluated as follows:(23)K=4π3E∗RiRjRi−Rj1/2.The value ofE∗ can be determined by(24)E∗=1−vi2Ei+1−vj2Ej,where νi, Ei and νj, Ej indicate Poisson’s ratio and Young’s modulus of each collision body.


### 2.1.2. Line Contact Model

The line contact area between the journal and bearing wall is a rectangular place, which appears between the cylindrical contact surfaces. The distributed loadP per unit length is expressed as P=F/l. In the Johnson model, the relative penetration δl is defined as(25)δl=PπE∗ln4πE∗ΔRP−1,where the parameter ΔR=Ri−Rj represents the difference between the cylinder radii. To guarantee that, equation (25) satisfies the physical requirements, and the following function must be greater than 1:(26)ln4πE∗ΔRP−1>1.Thus,P must follow the following condition:(27)P≤4πE∗ΔRe2,where e represents natural exponential. In the domain (0, 4πE∗ΔR/e2), the Johnson function linking P to δl increases monotonically. The calculations involved for the line force model using the Lambert function are excessively complicated [27]. Therefore, P has only one value in this domain. We apply the Newton iteration to determine P at a certain penetration δl. The target function is defined as(28)fP=PπE∗ln4πE∗ΔRP−1−δl.The derivative is obtained as follows:(29)f′P=1πE∗ln4πE∗ΔRP−1−1πE∗.The iteration algorithm can be applied iteratively to calculate(30)Pn+1=Pn−fPnf′Pn.The iteration error is set to less than 0.1 N, and the maximum iteration step is 7. Then, the final value of the loadP is obtained. The line contact force is expressed as [37](31)FN=Pl+Dδ˙r,where D is the damping coefficient.


### 2.1.3. Area Contact Model

When the journal and bearing axles are parallel to one another and the bearing penetrates the journal flange, area contact will occur between the journal flange and bearing base. The contact area is a ring and can be represented by(32)A=πRbase2−πRi2,where Rbase is the radius of the bearing base illustrated in Figure 3.The contact force is denoted by [44](33)FN=πE∗δaA.The normal impact forceFN is obtained as follows:(34)FN=FN⋅eN,where FN represents the point, line, and area contact forces mentioned previously and eN is the unit vector in equation (19).


### 2.1.4. Friction Model

When relative sliding occurs, the tangential friction forceFT affects the journal and bearing to prevent sliding, which can be described as [42](35)FT=0,vT≤v0,−cfvT−v0v1−v0FN,v0<vT<v1,−cfFN,vT≥v1,where cf is the friction coefficient, vT is the relative tangential velocity, and v0 and v1 are the given tolerances for the tangential velocity, which are selected as v0=0.001mm/s and v1=0.01mm/s in this case.The tangential friction forceFT is obtained as follows:(36)FT=FT⋅eT,where eT is the unit vector in equation (24).As the contact forces between the journal and bearing have been calculated in different scenarios, it is necessary to combine the interaction forces into the dynamic system. Thus, the forces acting on the journal can be expressed as(37)Fj=AjgFT+FN,Mj=Fj⋅rgQM,where rgQMk is the vector from the potential collision impact point Q to the mass center M, M1 and M2 are the mass centers of the cylinder barrel and rod, as illustrated in Figures 3 and 10, respectively, and M3 is the combination of the slider and moving support:(38)Fi=−Fj,Mi=Fi⋅rgPM,where rgPM is the vector connecting the potential collision impact point P on the bearing to the mass center M.Figure 10
Hydraulic cylinder and forces model.


## 2.1.1. Point Contact Model

The point contactFN is evaluated using the Lankarani–Nikravesh contact force model:(22)FN=Kδrn1+31−ce24δ˙rδ˙r−,where ce is the restitution coefficient, the exponent n is usually 1.5 for a collision between metal bodies, and δ and δ˙− are the initial velocities when the impact occurs. The stiffness K can be evaluated as follows:(23)K=4π3E∗RiRjRi−Rj1/2.The value ofE∗ can be determined by(24)E∗=1−vi2Ei+1−vj2Ej,where νi, Ei and νj, Ej indicate Poisson’s ratio and Young’s modulus of each collision body.


## 2.1.2. Line Contact Model

The line contact area between the journal and bearing wall is a rectangular place, which appears between the cylindrical contact surfaces. The distributed loadP per unit length is expressed as P=F/l. In the Johnson model, the relative penetration δl is defined as(25)δl=PπE∗ln4πE∗ΔRP−1,where the parameter ΔR=Ri−Rj represents the difference between the cylinder radii. To guarantee that, equation (25) satisfies the physical requirements, and the following function must be greater than 1:(26)ln4πE∗ΔRP−1>1.Thus,P must follow the following condition:(27)P≤4πE∗ΔRe2,where e represents natural exponential. In the domain (0, 4πE∗ΔR/e2), the Johnson function linking P to δl increases monotonically. The calculations involved for the line force model using the Lambert function are excessively complicated [27]. Therefore, P has only one value in this domain. We apply the Newton iteration to determine P at a certain penetration δl. The target function is defined as(28)fP=PπE∗ln4πE∗ΔRP−1−δl.The derivative is obtained as follows:(29)f′P=1πE∗ln4πE∗ΔRP−1−1πE∗.The iteration algorithm can be applied iteratively to calculate(30)Pn+1=Pn−fPnf′Pn.The iteration error is set to less than 0.1 N, and the maximum iteration step is 7. Then, the final value of the loadP is obtained. The line contact force is expressed as [37](31)FN=Pl+Dδ˙r,where D is the damping coefficient.


## 2.1.3. Area Contact Model

When the journal and bearing axles are parallel to one another and the bearing penetrates the journal flange, area contact will occur between the journal flange and bearing base. The contact area is a ring and can be represented by(32)A=πRbase2−πRi2,where Rbase is the radius of the bearing base illustrated in Figure 3.The contact force is denoted by [44](33)FN=πE∗δaA.The normal impact forceFN is obtained as follows:(34)FN=FN⋅eN,where FN represents the point, line, and area contact forces mentioned previously and eN is the unit vector in equation (19).


## 2.1.4. Friction Model

When relative sliding occurs, the tangential friction forceFT affects the journal and bearing to prevent sliding, which can be described as [42](35)FT=0,vT≤v0,−cfvT−v0v1−v0FN,v0<vT<v1,−cfFN,vT≥v1,where cf is the friction coefficient, vT is the relative tangential velocity, and v0 and v1 are the given tolerances for the tangential velocity, which are selected as v0=0.001mm/s and v1=0.01mm/s in this case.The tangential friction forceFT is obtained as follows:(36)FT=FT⋅eT,where eT is the unit vector in equation (24).As the contact forces between the journal and bearing have been calculated in different scenarios, it is necessary to combine the interaction forces into the dynamic system. Thus, the forces acting on the journal can be expressed as(37)Fj=AjgFT+FN,Mj=Fj⋅rgQM,where rgQMk is the vector from the potential collision impact point Q to the mass center M, M1 and M2 are the mass centers of the cylinder barrel and rod, as illustrated in Figures 3 and 10, respectively, and M3 is the combination of the slider and moving support:(38)Fi=−Fj,Mi=Fi⋅rgPM,where rgPM is the vector connecting the potential collision impact point P on the bearing to the mass center M.Figure 10
Hydraulic cylinder and forces model.


## 2.2. Hydraulic Cylinder Model

The rod moves in the cylinder barrel, causing the volumes of the forward and return cambers to change. When leakage in the hydraulic cylinder is neglected, the effective bulk modulusβe of the hydraulic oil causes pressure fluctuation with the change in volume [53]. The pressure p1 of the forward camber can be calculated as(39)p1=p10−βeA1zrodV10+A1zrod,where A1 is the area of the piston in the forward camber, V10 is the initial volume of the forward camber, p10 is the initial pressure of the forward camber, and zrod is the relative displacement of the piston from its initial position.The pressurep2 of the return camber can be determined in the same manner:(40)p2=p20+βeA2zrodV20+A2zrod,where A2 is the return tank flow, V20 is the initial volume of the return camber, and p20 is the initial pressure of the return camber.The force on the piston and rod is described by Newton’s equation [54]:(41)p1A1−p2A2−f=m2z¨rod−m2gsinθ+FR,where z¨rod is the piston acceleration; m2 is the rod mass; θ is the angle between the axis and the horizon, as shown in Figure 10; and FR is the decomposition of the impact force on the rod journal along the axial direction of the hydraulic cylinder. The impact force is evaluated as follows:(42)FK=AkgFkez,k=L,R,where AKg represents the transformation from the local coordinate system of the journal to the global coordinate system. ez is the Z-axle unit vector in the local coordinate system and is equal to [0, 0, 1]T. Moreover, Fk is the interaction force in the left or right joints when contact-impact occurs. The value f in equation (41) denotes the linear friction with stiction between the rod and cylinder barrel, as described in the following equation [55]:(43)f=FS⋅sgnvT,vT≤vS,FS−vT−vSvC−vS⋅FS−FCsgnvT,vS<vT≤vC,FC·sgnvT,vT>vC,where FS is the static friction and FC is the dynamic friction. In this study, FS = 385 N and FC = 340 N. Furthermore, vT is the rod velocity. The tolerances for the velocity in equation (43) were selected as vS = 0.01 mm/s and vC = 0.05 mm/s.


## 2.3. Dynamic Constraint Equations

The kinematical constraint equation in a multibody system can be expressed as(44)Φqr,t=0,where qr indicates the position and orientation vectors for describing the bodies in global coordinates and t represents the time. The motivation equations for a constrained multibody system in Cartesian coordinates can be represented as follows:(45)Mq¨+ΦqTλ=g,where M is the system mass matrix and the components of M denote the masses and moments of inertia of the bodies in the system. Moreover, ΦqT represents the transformation matrix, λ is the vector of Lagrange multipliers, and g represents the generalized forces, which include the external driving torques or forces, centrifugal force, Coriolis force, and impact forces of the clearance joints:(46)Φqq¨=−Φqq˙qq˙−2Φqtq˙−Φtt=γ.Equation (45) is combined with equation (44) and rewritten as the matrix of differential algebraic equations that describes the motion of a multibody system:(47)MΦqTΦq0q¨rλ=gγ.


## 3. Numerical and Experimental Examples

### 3.1. Numerical Example Schema

A schematic overview of the numerical example is presented in Figure11. The hydraulic cylinder is inclined at 5° in the anticlockwise direction, conforming to the installation of the rock-breaker in our previous study [53]. Table 1 presents the parameters of the numerical example and the material of each part is steel, which has a density of 7800 kg/m3. The other properties of steel are shown in Table 2. The study cases include different amplitudes and frequencies of input force, as well as varying clearance sizes.Figure 11
Schematic overview of numerical example with 3D clearance joints.Table 1
Parameters of hydraulic cylinder.
Item
Value
Mass of slider and moving support (kg)
325.2
Mass of rod (kg)
26.1
Mass of cylinder barrel (kg)
37.2
Moment of inertia of rod (kg·m2)
2.63
Moment of inertia of cylinder barrel (kg·m2)
4.5
Piston area of forward camber (mm2)
5675
Piston area of return camber (mm2)
3299
Initial volume of forward camber (L)
2.554
Initial volume of return camber (L)
2.144Table 2
Geometric and physical parameters of revolute clearance joints.
Item
Left joint
Right joint
Radius of bearing (mm)
25
25
Young’s modulus (GPa)
207
207
Poisson’s ratio
0.3
0.3
Friction coefficientcf
0.01
0.01
Restitution coefficient
0.9
0.9There are two similar 3D joints with axial and radial clearances on the bilateral side of the hydraulic cylinder. Table2 presents the geometrical and physical parameters of the 3D joints. A schematic of the joints is provided in Figure 2.In the experiment and simulation, the forceFinput is a sine wave in each cycle, which can be expressed as(48)Finputt=Fmaxsinωt,where Fmax is the maximum input force and was set to 1, 2, 3, and 4 kN, while ω is the period parameter for determining the frequency of the force and was selected as 5, 10, 15, and 20 Hz.


### 3.2. Test Rig

A photograph of the experimental test rig is presented in Figure12. The hydraulic cylinder was selected from the middle arm cylinder of the XE60D excavator, manufactured by the XCMG group. The Hall sensors were made by Hefei Bangli Electronic Co., Ltd. The signals form the Hall sensors were collected with a PCI-9111HR data acquisition card from ADLINK in the data acquisition system, with a sampling rate up to 100 kHz. Thereafter, all signals were calibrated and displayed in the LabView software. The vibration bench, from SuShi Co., Ltd, could generate horizontal excitation. The input force Finput of the excitation was a simple harmonic wave and limited to 5000 N. A fixed support was mounted on the ground on the opposite side of the vibration bench. Two connecting links were installed between the fixed support and vibration bench, as illustrated in Figure 12, to provide greater stiffness to the fixed support during the experimental process.Figure 12
Photograph of test rig.Five Hall sensors were installed around the outer surface of the right bearing, as illustrated in Figure13. Two Hall sensors were positioned in the orthogonal direction to the journal as a group for detecting one side of the ring center movement of the bearing shown in Figure 2. Two other sensors, which are mutually perpendicular to the first two sensors and the journal, were used to measure the displacement on the other side of the ring center. The fifth Hall sensor was mounted between the two sides for detecting the axial clearance distance and was incorporated with a steel block fixed on the bearing. Magnets were mounted on the back of the Hall sensors and provided magnetic fields. The intension of the magnetic field changed with the distance between the sensor and the steel parts. An additional five sensors were positioned in the same manner on the left joint. The 11th sensor was mounted on the rod to detect the hydraulic cylinder movement, as illustrated in Figure 10, which respected the movement of the rod in the cylinder.Figure 13
Installation positions of five sensors in the right 3D joint.


## 3.1. Numerical Example Schema

A schematic overview of the numerical example is presented in Figure11. The hydraulic cylinder is inclined at 5° in the anticlockwise direction, conforming to the installation of the rock-breaker in our previous study [53]. Table 1 presents the parameters of the numerical example and the material of each part is steel, which has a density of 7800 kg/m3. The other properties of steel are shown in Table 2. The study cases include different amplitudes and frequencies of input force, as well as varying clearance sizes.Figure 11
Schematic overview of numerical example with 3D clearance joints.Table 1
Parameters of hydraulic cylinder.
Item
Value
Mass of slider and moving support (kg)
325.2
Mass of rod (kg)
26.1
Mass of cylinder barrel (kg)
37.2
Moment of inertia of rod (kg·m2)
2.63
Moment of inertia of cylinder barrel (kg·m2)
4.5
Piston area of forward camber (mm2)
5675
Piston area of return camber (mm2)
3299
Initial volume of forward camber (L)
2.554
Initial volume of return camber (L)
2.144Table 2
Geometric and physical parameters of revolute clearance joints.
Item
Left joint
Right joint
Radius of bearing (mm)
25
25
Young’s modulus (GPa)
207
207
Poisson’s ratio
0.3
0.3
Friction coefficientcf
0.01
0.01
Restitution coefficient
0.9
0.9There are two similar 3D joints with axial and radial clearances on the bilateral side of the hydraulic cylinder. Table2 presents the geometrical and physical parameters of the 3D joints. A schematic of the joints is provided in Figure 2.In the experiment and simulation, the forceFinput is a sine wave in each cycle, which can be expressed as(48)Finputt=Fmaxsinωt,where Fmax is the maximum input force and was set to 1, 2, 3, and 4 kN, while ω is the period parameter for determining the frequency of the force and was selected as 5, 10, 15, and 20 Hz.


## 3.2. Test Rig

A photograph of the experimental test rig is presented in Figure12. The hydraulic cylinder was selected from the middle arm cylinder of the XE60D excavator, manufactured by the XCMG group. The Hall sensors were made by Hefei Bangli Electronic Co., Ltd. The signals form the Hall sensors were collected with a PCI-9111HR data acquisition card from ADLINK in the data acquisition system, with a sampling rate up to 100 kHz. Thereafter, all signals were calibrated and displayed in the LabView software. The vibration bench, from SuShi Co., Ltd, could generate horizontal excitation. The input force Finput of the excitation was a simple harmonic wave and limited to 5000 N. A fixed support was mounted on the ground on the opposite side of the vibration bench. Two connecting links were installed between the fixed support and vibration bench, as illustrated in Figure 12, to provide greater stiffness to the fixed support during the experimental process.Figure 12
Photograph of test rig.Five Hall sensors were installed around the outer surface of the right bearing, as illustrated in Figure13. Two Hall sensors were positioned in the orthogonal direction to the journal as a group for detecting one side of the ring center movement of the bearing shown in Figure 2. Two other sensors, which are mutually perpendicular to the first two sensors and the journal, were used to measure the displacement on the other side of the ring center. The fifth Hall sensor was mounted between the two sides for detecting the axial clearance distance and was incorporated with a steel block fixed on the bearing. Magnets were mounted on the back of the Hall sensors and provided magnetic fields. The intension of the magnetic field changed with the distance between the sensor and the steel parts. An additional five sensors were positioned in the same manner on the left joint. The 11th sensor was mounted on the rod to detect the hydraulic cylinder movement, as illustrated in Figure 10, which respected the movement of the rod in the cylinder.Figure 13
Installation positions of five sensors in the right 3D joint.


## 4. Calculation and Analysis of Results

The submodel of the joint with 3D clearance was one part of the studied mechanism, providing the basis of the simulation and experiment. In the following sections, the effect of the ellipse eccentricity on the contact angle error and its innovation are first investigated. Thereafter, the experimental verification of the simulation and dynamics is analyzed, following which the effects of other parameters are studied.


### 4.1. Effect of Eccentricity

This section presents the effect of the eccentricity on the contact pointQ, which was calculated by the gradient descent method. The θerr of the potential point Q in the domain (0, π) when the axial and radial clearances were the same at 0.5 mm is illustrated, while θerr in the domain (0, π) is demonstrated for the symmetric quadrants. As illustrated in Figure 14, the error θerr was nearly a complete sinusoidal curve. Moreover, the error reached maximum value at approximately π/4 and 3π/4. Thereafter, the error decreased to zero with the change in the angle. The error disappeared when the angles were 0, π/2, and π. In the domain (π/2, π), the absolute value of θerr was a mirror of the above curve based on the symmetry, but the signal was negative. The maximum amplitudes of θerr exhibited accelerated growth when the eccentricity of the ellipse increased, as indicated in Figure 14(b). Therefore, we can conclude that the angle error θerr is minuscule when the clearance is very small and the circular model and the elliptical model are very similar to each other. Despite the fact that the angle error θerr is not obvious, this study selects the elliptical model in the simulation for future cases with larger clearances.Figure 14
(a) Angle errorθerr with rotation angle and (b) maximum angle error with eccentricity.
(a)
(b)


### 4.2. Experimental Verification

The test rig model was a cosimulation with ADAMS and MATLAB. ADAMS used the HHT integral solver and calculated the position, velocity, angle, and angle velocity for implementation in MATLAB. MATLAB then translated the ADAMS data to determine potential contact points. If there was a contact point, MATLAB would provide the impact forces and torques to the parts in ADAMS, which changed the rigid body acceleration. Meanwhile, the hydraulic cylinder force was also calculated in MATLAB. The integration tolerance in the simulation was less than1×10−7, and the integral step was set to 1×10−5 s to provide sufficient accuracy. Finally, the experimental and simulation times were each set to 10 s to reach a steady state. The initial conditions were as follows: The input force is 3 kN, and the frequency is 10 Hz. The journal is at the center of the bearing at the start of the simulation, the initial velocity is 0, and the entire test rig is in a 9.8 m/s2 gravity field. The axial clearance and radial clearance were both 0.5 mm. After the simulation was complete, we used MATLAB to obtain the simulation results and compiled a few scripts to draw the figures. A comparison of the experimental results with the simulations in the steady state is presented in Figure 15, where two periods are selected as an example, from 5.35 to 5.55 s. These results are similar to those of a previous study [34, 53]. Figures 15(a)–15(d) illustrate the two-side ring center experimental trajectories in red lines, which are detected by sensors 1–4 and 6–9. The trajectory of ideal joint is green points at the centers. The experimental data exhibited similar trends to those of the simulation considering clearance. However, the rod displacement of the simulation considering clearance was significantly different compared to that obtained with the model with ideal joints, which obtained different magnifications and shapes of the rod displacement curve. Therefore, as the experimental data agree well with the simulation results, the simulations are verified.Figure 15
Trajectories and rod displacements in the numerical and experimental studies withca = 0.5 and cr = 0.5: (a) outside trajectories in the right journal, (b) inside trajectories in the right journal, (c) outside trajectories in the left journal, (d) inside trajectories in the left journal, and (e) displacement zrod of the rod.
(a)
(b)
(c)
(d)
(e)Figure16 depicts the two-side simulation clearance trajectories. This result was similar to that presented in the literature [27]. The axles are illustrated as black lines that are nearly parallel to the horizon. The angle β of the Euler translation represents the angle between the journal and bearing axles, as illustrated in Figure 17. The angle β was almost equal to zero most of the time and only changed in the free state, which means that the two clearance trajectories were the same at most times. Hence, the center point trajectory of the journal, namely, the middle of the two-side point trajectories, can demonstrate the detailed dynamics of the hydraulic cylinder system.Figure 16
Trajectories and bearing axles of 3D joints in simulation.Figure 17
Euler angleβ for the two joints in simulation and experiment.We investigated one period as an example. The journal trajectories can be classified into four modes: freedom, collision, contact, and occasional rebounding, as illustrated in Figure18. These results are similar to those obtained in previous studies [27, 35, 53]. The slider moved toward the inward direction, and the journal mounted on the slider moved in the same direction. The hydraulic cylinder was also involved in the shrink process. The journal slid on the inner surface of the bearing indicated by the red line, as illustrated in Figures 18(b) and 18(d). The two parts of joint were maintained in contact each other without separation owing to the coupling of the hydraulic cylinder gravity and contact force. This mode is defined as the slide contact mode. The slide contact force, indicated by the red line and points, increased gradually in the slide mode, as shown in Figures 18(a) and 18(e). Meanwhile, the force on the piston in Figure 18(c) was similar to the slide contact force and showed certain friction peaks. The contact mode demonstrated in Figures 18(b), 18(d), 18(g), and 18(i), in which the journal reached the bearing wall boundary and remained in contact with a small area of the bearing. The contact force of the contact mode, indicated by green, was larger than that of other modes, as shown in Figures 18(a), 18(e), 18(f), and 18(g); this was the same as the force on the piston. Subsequently, the journal moved toward the outward direction of the bearing wall. When the journal moved away from the bearing wall, the trajectories were inside the clearance circles and remained in the free mode. The contact force, indicated by blue, was equal to zero as illustrated in Figures 18(a), 18(e), 18(f), and 18(g), and the force on the piston showed some fluctuation as demonstrated in Figure 18(h). The journal did not contact the bearing wall, and no contact force was imparted to the two joint parts. The rod simultaneously moved in a simple harmonic state owing to the spring-mass system of the hydraulic cylinder, which was derived from the bulk modulus βe of the oil [53]. When the free mode was completed, the journal reached the bearing wall, and an impact force peak emerged, indicating the rebound state. Thereafter, the contact mode emerged as the contact area was limited to a small arc section.Figure 18
Contact forces, trajectories, and forces on the piston in the numerical study withca = 0.5 and cr = 0.5: (a, f) contact forces in the left journal; (b, g) trajectories in the left journal; (c, h) forces on the piston; (d, i) trajectories in the right journal; (e, j) contact forces in the right journal.
(a)
(b)
(c)
(d)
(e)
(f)
(g)
(h)
(i)
(j)The contact mode is analyzed in Figure19. The rod displacement zrod and the centers of bearing distance were adjusted near to zero. The clearances in the joints were 0.5 mm each. Thus, the sum boundary of the two clearances was 1 mm, as indicated by the red dashed line. When the rod reached the bearing wall and was maintained in a contact mode, the rod displacement added or subtracted the sum boundary of 1 mm (red dashed line) and then corresponded to the bearing center distance.Figure 19
Displacementzrod of rod in simulation with ca = 0.5 and cr = 0.5.Figure20 illustrates the journal center as a point cloud distribution, which is similar to the results obtained by Yan et al. [27]. One in every twenty points was selected for representation to enable easier viewing of the points. The points exhibited greater contact opportunities on the bottom and bilateral sides of the cylinder, whose radius was the radial clearance and the length was the axial clearance. Therefore, it can be speculated that the bilateral sides of the bearing wall will experience more serious wear.Figure 20
Point distribution of middle joint in simulation withca = 0.5 and cr = 0.5.Figure21 illustrates the axial displacement of the journal center for the entire experimental and simulation period. When the bearing contacted the journal flanges, the contact force emerged on the corresponding side. It can thus be concluded that the displacement was irregular and contact occurred occasionally. Moreover, the axial contact force was less than the radial contact force. Thus, it is speculated that the wear on the journal flanges is not serious.Figure 21
Displacements and contact forces between journal flanges and bearing bases in the axial direction during simulation withca = 0.5 and cr = 0.5.


### 4.3. Effect of Input Forces

Figures22 and 23 illustrate that the increased input force acting on the slider changed the rod movement; the maximum value of the input force Fmax was selected as 1, 2, 3, and 4 kN with an input force frequency of 10 Hz. The initial values of the radial and axial clearances were both equal to 0.5 mm. The trajectories of the free state were longer and closer to the clearance circle center. The contact state increased gradually, as indicated in Figure 24. Meanwhile, the slide state descended. In contrast, the proportion of the free mode increased. The zrod peaks became increasingly larger with increasing input forces. Furthermore, the distance zrod exhibited two peaks when the slider moved inward and the space within the hydraulic cylinder shrank. The kinetic energy of the moving parts increased with the change in the input force, thus generating the ascending peaks. Increasing the input force will therefore deteriorate the hydraulic cylinder dynamic response.Figure 22
Trajectories in numerical simulation and experiment with differentFmax values: (a) left trajectories at Fmax = 1 kN, (b) right trajectories at Fmax = 1 kN, (c) left trajectories at Fmax = 2 kN, (d) right trajectories at Fmax = 2 kN, (e) left trajectories at Fmax = 3 kN, (f) right trajectories at Fmax = 3 kN, (g) left trajectories at Fmax = 4 kN, and (h) right trajectories at Fmax = 4 kN.
(a)
(b)
(c)
(d)
(e)
(f)
(g)
(h)Figure 23
Displacementszrod in numerical and experimental studies with different Fmax values: (a) Fmax = 1 kN, (b) Fmax = 2 kN, (c) Fmax = 3 kN, and (d) Fmax = 4 kN.
(a)
(b)
(c)
(d)Figure 24
Statistical results of time in numerical and experimental studies with differentFmax values.


### 4.4. Effect of Frequency

Figures25–27 illustrate the effects of the input force frequencies 5, 10, 15, and 20 Hz on the rod movement. The maximum value of the input force Fmax was 3 kN. The initial values of the radial and axial clearances were the same as those used in the previous simulations (cr = 0.5 mm and ca = 0.5 mm). The trajectories and displacement zrod were more regular, and the free mode was closer to the clearance circle center. Moreover, the proportions of the free and contact modes increased with the frequency change corresponding to a reduction in the slide time, as indicated in Figure 26. When the frequency increased, the peak magnitudes of zrod decreased gradually. The reason for the descending peaks is that the slider movement time was shorter with the kinetic energy of the parts being reduced. One can thus conclude that the input force frequency reduces the deterioration of the dynamic response.Figure 25
Trajectories in numerical and experimental studies with different frequencies: (a) left trajectories at 5 Hz, (b) right trajectories at 5 Hz, (c) left trajectories at 10 Hz, (d) right trajectories at 10 Hz, (e) left trajectories at 15 Hz, (f) right trajectories at 15 Hz, (g) left trajectories at 20 Hz, and (h) right trajectories at 20 Hz.
(a)
(b)
(c)
(d)
(e)
(f)
(g)
(h)Figure 26
Displacementszrod in numerical and experimental studies with different frequencies: (a) 5 Hz, (b) 10 Hz, (c) 15 Hz, and (d) 20 Hz.
(a)
(b)
(c)
(d)Figure 27
Statistical results of time in numerical and experimental studies with different frequencies.


### 4.5. Effect of Clearance Size

It was necessary to analyze the influence of the clearance size because the clearance increases gradually during the servicing period of a hydraulic cylinder. The maximum value of the input forceFmax was 3 kN, and the frequency of the input force was 10 Hz. The clearance size cr was adjusted by varying the journal radius of the joints, selected as 0.2, 0.5, and 1 mm. The trajectories at cr = 0.2 mm were the most random for the small clearance owing to the rapid bearing collision, as shown in Figure 28. At cr = 0.2 mm, the displacement zrod in the contact mode was minimal and the least vibrations occurred in every period. When the clearance increased to 1 mm, the displacement zrod was maximal and the most vibrations emerged. The displacement and vibration when the clearance was 0.5 mm were within the two clearance values above, as illustrated in Figure 29. The contact state time decreased in the study, while the free mode time and the slide state increased, as shown in Figure 30. The reason for the above phenomenon is that expanding clearance provides more space for the kinetic energy of the slider and rod to be accumulated. This is also the cause for additional vibrations in the free mode and the shorter slide state time. The hydraulic cylinder dynamic response will therefore degrade with increase in the clearance size.Figure 28
Trajectories in numerical and experimental studies with differentcr: (a) left trajectories at cr = 0.2, (b) right trajectories at cr = 0.2, (c) left trajectories at cr = 0.5, (d) right trajectories at cr = 0.5, (e) left trajectories at cr = 1.0, and (f) right trajectories at cr = 1.0.
(a)
(b)
(c)
(d)
(e)
(f)Figure 29
Displacementszrod in numerical and experimental studies with different cr values: (a) cr = 0.2, (b) cr = 0.5, and (c) cr = 1.0.
(a)
(b)
(c)Figure 30
Statistical results of time in numerical and experimental studies with differentcr values.The results of this study demonstrate that the dynamics of the 3D hydraulic cylinder are affected by variations in the parameters. A comparison of the effects of the input force size, frequency, and clearance size in Figures24, 27, and 30 shows that the clearance size could drastically reduce the contact mode time in simulation and experiment results. However, the contact mode time increases with an increase in the input force size and frequency. The free mode exhibits a peak with a frequency change, which increases with the magnitude of input force and clearance size. The slide mode increases gradually with the input force and frequency. However, the time of sliding decreases when the magnitude of the input force and the frequency increase. These findings are of significance for the design and application of hydraulic cylinders.


## 4.1. Effect of Eccentricity

This section presents the effect of the eccentricity on the contact pointQ, which was calculated by the gradient descent method. The θerr of the potential point Q in the domain (0, π) when the axial and radial clearances were the same at 0.5 mm is illustrated, while θerr in the domain (0, π) is demonstrated for the symmetric quadrants. As illustrated in Figure 14, the error θerr was nearly a complete sinusoidal curve. Moreover, the error reached maximum value at approximately π/4 and 3π/4. Thereafter, the error decreased to zero with the change in the angle. The error disappeared when the angles were 0, π/2, and π. In the domain (π/2, π), the absolute value of θerr was a mirror of the above curve based on the symmetry, but the signal was negative. The maximum amplitudes of θerr exhibited accelerated growth when the eccentricity of the ellipse increased, as indicated in Figure 14(b). Therefore, we can conclude that the angle error θerr is minuscule when the clearance is very small and the circular model and the elliptical model are very similar to each other. Despite the fact that the angle error θerr is not obvious, this study selects the elliptical model in the simulation for future cases with larger clearances.Figure 14
(a) Angle errorθerr with rotation angle and (b) maximum angle error with eccentricity.
(a)
(b)


## 4.2. Experimental Verification

The test rig model was a cosimulation with ADAMS and MATLAB. ADAMS used the HHT integral solver and calculated the position, velocity, angle, and angle velocity for implementation in MATLAB. MATLAB then translated the ADAMS data to determine potential contact points. If there was a contact point, MATLAB would provide the impact forces and torques to the parts in ADAMS, which changed the rigid body acceleration. Meanwhile, the hydraulic cylinder force was also calculated in MATLAB. The integration tolerance in the simulation was less than1×10−7, and the integral step was set to 1×10−5 s to provide sufficient accuracy. Finally, the experimental and simulation times were each set to 10 s to reach a steady state. The initial conditions were as follows: The input force is 3 kN, and the frequency is 10 Hz. The journal is at the center of the bearing at the start of the simulation, the initial velocity is 0, and the entire test rig is in a 9.8 m/s2 gravity field. The axial clearance and radial clearance were both 0.5 mm. After the simulation was complete, we used MATLAB to obtain the simulation results and compiled a few scripts to draw the figures. A comparison of the experimental results with the simulations in the steady state is presented in Figure 15, where two periods are selected as an example, from 5.35 to 5.55 s. These results are similar to those of a previous study [34, 53]. Figures 15(a)–15(d) illustrate the two-side ring center experimental trajectories in red lines, which are detected by sensors 1–4 and 6–9. The trajectory of ideal joint is green points at the centers. The experimental data exhibited similar trends to those of the simulation considering clearance. However, the rod displacement of the simulation considering clearance was significantly different compared to that obtained with the model with ideal joints, which obtained different magnifications and shapes of the rod displacement curve. Therefore, as the experimental data agree well with the simulation results, the simulations are verified.Figure 15
Trajectories and rod displacements in the numerical and experimental studies withca = 0.5 and cr = 0.5: (a) outside trajectories in the right journal, (b) inside trajectories in the right journal, (c) outside trajectories in the left journal, (d) inside trajectories in the left journal, and (e) displacement zrod of the rod.
(a)
(b)
(c)
(d)
(e)Figure16 depicts the two-side simulation clearance trajectories. This result was similar to that presented in the literature [27]. The axles are illustrated as black lines that are nearly parallel to the horizon. The angle β of the Euler translation represents the angle between the journal and bearing axles, as illustrated in Figure 17. The angle β was almost equal to zero most of the time and only changed in the free state, which means that the two clearance trajectories were the same at most times. Hence, the center point trajectory of the journal, namely, the middle of the two-side point trajectories, can demonstrate the detailed dynamics of the hydraulic cylinder system.Figure 16
Trajectories and bearing axles of 3D joints in simulation.Figure 17
Euler angleβ for the two joints in simulation and experiment.We investigated one period as an example. The journal trajectories can be classified into four modes: freedom, collision, contact, and occasional rebounding, as illustrated in Figure18. These results are similar to those obtained in previous studies [27, 35, 53]. The slider moved toward the inward direction, and the journal mounted on the slider moved in the same direction. The hydraulic cylinder was also involved in the shrink process. The journal slid on the inner surface of the bearing indicated by the red line, as illustrated in Figures 18(b) and 18(d). The two parts of joint were maintained in contact each other without separation owing to the coupling of the hydraulic cylinder gravity and contact force. This mode is defined as the slide contact mode. The slide contact force, indicated by the red line and points, increased gradually in the slide mode, as shown in Figures 18(a) and 18(e). Meanwhile, the force on the piston in Figure 18(c) was similar to the slide contact force and showed certain friction peaks. The contact mode demonstrated in Figures 18(b), 18(d), 18(g), and 18(i), in which the journal reached the bearing wall boundary and remained in contact with a small area of the bearing. The contact force of the contact mode, indicated by green, was larger than that of other modes, as shown in Figures 18(a), 18(e), 18(f), and 18(g); this was the same as the force on the piston. Subsequently, the journal moved toward the outward direction of the bearing wall. When the journal moved away from the bearing wall, the trajectories were inside the clearance circles and remained in the free mode. The contact force, indicated by blue, was equal to zero as illustrated in Figures 18(a), 18(e), 18(f), and 18(g), and the force on the piston showed some fluctuation as demonstrated in Figure 18(h). The journal did not contact the bearing wall, and no contact force was imparted to the two joint parts. The rod simultaneously moved in a simple harmonic state owing to the spring-mass system of the hydraulic cylinder, which was derived from the bulk modulus βe of the oil [53]. When the free mode was completed, the journal reached the bearing wall, and an impact force peak emerged, indicating the rebound state. Thereafter, the contact mode emerged as the contact area was limited to a small arc section.Figure 18
Contact forces, trajectories, and forces on the piston in the numerical study withca = 0.5 and cr = 0.5: (a, f) contact forces in the left journal; (b, g) trajectories in the left journal; (c, h) forces on the piston; (d, i) trajectories in the right journal; (e, j) contact forces in the right journal.
(a)
(b)
(c)
(d)
(e)
(f)
(g)
(h)
(i)
(j)The contact mode is analyzed in Figure19. The rod displacement zrod and the centers of bearing distance were adjusted near to zero. The clearances in the joints were 0.5 mm each. Thus, the sum boundary of the two clearances was 1 mm, as indicated by the red dashed line. When the rod reached the bearing wall and was maintained in a contact mode, the rod displacement added or subtracted the sum boundary of 1 mm (red dashed line) and then corresponded to the bearing center distance.Figure 19
Displacementzrod of rod in simulation with ca = 0.5 and cr = 0.5.Figure20 illustrates the journal center as a point cloud distribution, which is similar to the results obtained by Yan et al. [27]. One in every twenty points was selected for representation to enable easier viewing of the points. The points exhibited greater contact opportunities on the bottom and bilateral sides of the cylinder, whose radius was the radial clearance and the length was the axial clearance. Therefore, it can be speculated that the bilateral sides of the bearing wall will experience more serious wear.Figure 20
Point distribution of middle joint in simulation withca = 0.5 and cr = 0.5.Figure21 illustrates the axial displacement of the journal center for the entire experimental and simulation period. When the bearing contacted the journal flanges, the contact force emerged on the corresponding side. It can thus be concluded that the displacement was irregular and contact occurred occasionally. Moreover, the axial contact force was less than the radial contact force. Thus, it is speculated that the wear on the journal flanges is not serious.Figure 21
Displacements and contact forces between journal flanges and bearing bases in the axial direction during simulation withca = 0.5 and cr = 0.5.


## 4.3. Effect of Input Forces

Figures22 and 23 illustrate that the increased input force acting on the slider changed the rod movement; the maximum value of the input force Fmax was selected as 1, 2, 3, and 4 kN with an input force frequency of 10 Hz. The initial values of the radial and axial clearances were both equal to 0.5 mm. The trajectories of the free state were longer and closer to the clearance circle center. The contact state increased gradually, as indicated in Figure 24. Meanwhile, the slide state descended. In contrast, the proportion of the free mode increased. The zrod peaks became increasingly larger with increasing input forces. Furthermore, the distance zrod exhibited two peaks when the slider moved inward and the space within the hydraulic cylinder shrank. The kinetic energy of the moving parts increased with the change in the input force, thus generating the ascending peaks. Increasing the input force will therefore deteriorate the hydraulic cylinder dynamic response.Figure 22
Trajectories in numerical simulation and experiment with differentFmax values: (a) left trajectories at Fmax = 1 kN, (b) right trajectories at Fmax = 1 kN, (c) left trajectories at Fmax = 2 kN, (d) right trajectories at Fmax = 2 kN, (e) left trajectories at Fmax = 3 kN, (f) right trajectories at Fmax = 3 kN, (g) left trajectories at Fmax = 4 kN, and (h) right trajectories at Fmax = 4 kN.
(a)
(b)
(c)
(d)
(e)
(f)
(g)
(h)Figure 23
Displacementszrod in numerical and experimental studies with different Fmax values: (a) Fmax = 1 kN, (b) Fmax = 2 kN, (c) Fmax = 3 kN, and (d) Fmax = 4 kN.
(a)
(b)
(c)
(d)Figure 24
Statistical results of time in numerical and experimental studies with differentFmax values.


## 4.4. Effect of Frequency

Figures25–27 illustrate the effects of the input force frequencies 5, 10, 15, and 20 Hz on the rod movement. The maximum value of the input force Fmax was 3 kN. The initial values of the radial and axial clearances were the same as those used in the previous simulations (cr = 0.5 mm and ca = 0.5 mm). The trajectories and displacement zrod were more regular, and the free mode was closer to the clearance circle center. Moreover, the proportions of the free and contact modes increased with the frequency change corresponding to a reduction in the slide time, as indicated in Figure 26. When the frequency increased, the peak magnitudes of zrod decreased gradually. The reason for the descending peaks is that the slider movement time was shorter with the kinetic energy of the parts being reduced. One can thus conclude that the input force frequency reduces the deterioration of the dynamic response.Figure 25
Trajectories in numerical and experimental studies with different frequencies: (a) left trajectories at 5 Hz, (b) right trajectories at 5 Hz, (c) left trajectories at 10 Hz, (d) right trajectories at 10 Hz, (e) left trajectories at 15 Hz, (f) right trajectories at 15 Hz, (g) left trajectories at 20 Hz, and (h) right trajectories at 20 Hz.
(a)
(b)
(c)
(d)
(e)
(f)
(g)
(h)Figure 26
Displacementszrod in numerical and experimental studies with different frequencies: (a) 5 Hz, (b) 10 Hz, (c) 15 Hz, and (d) 20 Hz.
(a)
(b)
(c)
(d)Figure 27
Statistical results of time in numerical and experimental studies with different frequencies.


## 4.5. Effect of Clearance Size

It was necessary to analyze the influence of the clearance size because the clearance increases gradually during the servicing period of a hydraulic cylinder. The maximum value of the input forceFmax was 3 kN, and the frequency of the input force was 10 Hz. The clearance size cr was adjusted by varying the journal radius of the joints, selected as 0.2, 0.5, and 1 mm. The trajectories at cr = 0.2 mm were the most random for the small clearance owing to the rapid bearing collision, as shown in Figure 28. At cr = 0.2 mm, the displacement zrod in the contact mode was minimal and the least vibrations occurred in every period. When the clearance increased to 1 mm, the displacement zrod was maximal and the most vibrations emerged. The displacement and vibration when the clearance was 0.5 mm were within the two clearance values above, as illustrated in Figure 29. The contact state time decreased in the study, while the free mode time and the slide state increased, as shown in Figure 30. The reason for the above phenomenon is that expanding clearance provides more space for the kinetic energy of the slider and rod to be accumulated. This is also the cause for additional vibrations in the free mode and the shorter slide state time. The hydraulic cylinder dynamic response will therefore degrade with increase in the clearance size.Figure 28
Trajectories in numerical and experimental studies with differentcr: (a) left trajectories at cr = 0.2, (b) right trajectories at cr = 0.2, (c) left trajectories at cr = 0.5, (d) right trajectories at cr = 0.5, (e) left trajectories at cr = 1.0, and (f) right trajectories at cr = 1.0.
(a)
(b)
(c)
(d)
(e)
(f)Figure 29
Displacementszrod in numerical and experimental studies with different cr values: (a) cr = 0.2, (b) cr = 0.5, and (c) cr = 1.0.
(a)
(b)
(c)Figure 30
Statistical results of time in numerical and experimental studies with differentcr values.The results of this study demonstrate that the dynamics of the 3D hydraulic cylinder are affected by variations in the parameters. A comparison of the effects of the input force size, frequency, and clearance size in Figures24, 27, and 30 shows that the clearance size could drastically reduce the contact mode time in simulation and experiment results. However, the contact mode time increases with an increase in the input force size and frequency. The free mode exhibits a peak with a frequency change, which increases with the magnitude of input force and clearance size. The slide mode increases gradually with the input force and frequency. However, the time of sliding decreases when the magnitude of the input force and the frequency increase. These findings are of significance for the design and application of hydraulic cylinders.


## 5. Conclusions

The present study was concerned with investigation of the dynamic responses of the hydraulic cylinder mechanism with a 3D joint model considering the radial and axial clearances. A simulation model for the hydraulic cylinder with 3D joints was presented, including the potential contact scenarios and points, line, and plate contact force models, and the addition of the effective bulk modulus of hydraulic oil, as depicted in our previous work. We also simplified the line contact force calculation. Subsequently, an experiment verified the simulation results. The dynamics of the hydraulic cylinder model considering 3D clearance in the simulation are closer to the experimental results compared to those obtained by the model with ideal joints. In contrast, the disadvantage is the increased amount of calculations. The main conclusions are summarized as follows:(1)
The effect of eccentricity on the contact point was thoroughly analyzed. The results demonstrate that the angle error in the domain (0,π) was nearly a complete sinusoidal curve. The error disappeared only when the angle = 0, π/2, π. The maximum amplitudes of the angle error ascended when the eccentricity increased. It can therefore be concluded that the angle error of the contact point is not ignored when the misalignment is large compared to a traditional 3D joint with clearances; in addition, more accurate collision points are provided.(2)
The dynamic response of the rod can be classified into four states: free, rebound, slide, and contact. The stiffness of the hydraulic cylinder caused the rod to vibrate rapidly in the free model. The journals slid on the bottom of the bearing wall because of the gravity and low velocity in the slide mode. When the journal reached the boundary together with the bearing at the lines, the joint was in the contact state.(3)
An increase in the input force and clearance size will degrade the dynamic response of the hydraulic cylinder, for example, by increasing the time of the free state and the peak of the rod displacement. The frequency of the input force causes the trajectories and displacement to be more regular and can reduce the deterioration of the dynamic response of the hydraulic cylinder.This work provides greater insight into the changing dynamics of the hydraulic cylinder with 3D clearance joints and offers theoretical support for further studies on the hydraulic cylinder.



---

*Source: 1204328-2019-12-17.xml*