A Review of Dynamic-Tribological Simulation Methods for Sliding Bearings in Internal Combustion Engines

The problem of friction reduction and wear resistance of sliding bearings is one of the key factors in determining the overall performance of internal combustion engines. This paper investigated and summarized the theoretical and simulation models of multi-body dynamics of crankshaft system, tribology of sliding bearings, and the wear calculation methods of the shaft-bearing friction pairs. Existing studies show that the dynamics model, hybrid lubrication model, and the friction and wear models request to be upgraded by comprehensively considering the material, structure, manufacturing process, working conditions, and etc. Based on the research status and existing problems of the above analyses, this paper summarizes the simulation models applicable to the field of dynamics and tribology of sliding bearings and presents the prospects for optimization of wear calculation methods for sliding bearings.


Background
As a power system with compact structure, wide power range and high adaptability, the internal combustion engine (ICE) is widely used in engineering and agricultural machinery, automobiles, motorcycles, national defense and other fields [1][2]. With the further development of the economy and national defense construction, the technical level and comprehensive performance of the engine put forward more comprehensive and strict requirements [3][4]. High power density, high compact structure, low fuel consumption, low emissions and high reliability become the common features and goals of the nest generation of engine design and development [5][6].
The crankshaft system is a mechanical system that converts the firing pressure of combustion gases inside the cylinder into cranking torque. Its power response, friction performance, and fatigue life of relevant components directly affect the ICEs' power transmission and output functions, which are also the key to the durability of the whole engine [7]. The sliding pairs, such as crankshaft and main bearing, crank pin and connecting rod big end bearing are important driving mechanism in the crankshaft system. Their working condition and service life are critical to the reliability of crankshaft system.
With the further increase of engine speed and power density, the working conditions of the sliding bearings in the crankshaft system become harsh and variable, such as high temperature, heavy load, impact, and corrosion, all of which have serious impact on the bearings during their service. As a result, the lubricating condition of the journal-bearing friction pair deteriorates with the contact and collision of the relative motion of two components, which intensifies the wear of the bearing and the dissipation of the friction force. The clearance of friction pair will increase simultaneously during the wear process, affecting the dynamic performance of the friction pair [8].
In summary, the friction-reduction and anti-wear problem of sliding bearings is of great importance to optimize the dynamic response and fatigue properties of crankshaft, connecting rod, piston and other components, and to improve the durability of the ICE mechanical system. Consequently, the mechanism research and simulation analysis of the dynamics, lubrication and wear performance of engine sliding bearings become the research interests.
The primary work of this paper is to summarize the relevant theories and research status in the fields of multi-body dynamics of crankshaft system and tribology of engine sliding bearings, with the aim to lay the theoretical foundation for the establishment of a more comprehensive and accurate bearing wear calculation model. In addition, combined with the existing wear calculation method of journal-bearing pair, the numerical model and simulation method of engine sliding bearings are also analyzed and summarized, providing a basis and reference for the application of wear simulation methods in engineering.

