Multi-body Dynamic Modelling and Error Analysis of Planar Flexible Multilink Mechanism Considering Clearance and Thermal-mechanical Coupling Effect of the Crankshaft-bearing Structure

: In order to study the dynamic position accuracy of bottom dead point (BDP) for multilink high- speed precision presses (MHSPPs), it’s essential to develop a dynamic model of planar multilink mechanism with clearance and spindle-bearing structure. Traditional models always neglect the effect of thermal characteristics of spindle-bearing structure, which reduces the prediction accuracy of dynamic model for multilink transmission mechanisms. To overcome the shortcomings of the previous models, a thermal network model (TNM) of the crankshaft-bearing system is established firstly considering the effects of thermal contact resistance and variable stiffness of bearing concerning the temperature rise. Then, dynamic model of the crankshaft-bearing system is built through the finite element method, which includes rigid disk, Timoshenko beam and quasi-statics model of ACBB. On this basis, an improved dynamic model of planar flexible multilink mechanism with clearance considering the thermal-mechanical coupling effect of the crankshaft-bearing structure is developed and the corresponding dynamic error dimension chain between slider and crankshaft is constructed in this work. Compared to the simulation from traditional models, the simulated slider’s BPD position error from the improved model agrees better with experimental data, which verifies the correctness of the proposed model. It’s demonstrated that the punching force and thermally induced variable stiffness of bearing lead to a significant increase of slider’s BDP pos ition error, which reduces the machining precision of MHSPP. Furthermore, the influence of crankshaft speed, contact angle of bearing and clearance size on the slider’s BDP position error is also investigated.

