A numerical coupling model for a multibody system with multiple lubricated clearance joints

In machines and mechanisms, most of the clearance joints are designed to operate with lubricant fluid, which is an effective way of ensuring better performance of the mechanical systems. The differently located lubricated joints interact with each other, influencing the total performances of the mechanical system. In this study, a numerical coupling model for the dynamics and lubrication analyses of the multibody system with multiple lubricated clearance joints is established. This approach couples the dynamics model of the multibody system with the lubrication models of translational and revolute clearance joints. The dynamics model is established by the Lagrange’s method, while The lubrication models are solved using the finite element method. The hydrodynamic forces built up by the lubricant fluid are evaluated from the motion states of the system and then input into dynamics model of the multibody system. Finally, the proposed approach is applied to a four-stroke gasoline engine with the lubricated piston skirt-liner subsystem and lubricated bearing at the big end of the connecting rod.


Introduction
Most of the joints in machines and mechanisms are designed to operate with lubricant fluid, which is an effective way of ensuring better performance of the mechanical systems [1][2][3]. The differently located joints with lubricated clearances depend on and work on each other, influencing the total performances of the system. Therefore, an improved model of the multibody system with the consideration of the interactions among the multiple joints is required to achieve a better understanding of the dynamics and tribology performances of the mechanism.
Internal combustion engines (ICEs) are important applications, and the lubricated joints in them play crucial roles in maintaining their stability of motion and reducing the friction power loss. The translational joints at the piston assembly and the revolute joints at the bearings interact each other and cooperate to affect the dynamics and tribological performances of engines.
Over the last several decades, considerable studies have been carried out on the lubrication of the piston skirt-liner system and the performances of the dynamically loaded bearing, respectively. To illustrate for the investigation on the piston skirt-liner system, Li et al. [4] modeled the dynamics of the piston by incorporating the hydrodynamic lubrication model into the dynamics equations of the piston. Zhu et al. [5,6], Keribar and Dursunkaya [7], Wong et al. [8] and Liu et al. [9] established more perfect dynamics models of the piston based on Patir and Cheng's average Reynolds equation [10,11].
By taking dynamically loaded bearings as the lubricated revolute clearance joints, many endeavors have also investigated the performances of bearings. Flores et al. [12] presented hydrodynamic lubrication models for lubricated revolute joints in constrained rigid planar multibody systems by considering the infinitely long and short journal bearing conditions. Daniel et al. [13] proposed a model to analyze the dynamics performance of a slider-crank mechanism with lubricated clearance joint. Machado et al. [14] studied the effect of the lubricated revolute joint parameters and hydrodynamic force models on the dynamic response of planar multibody systems. Zhao et al. [15,16] established a dynamics model of the piston-connecting rod-crank system in a framework of the multibody mechanical system with lubricated joint.
The above investigations just focused on the mechanisms with single clearance joint. The investigation on the mechanisms with lubricated multiple joints, including the translational and revolution clearance joints, have been rarely reported.
The current study proposed a new numerical coupling method for the multibody mechanical system with lubricated translational and revolute clearance joints. This approach is established by coupling the dynamics model of the multibody system with the mixed lubrication models of multiple clearance joints. The Reynolds equation is solved based on the finite element method (FEM) due to its advantage of the flexibility in dealing with irregular domains [17]. In the end, the proposed approach is applied to a four-stroke gasoline engine with lubricated piston skirt-liner system and lubricated revolute clearance joint at the big end of the connecting rod.

The lubrication models of lubricated clearance joints in ices
In this section, the piston skirt-liner systems in ICEs can be considered as lubricated translational joints with clearance, while the bearings at the big end of the connecting rod can be considered as lubricated revolute clearance joints.