Multi-Body Dynamics of the Engine Crankshaft System
Multi-body dynamics is a field that studies of the relative laws between the kinematics and dynamics of multi-body system, including both multi-rigid-body system and multi-flexible-body systems. The multi rigid body model is difficult to reflect the large range of nonlinear motion and elastic deformation of the components, due to the complex structure of crankshaft, bearing, body and other components, and also due to the nonlinear and high load operating conditions. As a result, it is necessary to introduce multi flexible body model to analyze the dynamics of crankshaft system [9].
Since the 1960s, the multi-body dynamics analysis commonly assumed that the subject structure is a rigid body for multi-body system modeling. Since the above methods still cannot accurately describe the mechanical properties of flexible components in complex multi-body systems, domestic and abroad scholars conducted comprehensive and pervasive discussions and studies on the dynamic modeling theory and numerical algorithms of multi-flexible-body system in recent years. The dynamics theory of rigid-flexible coupling systems has been gradually established, which can not only accurately describe the coupling dynamic stiffness of the flexible components, but also reflect the influence of the coupling relation on the damping characteristics of the system. At present, finite element method (FEM) [10], Craig-Bampton modal synthesis method [11], and the substructure method [12] are commonly used in modeling the kinematics and dynamics of flexible bodies and couplings in the crankshaft system.
The establishment of dynamical equations for multi-flexible-body systems is far more complicated than that for multi-rigid-body system. Considering the huge computing scale and arduous calculation work, both domestic and foreign scholars proposed the algorithms that helped to simplify the multibody dynamics equations.
Cao et al. [13] from Tsinghua University proposed an integrated solution strategy that combines the OpenMP concurrent computing system and the parallel sparse linear equation solver Pardiso. Wu et al. [14] from Dalian University of Technology proposed an alternating conservative perturbation iterative algorithm for the coupling function between rigid and flexible bodies in the dynamical equations to handle the low and high frequency vibrations of the system separately, thus increasing the time increment. Guo et al. [15] from Nanjing University of Science and Technology compared several typical numerical integration methods for flexible multi-body dynamics and concluded that the computational efficiency of the display algorithm is higher than that of the implicit algorithm, but increased the dependence on the time increment. Fang et al. [16] from Shanghai Jiao Tong University adopted a linear implicit multi-step based high-order A-stable MEBDF integration algorithm for the time integration of the multi-body dynamics equations, which greatly improved the overall computational efficiency. The engine crankshaft system consists of piston assembly, connecting rod assembly, and crankshaft flywheel set, including the piston, the cylinder liner, the connecting rod, the crankshaft, the bearings and other components. In multi-body dynamic simulations of crankshaft systems, the piston and cylinder liner are usually simplified as a set of rigid nodes. The master nodes and master degrees of freedom (MDOF) of components such as the connecting rod, the crankshaft and the engine body are extracted by substructure condensation to retain the mass and linear stiffness of each component. Considering the output forces and torques of the crankshaft under complex multi-axial loads, which need to be handled flexibly, a modal neutral file (MNF) is created by the simulation software, which contains information about the materials, nodes, elements and modes. Multi-body dynamics and finite element theory are applied to establish the non-linear connections between components, such as the oil film and contact between the shaft and the bearing. The internal structure of crankshaft system and the relative movement of each component are simulated efficiently and precisely for the subsequent dynamic simulation of the components and the obtainment of stresses, displacements and velocities [17].

Tribological Simulation of Engine Sliding Bearings
Tribology is the science and engineering of studying interacting surfaces in relative motion, including the investigation and application of the principles of friction, lubrication and wear. The ICE sliding bearing mainly includes the main bearing, connecting rod big-end bearing, connecting rod small-end bearing, and end thrust bearing. As the main moving parts of engine sliding pairs, their tribological properties provide the boundary conditions for the multi-body dynamics simulation of crankshaft systems.
In general, it is assumed that the bearings are remained in a hydrodynamic lubricated condition. However, variations in speed and firing pressure may lead to contact and collision in certain areas under operative conditions. During the star-stop process, the oil film thickness may be reduced to the same order of magnitude as the roughness, leaving the bearing in a state of mixed lubrication, boundary lubrication, or even dry friction.