When the press operates under high-speed punching conditions, flexible spindle structure and variable stiffness of bearing cause vibration error between ideal and actual axis trajectory, which affects the machining accuracy and its stability seriously. Meanwhile, temperature rise and thermal expansion due to heat generation of bearings lead to reduction of relative position accuracy between upper and lower molds. Therefore, it's necessary to develop the dynamic model of planar flexible multilink mechanism considering thermal-mechanical coupling effect of rotor-bearing structure to study the dynamic error of bottom dead point (BDP) position for multilink high-speed precision presses (MHSPPs).
To study the transmission error of multilink mechanism, it's key to model the clearance of revolute joints and flexibility of components. According to the relative motion state between bushing and pin, models of revolute joints can be classified as one-state (continuous contact), twostate (contact-separation) and three-state (collision-contact-separation), which have been widely used in the dynamic analysis of mechanism [1][2][3][4][5][6]. To simplify the numerical calculation, a continuous contact force model combining Hertzian elastic contact component with energy dissipation is presented to describe the collision/contact state uniformly [7], such as Lankarani- Nikravesh model [8], Flores model [9] and Hu-Guo model [10]. It's noteworthy that one-state model is simple but with largest error, inversely the most complex three-state model has the highest accuracy. Furthermore, the elastic foundation model (EFM) [11] and finite element model (FEM) [12] have also been proposed to improve the predictive precision of contact pressure between two joint elements.
The presence of friction phenomena during motion between bushing and pin affects the dynamic characteristics of mechanism significantly. To simulate the phenomena of sliding friction, viscous slip, viscous friction and Stribeck between pin and bushing in the revolute clearance joint, static and dynamic friction models are formed [13]. Static friction models are mainly used to describe the steady-state behaviour of friction force, while dynamic friction models are applied to capture the dynamic characteristics using additional state variables. As a typical representative of the static friction model, the Coulomb friction model has such problems as zero velocity discontinuity and strong nonlinearity, and needs to be modified to obtain stable numerical results [14,15]. In order to overcome the shortcomings of static friction models in capturing friction phenomena such as pre-sliding displacement or friction hysteresis, dynamic friction models based on "state variables", such as Dahl model [16], bristle model [17] and LuGre model [18], etc, are developed to study the dynamic characteristics of mechanism.
Compared with clearance of revolute joints, elastic deformation of components has more significant impact on transmission accuracy of mechanism. FEM is mainly used to model the flexible components in order to study the dynamic characteristics of mechanism. Based on FEM and absolute node coordinate method (ANCF), Wang et al. [19,20] studied the dynamic response of flexible slider-crank mechanism with clearance and proposed an improved delayed feedback control strategy to stabilize its chaotic motion. Wang and Liu established a dynamic model of flexible five-bar linkage mechanism with clearance, and studied the influence of flexibility and clearance on its dynamic performance [21]. The authors also established the dynamic models of flexible slider-crank and multilink mechanism with clearance to investigate the influence of clearance and flexibility on the working accuracy of mechanism [22][23][24].
The existing methods for modelling and dynamic analysis of the spindle-bearing structure for machine tools are composed of the transfer matrix method (TMM) and FEM. The matrix order based on TMM won't rise with the increase of system's degree of freedom, which has been applied to model the spindle-bearing system. Aleyaasin et al. [25] established the flexural vibration model of a rotor-bearing system and deduced the general formula of system's tridiagonal division matrix through TMM to determine the unreasonable factors of model. Based on the concept of continuous system, Hsieh et al. [26] analyzed the coupled lateral and torsional vibrations of a symmetrical rotor-bearing system using TMM. Jiang and Zheng built a dynamic model of high-speed electric spindle system and analyzed the critical speed and dynamic stiffness characteristics of system [27,28]. When TMM is used to compute high-order complex rotor systems, it's easy to cause the problems of numerical instability and "leakage roots". One solution is proposed to overcome the shortcomings of TMM above, which uses FEM to deal with complex rotor coupling systems considering the effect of surrounding structure. Ruhl [29] developed a FM model of spindle rotor system with translational inertia and bending stiffness of bearings considered to predict its stability and unbalanced responses, which neglects the effect of rotational inertia, gyroscopic moment, shear deformation, axial load, axial torque, and internal damping.
Based on five-degree-of-freedom model of bearings and rigid shaft model, Aini [30] analyzed the dynamic characteristics of the spindle system supported by a pair of angular contact ball bearings (ACBBs). Considering the effect of radial clearance of rolling bearings and ball-race contact, Sinou [31] developed a FM model of flexible spindle supported by ball bearings and analyzed its nonlinear unbalanced response and rotor trajectory. When bearing is regarded as linear spring, Nelson [32] developed a dynamics model of the spindle system for machine tools, which considers the effect of rotational inertia, gyroscopic moment, shear deformation and axial load. Wang et al.
[33] established a vibration model of the spindle-bearing structure for crank press system and analyzed its vibration characteristics.
In the process of operating the spindle system, the thermal error caused by the heating of the bearing and other components cannot be ignored. Thermal analysis methods of the spindle rotorbearing system can be currently divided into two types: thermal network method (TNM) and FEM.
TNM has the advantages of small computation and is widely applied to predict the temperature field distribution of the spindle system. With lubrication and assembly form of bearings The existing models only consider the factors of revolute clearance joints, linkage flexibility and spindle-bearing structure dynamic, neglecting the effect of thermal characteristics of rotorbearing structure, which reduces the prediction accuracy of dynamic model for multilink mechanism of MHSPP. In this work, an improved dynamic model of planar flexible multilink mechanism considering clearance and thermal-mechanical coupling effect of rotor-bearing structure is developed and the corresponding dynamic error dimension chain between slider and crankshaft is constructed. The correctness of the proposed improved model is verified through comparison between simulation and experiment. Furthermore, the influence of crankshaft speed, contact angle of bearing and clearance size on slider's BDP position error is also investigated.

Traditional model of multilink mechanism with crankshaft-bearing structure
The solid model of the multilink mechanism for MHSPP is shown in Fig. 1. The mechanism is composed of the crankshaft-bearing structure and toggle mechanism in series. When crank shaft 2 rotates around point O, slider 4 is driven to execute reciprocating motion from left to right along the guide rail through connecting rod 3. Pendulum rods 5 and 6 are hinged with slider 4, which forces lower and upper sliders (7,8) to move back and forth in the vertical direction.
The dynamic accuracy of the MHSPP is mainly determined by the BDP position accuracy of lower slider 7. In order to predict and evaluate the performance of MHSPP, a dynamic model of planar flexible multilink mechanism with clearance and crankshaft-bearing structure is established, as shown in Fig. 2. The deformation of pendulum rod 6 due to the punching force between upper and lower moulds, gravity and inertia force of slider 7 is very significant and the flexibility of pendulum rod 6 and imperfect revolute joint between slider 7 and pendulum rod 6 are considered in this work. Furthermore, the effect of clearance between crank shaft 2 and connecting rod 3, and flexibility of crank shaft is also considered.