Hydrodynamic lubrication model of the piston skirt-liner system
In order to take the influences of the roughnesses of mating surfaces into account, the oil film pressure between the skirt and liner will be evaluated based on the average Reynolds equation [10,11] where P is the dynamic viscosity of the oil, and p u is the relative velocity between the piston and liner. The variables h and p are the oil film thickness and pressure. The local coordinate axis x and y are respectively along the circumferential direction and axial direction of the piston skirt as shown in Fig. 1 where 0 e and J denote the lateral displacement and tilt angle of the piston, respectively.The function skirt ( , ) h y D is used to define the piston skirt profile ( Fig.   1 (d)), and D is the angular coordinate of the piston skirt from the thrust side ( Fig. 1 (c)). According to Eq. (2), the squeeze part of the Reynolds equation can be defined as The dimensionless process for Eq. (1) can be written as where R is the radius of the piston and 1 L is the length of the crank. In this study, the angular velocity of the crank is chosen as the reference angular velocity c Z .
The dimensionless Reynolds equation can be written as Considering that the distribution of oil film pressure is symmetrical about the xoy plane, so only the zone from 0 D toD S is solved. The boundary conditions of the piston skirt lubrication model can be expressed as [21,22]  The eccentricity vector connecting the centers of journal and bearing can be calculated as where P i r and P j r are respectively the position vector of i P and j P in the global coordinate frame XOY .
The magnitude of the eccentricity vector can be obtained by The unit eccentricity vector along the eccentricity direction can be expressed as e e n = . ( The eccentricity ratio can be defined as The time rate of eccentricity ratio variation H H can be expressed as b c H c H en (11) where e is the velocity of j P relative to i P and can be obtained by differentiating Eq. (7) with respect to time. The average Reynolds equation is also adopted to evaluate the hydrodynamics oil forces at the joint and can be expressed as [15] 3 3 12 12 where, the parameter b where t is the unit tangential vector defined by rotating the vector n by an angle of 90° in the counterclockwise direction as shown in Fig. 2 (b). The dimensionless process for Eq. (12) can be written as where b is the width of the bearing.
The dimensionless Reynolds equation for Eq. (12) can be written as  I  I  T  T   I  I  Z I H T V  TH H TJ  T  It is assumed that the pressure gradient is equal to zero in the middle section due to the symmetry, and the pressure at both edges will decrease to zero. Thus, the boundary condition in the axial direction can be expressed as [15,23] , 0 2 In the circumferential direction, the Reynolds boundary condition is adopted due to its accuracy and can be written as [15,23] where * T is a floating boundary [23].

Finite Element Method for solving the lubrication models
Equations (5) and (17) are nonhomogeneous partial differential equations of elliptical type. Such equations will be solved currently by way of using FEM due to its flexibility in solving the elliptic equation with irregular domains [17]. Eqs. (5) and (17) can be rewritten with the same form as where for the piston skirt-liner system, 3 while for the lubricated revolute joint, According to the Galerkin theory, the integration over their respective whole lubrication domain : shown in Fig. 3 can be built as where 1 is the shape function of discrete elements.
After the integration using the Green's theorem, Eq. (23) can be expressed as where * is the boundary of each integral domain. In order to solve the above equation, the domain is discretized into a number of triangular elements (Fig. 3).
The shape function of the triangular element can be expressed as For the boundary elements with the second kind boundary conditions stated in Eqs. (6), (18) and (19), it can be easily proved that the integrations on the boundaries are equal to zero. While the first kind boundary conditions can be considered by setting the oil film pressure on the corresponding node to zero. Based on the above statement and by assembling the corresponding items of the element in Eq. (27), Eqs. (5) and (17) can be transformed into: The hydrodynamic oil film pressure distribution on the both skirt and bearing can be obtained by solving this equation with the consideration of the boundary conditions.

Asperity contact Forces
When the oil film thickness ratio / h V between two mating surfaces is below 4.0, the effect of solid-to-solid contact in mixed lubrication could be considered [24,25]. The asperity contact pressure ( ) c p x,y is given as where ' E is the composite elastic modulus. The dimensionless parameter K is assumed to be a constant with the value of 4 1.198 10 u [25,26]. The switch function 5 2 F H can be written as [24] The friction forces between two mating surfaces result from the lubricant shear and asperity contacts in mixed lubrication. The contact friction can be simply calculated with a boundary friction coefficient. While the average shear stress can be expressed as [10,11] 2 where f ) , fs ) and fp ) are shear stress factors [10,11].
3 Dynamic analysis of the piston-connecting rod-crank system Figure 4 depicts a piston-connecting rod-crank multibody system with the lubricated revolute and translational joints at the piston and big end bearing of the connecting rod, respectively. All the bodies are assumed to be rigid, and other joints are ideal without considering the existence of the friction and clearance. Subsequently, the dynamics equations that govern the of system is first formulated.
The generalized Cartesian coordinates of the system can be written as  The ideal joints in the system provide kinematic constraints to the connected components, while the clearance joints provide the dynamics constraints instead of the kinematic constraints [27]. The kinematic constraints Φ in this system due to the mechanical joints and specified motion trajectories can be written in a compact form as where r L determines the location of mass center on crank.
where q Φ is the Jacobian matrix of the constraint equations. q and q are vectors of the system velocity and acceleration, respectively.
The dynamics equations for the rigid bodies in a constrained multibody system can be written as [27,28] where the vector O denotes the Lagrange multipliers.
The term can be used to express the joint reaction forces [29]. The system mass matrix M is a diagonal matrix with diagonal elements as  T  1  1   2  2  3  3   ,  , , , , , , x y  1). In multibody mechanical systems, a lubricated clearance joint does not produce any kinematic constraint. Instead, it acts as a force element producing time-dependent forces. Thus, the lubricated joint produces, the so-called, force constraints. The hydrodynamic forces built up by the lubricant fluid are evaluated from the knowledge of the system variables as stated in Section 2 and then introduced as external generalized forces into the dynamics equations of the mechanical system.
Combining Eqs. (34) and (35), the dynamics equations of a system subjected to constraints can be stated in the form [1,[27][28][29][30][31]: ) ) ) . (38) According to Eq. (37), q and O can be solved with the proper initial conditions. However, the equation does not use the explicit position and velocity equations associated with the kinematic constraints, which implies that during simulation, the original constraint equations will be violated [1,14,30,32]. To control the constraints violation during numerical integration, the Baumgarte stabilization technique is adopted to damp out the acceleration constraint violations by feeding back the violations of positions and velocity constraints [1,14,30,[32][33][34]. Eq. (37) can be modified as whereand 9 are feedback parameters and should be adequately chosen [34,35]. The proposed model is performed on a four-stroke gasoline engine in this section. Some necessary input data for the simulation are listed in Table 1. Fig. 5 shows the combustion gas pressures acting on the piston crown respectively at the speeds of 1000, 2000, 4000 and 6000 rpm, which were obtained from experiments. The simulation parameters for this system are shown in Table   2. Here,the parameter : of the bearing. In the following calculations, the solving efficiency for hydrodynamic lubrication with the above mesh density is improved by about 65% compared with the case with the doubled element numbers, but the results remain about the same. Therefore, it can be proved that such mesh density is convergent and can guarantee a good accuracy. Convergence criterion 1%

Results and discussion
The case with the constant crank velocity of 6000 rpm is first simulated. Several numerical results are presented and analyzed to demonstrate the computational implementation of the proposed methodology, as well as to reveal the dynamics and lubrication performances of the system. In order to reveal the interactions between the lubricated piston skirt-liner system and lubricated bearing at connecting rod big end (LBRBE), several other previously proposed models are also introduced: the first model is named as "Only-Skirt model" in which only the piston skirt-liner joint is taken into account and other joints are considered to be ideals without clearances and friction as proposed in study [21]. The second model just considers the LBRBE in the multibody system as stated in [15,16] and is called as "Only-LBRBE model". The current coupling model considers both the lubricated piston skirt-liner system and LBRBE and is named as "Skirt-LBRBE model" for short. Figure 6. The kinematics results of the piston and connecting rod. Figure 6 shows the kinematics results of the piston and connecting rod. All the results presented are relative to two full crank rotations after steady-state has been reached. It can be found that the translation accelerations of the piston (i.e., 0 e 0 e and y y ) obtained from the Skirt-LBRBE model have more obvious oscillations than those from the Only-skirt model, especially in the horizontal direction. The reason may be that the LBRBE has influences on the oil film force acting at the piston skirt, or the reaction forces between the piston and connecting rod.
In order to better explain the influence of the LBRBE, Figs. 7 (a)-(c) shows the hydrodynamic oil film forces acting on the piston skirt. It is obvious that the plots of oil film forces from both Skirt-LBRBE and Only-Skirt models, including total normal thrust force, friction force as well as their moments relative to the center of the piston pin, overlap very well. This means that the LBRBE almost has no influence on the lubrication performance of the piston skirt-liner system. The results are calculated based on three different models: the Only-Skirt model, the Only-LBRBE model, and the Skirt-LBRBE model. By way of comparison, it can be found that compared with the Only-Skirt model, the other models with the LBRBE will introduce the extra degrees of freedom for the connecting rod, which will result in the oscillation of the reaction forces at the revolute joint, and thus affect the accelerations of the connecting rod (Figs. 6 (d)-(f)), as well as the reaction forces between the rod and piston.
By comparing the Skirt-LBRBE model with the Only-LBRBE model, it can be found that the absence of the clearance at piston skirt-liner will decrease the accelerations oscillation of connecting rod and thus the reaction forces at bearing (Figs. 6 and 7). In addition, from the comparison between the Skirt-LBRBE model with the Only-Skirt model, it can be found that due to the fact that rod and piston are connected by the piston pin, although the LBRBE will result in the oscillation of the connecting rod acceleration and reaction forces at the joint of the big end of the connecting rod, no extra moment about the piston pin will be introduced into the piston. This can also be observed in Figs. 6 (c) and 7 (c).  It is clear that the existence of the LBRBE almost has no effect on the lubrication of the piston skirt as observed in Fig. 7. However, for the lubrication performances of the LBRBE, there are obvious differences between the Skirt-LBRBE model and the Only-LBRBE model. Fig. 8 (d) shows that for most period of the time, MOFTs obtained from the both models agree well with each other. While at some special moments when the MOFT has very small values (i.e., at 62, 252 and 562 °C rank Angel (CA)), the MOFTs derived from the Only-LBRBE model have the smaller values.
The above phenomenon can also be inferred from Fig.  9 which presents the center trajectory of journal relative to the center of the bushing. It can be found that except for the moments of 62, 252 and 562 °CA where there are some obvious differences, the trajectory from the both models almost overlap. In addition, the journal will quickly pass by the center of bearing at about 361°CA, which causes the MOFT a transient sharply increase as shown in Fig. 8 (d).   Fig. 10 displays asperity contact pressures, as well as the oil film and asperity friction stresses at these two moments. It can be found that due to the smaller oil film thickness in the Only-LBRBE model, the higher oil film shear stress and more surface asperity contacts will be obtained, which will generate the larger friction force and FPL (Fig. 8 (e)).
From the comparison between Figs. 10 (c) and (f), it can be found that although the MOFT at 62 °CA is larger than that at 562 °CA, the friction from the oil film shear stress is smaller at the latter moment. The reason may be that the relative slip velocity between the surfaces of the journal and bearing at 562 °CA is smaller than that at 62 °CA , which limits the contribution of the lubricant shear stress to the friction as stated in Eq. (31). Furthermore, the lower relative slip velocity at 562 °CA will also result in the lower PLF ( Fig. 8 (e)). system is very limited. However, the lubricated skirt-liner system is obviously conducive to improving the lubrication performances of the LBRBE (i.e., the MOFT is increased and the FPL is reduced). This conclusion is consistent with that shown in Fig. 7. Subsequently, the dynamics and tribological performances for different crank speeds are investigated based on the Skirt-LBRBE coupling model. Fig. 11 presents the comparisons of the secondary motions of the piston for different speeds. It can be found that the amplitudes of piston secondary velocity (i.e., 0 e 0 e and J J ) tend to increase with the speed, and both reach the maximum values at about 387 °CA. The reason may be that in the case with higher speed, a smaller period of motion is required to complete the secondary motion of the piston, and thus a larger secondary velocity of the piston will be produced. Moreover, slap noise has a direct bearing on the secondary movement velocity of the piston. This means that the increase of the speeds will cause higher slap noise. With the increase of speed, the ranges of secondary motion of the piston (i.e., 0 e and J ) will also increase for the most of the time. While during the period from about 270 to 450 °CA, the amplitudes of secondary motion will decrease with the increase of the speed. In fact, 0 e and J depend on the combined effects of both the piston secondary velocity and the thrust forces acting on the piston skirt. Moreover, the higher secondary velocity and smaller thrust force will result in the decrease of secondary motion. Figure 12 (a) presents the thrust force normally acting on the piston skirt. By way of comparison, the same phenomenon as that of 0 e and J can also be observed. This is due to the fact that the thrust forces are mainly used to balance the forces which act on the piston in the horizontal direction and are mainly induced by the combustion force and the inertias of the piston and connecting rod (Fig. 13 (a)). Compared with the forces induced by the combustion force, the forces induced by the inertias of the piston and the connecting rod are dominant and increase with the crank speed for most of the time.  However, during the period from about 270 to 450 °CA, the combustion forces have high values and the thrust forces are mainly used to balance the forces induced by them. Simultaneously, the directions of the forces induced by the inertias of the piston and connecting rod are horizontally opposite with the directions of the forces induced by combustion forces, which will weaken the influence of the combustion force. With the increase of the crank speeds, this weakening effect will increasingly obvious and result in the decrease of the thrust forces, especially as the speed increase from 2000 to 6000 rpm. Figure 12 (b) lays out the comparison of the MOFT for different speeds. It can be found that the higher speed is accompanied with the smaller MOFT for most of the time. As stated in Eq. (31), the decrease of the oil film thickness and the increase of the relative velocity can increase the lubricant shear stress and thus leads to a rise in the friction force and friction power loss (Figs. 12 (c) and (d)). However, in the range from about 270 to 450 °CA, a reserve trend of the MOFT can be found, and the asperity contacts will occur in a short period in the cases with the speeds of 1000 and 2000 rpm. Although the direct contact among micro-asperities will dramatically increase the friction force and friction power loss, the higher average lubricant shear stress induced by the higher relative velocity as stated in Eq. (31) is still the main source of the friction. Therefore, the friction forces as well as the friction power loss still increase with speeds during this period. to reduce, and the dynamic behavior of the LBRBE is more regular and stable. This also means that the higher running speed will result in the decline of the dynamics performance of the LBRBE. In the horizontal direction, both the thrust forces acting on the piston ( Fig. 12 (a)) and the reaction forces at the LBRBE are also used to balance the horizontal inertia forces of the piston and connecting rod ( Fig. 13 (a)). It can be found that the horizontal reaction forces increase with the speeds.
In the vertical direction, it is obvious that the plot of the reaction forces has the same tendency with that of combustion pressures when the system has a small speed. However, with the increase of the speed, the tendency of the reaction forces comes to resemble that of the inertial forces ( Fig. 13 (b)), especially in the case with the speed of 6000rpm. Therefore, it can be concluded that the reaction forces at the LBRBE result from the joint effects of both combustion and inertial forces. In the case with the higher crank speed, the higher reaction forces at the LBRBE will be generated for most of the time (Figs. 15 (a) and (b)). Thus, the smaller MOFTs will be required so as to produce higher oil film pressure to separate the journal from the bushing. While during the period from about 270 to 450 °CA, the MOFT will depend on the joint effects of both combustion forces and inertial forces rather than decrease with the increase of the speeds.  The main reason may be that the higher speed will generate higher lubricant shear stress and thus higher friction force and friction power loss. Although in other moments, the cases with the lower speeds will have the thinner oil film thickness and thus more asperity contacts (Fig. 15 (c)), the lubricant shear stress induced by the relative velocity as stated in Eq. (31) is still the main source of the friction. Therefore, the friction forces as well as the friction power loss still increase with speeds during this period.

Conclusions
This study develops a new numerical model for the dynamics and lubrication analyses of a multibody system with multiple lubricated clearance joints. This method is established by coupling the dynamics model of the multibody system with the mixed lubrication models of translational and revolute clearance joints. The multibody dynamics model is established by Lagrange's method, while the lubricated clearance joints are built based on the average Reynolds equation and solved using the finite element method. The hydrodynamic forces built up by the lubricant fluid are evaluated from the motion state of the system and then input into the dynamics equations of the multibody systems.
The proposed approach is applied to a four-stroke gasoline engine with the lubricated piston-liner subsystem and bearing at the big end of the connecting rod. Both the dynamics and lubrication performances are also analyzed. From the results obtained, it can be found that there are obvious interactions between the piston skirt-liner system and the LBRBE. In the end, the performances of the system with different crank speeds are also analyzed. With the increase of the speeds, the dynamics and lubrication performances of the system will decrease.