Lubrication Theory and Model of ICE Bearing
Based on the fluid differential unit, Reynolds equation describes the mechanism of hydrodynamic pressure and becomes the basis of modern hydrodynamic lubrication theory. Due to the lack of computing devices in the early days, the research related to fluid lubrication mainly focused on the simplified calculation method of Reynolds equation.
For the radial sliding bearing, considering its axial constraints and the incompressibility of the fluid, the simplified form of Reynolds equation is shown in equation (1). 3 3 Where h and p are the thickness and pressure of oil film, respectively,  is the dynamic viscosity of lubricant, and U is the component of fluid velocity along the x direction [18].
With the increasingly development of computer simulation technology and numerical calculation methods, the simulation work on hydrodynamic pressure lubrication of bearings has developed rapidly. The thermal effect, elastic deformation of bearing and other important factors have been included in the modification version of the Reynolds equation and boundary conditions, and the lubrication model of elasto-hydrodynamic (EHD) considering the elastic deformation of the bearing, thermohydrodynamic (THD) considering the thermal effect and thermo-elastic-hydrodynamic (TEHD) considering multiple factors of the crankshaft system and bearing.
The extended Reynolds equation of elasto-hydrodynamic lubrication is shown in equation (2), which considered the influences of dynamic pressure effect, expansion effect and extrusion effect of oil film in the sliding direction, which is able to conduct a more accurate analysis of the bearings under various lubricating state [19].
Where  is the angular position of z-axis (radial direction) and R is the bearing radius, u is the surface velocity ( In the simulation of thermo-hydrodynamic lubrication, it is assumed that the dynamic viscosity of lubricant stays the same along the direction of oil film thickness and that the density of lubricant ignores the temperature. Then the two-dimensional adiabatic energy equation commonly used in thermo-hydrodynamic lubrication can be obtained, as shown in equation (3), also known as Cope theory [20].
Where x q and z q are the quantity of flow in x and z directions, respectively, T is the temperature,  is the density of lubricant, and c is the radius clearance of bearing.
The Cope theory assumes that the viscous friction heat is all dissipated by the lubricant, and that the friction heat conduction in the direction of oil film thickness is ignored. Hence, large differences exist between the calculations and test results.
In the year of 1962, Dowson et al.
[21] established a generalized Reynolds equation considering the variation of viscosity in the radial, axial and circumferential directions, as shown in equation (4).
Where U and V are the radial and axial linear velocities of the surface, respectively. In recent years, the factors of surface roughness and texture, the manufacturing and processing technology of bearings, the journal misalignment, the assembly error of journal-bearing, the rheological properties of lubricant, the cavitation effect and the elastic and plastic deformation of bearing materials, and their influence on the lubricating properties of bearings have gained more attentions in the field of engine sliding bearing research.
Among them, Li et al. [22] from Hefei University of Technology considered the axial movement of crankshaft in the actual engine working condition. The bearing load, crankshaft deformation, and journal inclination under different working conditions were calculated and the influence law of axial movement of crankshaft on the lubricating property of bearings were studied. Zhou et al. [23] from Beijing Institute of Technology studied various nonlinear factors (oil pressure, structure and dimensions of oil groove, lubricant quality, surface roughness, and etc.) in the service process of crankshaft and the law of influence on the dynamic property of crankshaft, deformation coordination and lubricating properties of bearings. According to the basic theory of elasto-hydrodynamic lubrication, the influence mechanism of temperature and geometric form error were also analyzed. Zhang et al. [24] from Shanghai University studied the influence mechanism of non-Newtonian properties of lubricant and surface topography effects on the lubricating properties of bearings under dynamic load. Liu et al. [25] from Shandong University established a thermo-elasto-hydrodynamic mixed lubrication model of the connecting rod large end bearing considering the effect of cavitation and microscopic elasto-hydrodynamic lubrication, proposed a method of recognition of cavitation position, and obtained the law of heat distribution of bearing friction loss.

Friction Theory and Modeling of the Engine Bearing
For ICE sliding pairs, the friction model is generally used to calculate the tangential forces on the component surfaces in contact and collision states determined by the lubrication boundary conditions. Coulomb's law of friction, one of the classical friction theories, defines that the friction between two surfaces is proportional to the normal pressure in the direction opposite to the direction of relative motion.

Static Friction
Mode. The Coulomb model is a static friction model based on Coulomb law of friction. The relationship between friction force f and relative sliding speed v is shown in equation (5).
Where sgn (v) is a symbolic function, representing the mutation in the direction of friction force f [26].
The Coulomb model simply describes the sliding friction, but ignoring the deviation between the maximum static friction and the sliding friction. The tribological properties of the "zero velocity zone" shown in figure 1 are also ignored, resulting in a discontinuous description of the friction state transition. Therefore, it is barely used in mechanical systems where adhesion and inversion occur.
Considering that the maximum static friction is larger than the sliding friction, Doebelin, Li Chun Bo et al. incorporated the effect of static friction into the Coulomb model to simulate the nonlinear characteristics of friction in the "zero speed zone".
Doebelin et al. [27] used a straight line with a certain slope to display the friction characteristics of the "zero speed zone", with the slope determined by the velocity threshold and Coulomb friction. The approximation method is not only convenient for modeling but also has a relatively simple algorithm, and the solving accuracy of the model can be adjusted by the variation of V D . However, when the external force is less than the maximum static friction, the object will acquire a certain acceleration according to the algorithm of this model, which is different from the actual phenomenon.
Li et al. [28] considered the influence of negative viscous damping on the friction variation in the low-speed region and used an exponential model to describe the Stribeck effect, as shown in equation (6).
Where s f is the maximum static friction, s v is the Stribeck velocity,  is the empirical constant.
The Stribeck friction model, as shown in figure 2 and equation (7), can be obtained when the Stribeck effect is included in the Coulomb model. f   In the process of simulation based on the static friction model, it is difficult to estimate the relative motion of the two colliders because the friction in the "zero velocity region" is difficult to solve, resulting in distortion when calculating the nodal values of the distributed contact forces. Hence the dynamic friction model is needed to detect the relative sliding and stationary states between the pairing elements to avoid energy increment caused by collision during the simulation.

Dynamic Friction
Model. The static and dynamic friction models can be distinguished according to whether they can be described by differential equations. The former describes friction as a function of relative velocity, while the latter describes friction as a function of relative velocity and displacement [29]. In the dynamic friction model, the interface is deformed in both the normal and tangential directions, both in the stationary phase and in the sliding phase. The most common dynamic friction models are the Dahl model and the LuGre friction model. The Dahl model considered the small displacement (pre-sliding displacement) and the hysteresis of the interfacial friction before the friction force reaches its maximum static friction value [30]. The presliding displacement has a certain amount of elasticity that allows the object to move in the opposite direction and reach its original position after the friction is removed, while retaining a certain residual displacement, as shown in figure 3 and figure 4. The relation between stress and strain on the friction surface is described by the partial differential equation as shown in equation (8). Where x is the deformation amount, v is the relative velocity,  is the stiffness coefficient,  is the shape coefficient. The curvature of curve increases with . The LuGre friction model is an extension of the Dahl model, to which the bristle model is added [31]. The rough contact surface is equivalent to a large number of elastic bristles with random behavior in the model. The average deformation of the bristles and the friction force caused by deflection are shown in equations (9) and (10). g Where 0  is the bristle stiffness, 1  is the microscopic damping coefficient, v is the relative sliding velocity between the two surfaces, and s v is the Stribeck velocity.
The LuGre friction model provides a comprehensive account of the Coulomb friction, pre-sliding, variable static friction, Stribeck effect and friction hysteresis. As a continuous model with a smooth transition between different friction states, it has been widely used to solve dynamic friction in various  is the major difficulty of this model [32].
In the field of tribology of sliding bearings, lubrication and friction models are the basis and prerequisite for wear simulation and prediction. In view of the high speed and heavy load conditions of the engine, thermal effects, elastic deformation of the bearing, surface roughness and other factors need to be considered when building a mixed lubrication model of the engine bearing.
At the same time, the relative movement and the lubricating state of the journal-bearing pairs need to be determined in order to select the corresponding lubrication and contact models and thus determine the boundary conditions for the multi-body dynamics model of the crankshaft system. When performing the coupled simulation of dynamics and tribology of the crankshaft system, the kinematics and dynamics parameters of each node on the bearing surface are required for the following wear calculations of the contact area.

Wear Theory of Engine Bearing.
Wear is a phenomenon of material loss caused by friction on the surface of contacting objects, which represents the inherent mechanical properties of the material, and is inevitable during engine operations. As a reflection of the tribological properties of the system, different degrees and types of wear would occur under different working conditions with different materials. The main parameters of wear property are wear quantity and wear rate, including linear wear quantity and rate, volume wear quantity and rate, and weight wear quantity and rate.
Typically, the friction pairs in the engine mechanical system go through three stages from start-up to failure, i.e., the run-in stage, the stable wear stage and the rapid wear stage. The three-stage wear pattern is generally expressed as a curve relating wear volume and rate to time. The curve of wear rate v.s. time is also known as the "bathtub curve", as shown in figure 5. Due to the roughness of the surface, the external load is borne by the randomly distributed asperities during the interaction process, and the contact pressure per unit will be relatively high. In view of some factors induced by fabrication technology and fitting process, the clearance between the surfaces is not uniform, resulting in unstable oil film. Therefore, in the run-in stage, the bearing is generally in a mixed lubrication state.
Toward to the end of the run-in stage, the asperities are gradually worn. Then the contact pressure per unit decreases with the increase of the contact area, the clearance tends to be uniform, and a stable oil film is established between the friction pairs. Finally, the lubricating state of the bearing gradually converts from irregular to the mixed state, entering the hydrodynamic lubrication and stable wear stage. The stable wear stage is the normal working stage of the bearing, with its length to be related to the working conditions and run-in quality of the interface in the run-in stage.
After a sufficiently long period of stable wear phase, the wear may be aggravated because of load fluctuations, lubricant film failure, fatigue damage of the surface material caused by long-term alternating stress, and etc. The severe wear stage is the accumulative effect of long-term wear and tends to be abrupt. Once it happens, both the curves of wear quantity and rate would increase sharply.
The main forms of bearing wear include adhesive wear and abrasive wear. Adhesive wear refers to the wear caused by the extrusion of the adhesive lubricant film and the instantaneous contact between the interface asperities. Under the relative sliding of the friction pairs and firing pressure load, plastic deformation will occur at the concave and convex parts of the contact surfaces. As the temperature increases, the metal material may melt, causing adhesion and welding-on between the interface, and the asperities may rupture under shear forces.
Abrasive wear refers to the wear phenomenon caused by hard particles and wear debris that cannot be discharged from the lubricating system. According to statistics, the loss caused by abrasive wear accounts for about half of all kinds of wear losses.
There are two kinds of abrasive wear: the one caused by brittle micro-peeling is called two-body abrasive wear; and the other caused by the contact between the asperities and particle shattered in the clearance of friction pair is called three-body abrasive wear. In accordance with different material properties, the failure modes include scratch, scuffing, ploughing and micro-cutting.
Engine bearings are in hydrodynamic lubrication most of the time, which are likely to enter boundary lubrication during starting, stopping, and switching. The asperities of the journal and bearing surface chip off due to adhesive effect and adhesive wear mainly occurs on the bearing surface. With the increase of the size and number of abrasive particles, the major wear form of bearing becomes abrasive wear.

Wear Modeling of Journal-Bearing Pair.
In 1946, Holm et al. [33] proposed a wear theory based on interatomic interactions of the contact surface between two solids at the microscopic level. If the distance m between the atoms on each surface is smaller than the distance d between the atoms inside the components of the friction pair, a large molecular force will be generated between the atoms on the interface, resulting in serious adhesive effect. Wear occurs when the adherent surfaces are pulled apart and a certain number of atoms will be removed.
Assuming that there are z atoms being removed for each contact and that the volume of atom is 3 d , the total wear volume can be calculated, as shown in equations (11) and (12).
Where L is the sliding distance, r A is the contact area, N F is the normal load, and z  is the contact stress. Soon after, Archard [34], a British scholar, proposed the adhesive wear theory and the corresponding model in 1953, which described a pair of hemispherical asperities with the same radius in contact with each other on the interface. Given that the probability of particle generation in a relative motion is uncertain, the asperities between the two surfaces are numerous, assuming that there are n hemispherical asperities between two surfaces, the probability of particles generation in a relative motion is K and the total wear quantity generated when the relative sliding distance between two surfaces is L, as shown in equation (13).  (13) In order to upgrade the theoretical system of wear, Archard et al. [34] extended the model to the application of abrasive wear, fatigue wear and corrosive wear in 1980. The general formulae of the model regarding the above wear mechanism are the same as for the adhesive wear model, in which only the definition of wear coefficient K varies according to the different wear mechanisms. In the abrasive wear model, the wear coefficient K is the product of the geometric parameters of asperity and the probability coefficient, as shown in equation (14).
In the late 1950s, based on the micro characteristics of wear, Raqinowicz [35] reintroduced the atomic forces into the wear model, considering the atomic properties on the surface and the mechanical properties of the material in the form of surface energy. Then, many scholars have successively committed to some new theories based on atomic properties. For example, Suh [36] and Vijh [37] proposed two microscopic wear models by considering the dislocation behavior and cohesive bond energy of atoms, respectively.

Calculation Method Based on the Archard Model
The theoretical model of adhesive wear proposed by Archard [34] has been cited and validated by many scholars in the field of tribology. The Archard wear model is a numerical model based on the mechanism of adhesive wear, which simulates the removal process of material caused by the relative sliding of the friction pair. It is widely used for the numerical calculation of wear of various mechanical parts. The original form of the model is the formula for calculating the amount of wear, as shown in equation (15).
The differential value of wear depth dh can be obtained by taking the derivative of both sides of equation (15) and dividing both sides by the contact area A, as shown in equation (16).
where, / N FA  is the contact pressure between the nodes, whose distribution along the sliding direction can be expressed as p(s), a function of the relative sliding distance s.
In the substructure condensation model of the bearing, the geometric parameters of each node contain a set of initial values of the radial position. The position will be in flux during the calculation process because of the dynamic nature of the system. The wear depth at each node will be accumulated simultaneously. Since the wear depth at each node determines the topography of the bearing surface and the values of the contact pressure and sliding distance change at each step, it is necessary to calculate the wear depth and wear profile repeatedly, as shown in equation (17). 10 node at the moment j, and j s  is the relative sliding distance between the two nodes on the interface from the moment j to j-1.
According to equation (17), the wear coefficient K, the dynamic contact pressure j p , and the relative sliding distance j s  are key parameters for calculating the wear depth h. The relative sliding distance and the distribution of dynamic contact pressure can be obtained from the simulation results of kinematics and dynamics analysis of the crankshaft system.

Methods for Wear Coefficient Determination and Modification
The wear coefficient is a characteristic of bearing materials that holds a variety of different definitions on different scales, depending on the object and the method. For example, it is generally used as a nodal parameter in finite element simulations and can be used as a characteristic parameter of the system in multi-body dynamics simulations and wear calculations. It can be affected by the physical properties of the material, surface roughness, working condition parameters and many other factors. In addition, the influence law of each factor is complex.
In engineering applications, the wear coefficient of material is commonly calculated with the wear quantity which can be measured in the frictional wear test of a standard specimen. The calculation formula is shown in equation (18) and (19).
Where h is the wear depth of the specimen, T and M are respectively the initial thickness and weight of the specimen, M is the total wear quantity, v is the sliding speed measured in the test, F is the applied load, and w is the wear rate of the specimen.
In engineering applications where accuracy is not highly demanded, the span of the wear coefficient can be roughly determined based on empirical values under different operating conditions. It can also be determined by the wear spectrum proposed by Czichos H [38], which includes a variety of wear coefficients for different lubrication conditions. The wear spectrum for different speeds and load conditions is shown in figure 6. Under light load conditions, the wear coefficient is mainly distributed in the lower region of the spectrum, while the wear coefficient under heavy load conditions is mainly located in the upper region. The applicable lubricating state only includes mixed lubrication and dry friction and the range of wear coefficient is between 10 -10 and 10 -2 mm 3 /Nm. Meanwhile, based on the theory of thermal desorption and load sharing, Sang Myung Chun et al. [39] modified the value by introducing the defect coefficient of oil film and the load proportionality coefficient in the calculation of wear coefficient. He obtained the wear coefficient under the mixed lubrication state as shown in equation (20) and (21).

Difficulties and Outlooks
The factors causing wear-out failure of engine bearings are complex because of the complicated structure and variable working conditions. It is necessary to comprehensively consider the influences of material, structure, manufacturing, working conditions and other factors on the friction and wear properties. Among them, in the dynamics and tribology coupling simulation, the factors affecting wear properties include speed, load, lubricating state, surface condition, and etc., which comprehensively determine the sliding speed, contact pressure and nominal clearance of the nodes. In the process of dynamics and tribology simulation of engine sliding bearing, the major difficulty is to consider the effects of the above factors when modifying the lubrication, friction and wear models. Since the calculation of dynamic wear requires continuous renewal of the wear profile and dynamic response, more targeted and efficient models and algorithms are also needed. Meanwhile, the judgment and conversion of oil film bearing and asperity contact under mixed lubrication state is worth further exploration.
In recent years, some scholars have optimized the wear coefficient in the local area of bearing surface by introducing fatigue model, damage model, thermal desorption theory, load sharing theory, and etc.
In general, the wear process of sliding bearings is a complex engineering problem involving multibody dynamics, tribology, thermodynamics, fatigue strength of components and parts, fatigue damage properties of materials, and other multi-scale and multi-disciplinary fields. Considering that the wear issue also involves microscopic parameters and their influence on the wear resistance of the material, the coordination between the dynamic variation of macro-structure and kinematic parameters is also in urgent need of further research.