Model of revolute clearance joint
The centre of the pin coincides with that of the bushing in the ideal joint and only relative rotation occurs between the bushing and pin. In the actual revolute joint, there exists clearance between the pin and the bushing, which causes the centre deviation between two joint elements.
When the pin keeps contact with the bushing, it produces impact and the corresponding contact force. The model of revolute clearance joint is shown in Fig. 3 and the eccentricity vector that connects the pin and bushing centre can be given by Where m R represents the radius of the pin or bushing, iP RR = , jB RR = .
When the pin collides with the bushing, the penetration depth can be written as Where c is the radius difference between the pin and bushing, The relative contact velocity between the pin and bushing can be defined as Where n v and t v are the normal and tangential velocities of the contact point, respectively.
The angle between the eccentricity and the horizontal directions can be given by When the penetration depth is less than zero, contact force will disappear between the pin and bushing. On the contrary, contact force generates between two joint elements. According to the contact force model proposed by Lankarani and Nikravesh, the contact force can be divided into elastic and dissipative components [8]. Based on the Hertz contact theory, the contact force model combining with the hysteresis damping function can be defined as The contact stiffness coefficient n K can be expressed as The damping coefficient n D can be given by Where e c is the collision recovery coefficient, ()  − is the initial collision velocity before collision.
To overcome the numerical difficulties due to the saltation of friction force that results from the sudden change of tangential velocity, the Coulomb friction force is modified as [14,15] Where f c is the friction coefficient, d c is the dynamic correction factor.
The dynamic correction factor can be expressed as Where u v and l v represent the given upper and lower limits of the relative tangential velocity between the pin and bushing.
The multi-body dynamic equation of rigid multilink mechanism can be given by Where q represents the configuration coordinate, r M is the mass matrix with respect to the configuration coordinate, r F is the external force matrix, and λ is the Lagrange multiplier,

Model of flexible linkage
When the MHSPP operates, especially under punching condition, the influence of linkage's deformation on the dynamic accuracy of slider's BDP position cannot be ignored. The original rigid linkage must be described as a flexible model, as shown in Fig. 4. In this work, ANCF is used to establish the model of flexible linkage 6 and its dynamic equations can be given by The multi-body dynamic equations of the planar flexible multilink mechanism with clearance can be expressed as Where   =   

Model of the crankshaft-bearing system
The structure of the crankshaft-bearing system for MHSPP is shown in Fig K is the stiffness matrix due to the axial load,{} b F is the discrete and concentrated force vector.

(3) Model of the ACBB
The stiffness of ACBB has a significant impact on the dynamic responses of the crankshaftbearing system for MHSPP. It's necessary to establish a mathematical model of the ACBB to calculate its stiffness. As shown in Fig. 7 Where N is the number of rolling element for the ACBB.
The curvature centres of inner and outer raceways are initially coincided with the centre of rolling element. As shown in Fig. 8, the distance between curvature centres of inner and outer raceways can be calculated as Where D is the diameter of the rolling element.
When the press system operates, the external load acting on the ACBB will cause its contact angle to change. It is assumed that the outer raceway of ACBB is fixed and the distances between the curvature centre of the inner ring groove and the final position of the ball centre, between the curvature centre of the outer ring groove and the final position of the ball centre can be given by Where ik  and ok  represent the deformation between inner raceway and the rolling element, between outer raceway and the rolling element. Based on Pythagorean theorem, the position geometric equations of ACBB can be given by  and ok  represent the contact angles of the inner and outer raceways, ik Q and ok Q are the contact forces of the inner and outer Ek can be calculated as Based on Hertz theorem, the contact forces ik Q and ok Q can be given by Where i K and o K are the contact stiffness of the inner and outer raceways, which can be written as Where 0.636 When the inner raceway keeps in contact with the rolling element, x R and y R can be given by When the outer raceway contacts with the rolling element, x R and y R can be calculated as It is assumed that the initial variables and the corresponding errors of Eqs.18 and 19 are respectively. The position geometry and the force balance equations can be rewritten as To obtain four parameters k U , k V , ik  and ok  , Newton-Raphson method is used for iteration.
Where [] ij a represents the matrix element, / ij   .
Based on the contact force between the rolling element and inner raceway, the external load acting on the inner raceway can be given by  , the equivalent stiffness matrix of the ACBB can be obtained as Combining the rigid disc, Timoshenko elastic beam and ACBB models, the nonlinear dynamic equations of the crankshaft-bearing system can be expressed as T is the period of the pulse function, is a variable related to the pulse width, max F is the amplitude of the blanking force. The dynamic simulation parameters of the crankshaft-bearing system are listed in Table 1.

Improved model of multilink mechanism with the crankshaft-bearing structure
When the press is working, heat generation produces between rolling elements and inner/outer raceways as a result of friction. Due to the influence of heat conduction and convection, the temperature rise of the bearings and adjacent components increases significantly and thermal deformation brings forth. The axial and radial displacements caused by the thermal expansion of ACBB affect the stiffness of the bearings and vibration response characteristics of the crankshaft-bearing system. In order to predict the position error of BDP for MHSPP accurately, it's essential to develop a dynamic model of planar flexible multilink mechanism with clearance considering the thermal-mechanical coupling effect of the crankshaft-bearing structure.

TNM of the crankshaft-bearing system
To simplify the complexity of the model, only variable bearing stiffness as a result of thermal expansion is considered in this work, while the effect of thermal deformation of crankshaft on the dynamic characteristics of crankshaft-bearing system is neglected. According to the physical structure of the crankshaft-bearing system, it is discretized into 70 nodes and the corresponding TNM is developed, as shown in Fig. 13.
The temperature on the centroid of each node is regarded as the average temperature of the node. Heat transfer of each node is carried out by means of heat conduction and thermal contact resistance between two adjacent contact nodes exists. Nodes that are contact with environment directly transfer heat with air through convection and there exits convection heat transfer resistance. The energy transfer relationship among nodes in the TNM is shown in Fig. 14. Based on the principle of energy conservation, the finite difference (FD) equation is used to describe the temperature field of nodes.
Where i T is the temperature of node i , n T is the temperature of node n ,which is adjacent to in R is the thermal contact resistance between nodes i and n , i q is the friction heat at the contact zone of rolling element, i C is the thermal capacity of node i , i m is the mass of node i .

(1) Heat generation of bearing
The heat of the crankshaft-bearing system mainly results from the friction between rolling elements and inner raceway, rolling elements and outer raceway. The total heat generated from the bearing due to friction torque can be calculated as Where M is the total friction torque of the bearing,  (2) Thermal contact resistance Thermal contact resistance between the rolling elements and outer/inner raceways for the ACBB can be written as Where b k represents the thermal conductivity between the balls and outer/inner raceways, Thermal contact resistance between inner raceway and crank shaft can be given by Where c h is the thermal contact conductance between inner raceway and crank shaft, which can be determined as [46]  Thermal convection occurs between the crankshaft-bearing system and surrounding air, including forced convection between rotation bodies and air, and free convection between stationary surfaces and air. The convective heat transfer coefficient can be calculated as The Nusselt number of free convection between the bearing support and air can be given by

(3) Model verification and temperature field analysis
The parameters of TNM for the crankshaft-bearing system are listed in Table 2. When the rotation speed of MHSPP is 200 rpm with the initial temperature of the surrounding air 25 0 C, the simulated and experimental temperature rise of four bearings for the crankshaft-bearing system under no-load condition is shown in Fig. 15.  It can be seen from Fig. 15 that the simulated temperature rises of the inner and outer raceways for four bearings are generally consistent with the experimental data, which proves the correctness of TNM for the crankshaft-bearing system. It's demonstrated that the temperature rises of four bearings tend to be stable at about 90 min under no-load condition, and the stable temperature of the inner raceway is larger than that of the outer raceway.
The steady-state temperature distribution and temperature rise variation of the network nodes for the crankshaft-bearing system are shown in Figs. 16

Thermal-mechanical coupling model of the crankshaft-bearing system
Positioning mode is adopted to perform the preload of the bearing in the crankshaft-bearing structure of MHSPP, as shown in Fig. 18. The relative axial displacement between inner and outer raceways under static condition is caused by preload acting on the bearings, which can be divided into two stages. In the first stage, the axial clearance 0.5 e P is eliminated and a free contact angle 0  is produced. Then, a new axial displacement p  causes the internal contact deformation during the second stage. For the positioning preload mode, the axial displacement between the inner and outer raceways is limited and the total displacement  during two stages remains constant, 0.5 ep P  =+ .

(1) Radial thermal displacement of inner ring raceway
The radial thermal displacement bi u of the inner ring raceway considering the effect of its temperature rise for the bearing can be calculated as The radial thermal displacement of the inner ring raceway due to the thermal expansion of the crank shaft can be written as The radial thermal displacement r u between the inner and outer ring grooves can be given by Temperature variation of bearings will lead to the geometric dimension changes of rolling elements and inner/outer raceways. It is assumed that the curvature centre position of outer raceway is fixed and the geometrical relationship of ACBBs considering the thermal-mechanical coupling effect is shown in Fig. 19. The distances between the final curvature centre of the inner ring groove and ball centre, between the final curvature centre of the outer ring groove and ball centre can be given by .19 Geometric relationship of ACBBs considering the thermal-mechanical coupling effect The geometric equations 27(a) and 27(b) can be updated as Where ik   and ok   are the contact angles between the rolling element and inner/outer raceways considering the effect of thermal expansion, respectively. These can be calculated as Similarly, the equivalent stiffness matrix of the ACBB considering the effect of the thermal expansion can also be obtained through Eq. (30). When the crankshaft speed is 200 rpm, the preload and radial stiffness of the ACBB considering the effect of the thermal expansion are shown in Fig. 20.  stiffness of ACBB reduces from 3.33×10 7 to 3.25×10 7 N/m after 4 min due to the thermalmechanical coupling effect of bearing, as shown in Fig. 20(b). It's demonstrated that the preload rises significantly, while the radial stiffness decreases due to the bearing's thermal expansion.
The dynamic simulation flow chart of the crankshaft-bearing system considering the thermalmechanical coupling effect is shown in Fig. 21.

Model validation and results discussion
The position accuracy of slider's BDP for MHSPP is mainly affected by the transmission errors of multilink mechanism and the crankshaft-bearing structure. As shown in Fig. 2, the dimension chain of the flexible planar multilink mechanism with clearance and crankshaft-bearing structure can be given by  In order to verify the effectiveness of the improved model, the dynamic responses of slider for the multilink mechanism are measured, based on the experimental platform of multilink mechanism shown in Fig. 24(a). The dynamic test system of slider for the multilink mechanism, is mainly composed of position magnet, L-shaped fixed frame, waveguide and electronic chamber (strain pulse detector). Figure 24 The parameters of flexible multilink mechanism with clearance are listed in Table 3   It's inferred that the stable time of slider's BDP position error is not sensitive to the clearance size.

Conclusion
(1) TNM of the crankshaft-bearing system is established considering the effects of thermal contact resistance and variable stiffness of bearing concerning temperature rise and the effectiveness of the developed TNM is verified through experiment. The simulation results show that the maximum temperature rise of the crankshaft-bearing system is concentrated on the nodes near the rolling elements of bearing and the corresponding average temperature rise reaches 25 0 C, while the average temperature rise of the remaining nodes is only 5 0 C.
(2) Thermal-mechanical coupling model of the crankshaft-bearing system is built. The preload of ACBB is increased significantly, while the radial stiffness is reduced as a result of bearing's thermal expansion. The vibration amplitude of each shaft section under blanking condition is larger than that under no-load condition and the zigzag shapes of the shaft centre trajectories are obvious due to the thermal expansion of ACBBs and the eccentric mass of the synchronous pulley.

Conflicts of Interest
The authors declare that they have no known competing financial interests or personal relationships that could have appeared in this work. Table 1 Dynamic simulation parameters of the crankshaft-bearing system Table 2 Parameters of TNM for the crankshaft-bearing system Table 3 Parameters of flexible multilink mechanism with clearance The list of figure captions: