State of the art of friction modelling at interfaces subjected to elastohydrodynamic lubrication (EHL)

Elastohydrodynamic lubrication (EHL) is a type of fluid-film lubrication where hydrodynamic behaviors at contact surfaces are affected by both elastic deformation of surfaces and lubricant viscosity. Modelling of contact interfaces under EHL is challenging due to high nonlinearity, complexity, and the multi-disciplinary nature. This paper aims to understand the state of the art of computational modelling of EHL by (1) examining the literature on modeling of contact surfaces under boundary and mixed lubricated conditions, (2) emphasizing the methods on the friction prediction occurring to contact surfaces, and (3) exploring the feasibility of using commercially available software tools (especially, Simulia/Abaqus) to predict the friction and wear at contact surfaces of objects with relative reciprocating motions.


Introduction
A machine transfers the motion and force from input to output component. Therefore, relative motions exist among parts or components in implementing motion transformations. The surface contact of two parts with a relative motion involves the friction force and wear over parts. The tribological behaviors of the parts such as bearings and seals affected the efficiency and reliability of machines [1,2], and the failure mode of excessive wear at sliding contact was one of the most pervasive surface failure modes of machine elements [3]. In addition, Elastohydrodynamic lubrication (EHL) is a type of fluid-film lubrication where hydrodynamic behaviors at contact surfaces are affected by both elastic deformations of surfaces and lubricant viscosity. The EHL lubrication conditions existed in many critical applications such as highly stress gears, bearings, cams, and soft-bearing elements including elastomeric bearings, seals, and synovial joints [4]. Furthermore, there are studies on using friction as an active role in the device, e.g., inertiafriction actuation principle [5−10]. There are studies on complex dynamics involving friction, e.g., in the medical device [11,12].
Friction causes wear and damage of two contacting parts with a relative motion. Wear refers to the material loss at the contact surface [13]. Wear is directly related to the fatigue life of part [14]. Bayer [15] found that wear was the major factor to define and limit the life of a machine element such as dies or molds; it costed 1% to 4% of the gross national product (GNP) of the nation. The majority of machine elements including bearings, gears, artificial hips or knee joints, or seals were failed mainly due to wear and friction [3].
The studies have shown that friction is related to numerous factors such as relative displacement, velocity, and acceleration of motion, contact force, surface finish, and the properties of materials of parts [16]. In particular, friction varied during the contact when the lubrication condition changed. Tribology is to investigate the relation of friction and lubrication [17].
The experimental approach has been the most reliable way to investigate the tribological behaviors of machine elements; however, it is not always an effective way to take into consideration of a large number of design factors. Various theoretical approaches were developed to evaluate different design alternatives [18]. However, predicting the friction using a theoretical or numerical model is not a trivial task. Due to the complexity and uncertainties, developing a friction-prediction model or tool turned to be an elusive goal for researchers in 50 years [13,19].
Bi and Mueller [20] have had a thorough discussion on available simulation tools to model tribological behaviors at pin-plate interfaces. This article makes an extension by modeling the friction of sealing parts and investigating software tools which are capable of predicting the friction subjected to given operating and lubrication conditions. The rest of the article is organized as follows: (1) The friction characteristics are discussed to identify the major factors that affect tribological behaviors in Section 2; (2) existing literature on mixed lubrication frictions are surveyed and summarized in Section 3; (3) different techniques for numerical simulation, as well as the applications in solving different engineering problems, are described and compared in Section 4; (4) the practical challenges in developing a reliable modelling tool for friction prediction are discussed in Section 5.

Characteristics and factors of friction
Friction is the resistance to the relative motion at the interface of a solid to other solid or fluid. Friction was classified in terms of lubricant conditions [21]. Friction causes wear of solids in different forms such as fatigue, abrasion, corrugation, ploughing, cavitation, and erosion. Friction varies with the lubrication condition, which can be characterized by a Stribeck curve. As shown in Fig. 1, a Stribeck curve has three regions, i.e., boundary regime, mixed regime, and hydrodynamic regime [23]. The amount of friction is measured by the coefficient of friction (CoF), and three main factors to determine CoF is pressure, velocity, and viscosity. For normalization, a traditional Stribeck curve used a dimensionless Hersey number, which was calculated from (viscosity ) × (speed N) / (Load P) [24,25]. If design variables other than viscosity, speed, or load have to be considered, Stribeck curves can be extended in three or even high dimensions. For example, Wang [26] and Wang et al. [27] defined 3-D Stribeck surfaces to evaluate the impact of loads and velocities in journal bearings, respectively. In a three-dimensional Stribeck surface, x-and y-axes were for the changes of the load and velocity, and the z-coordinate was for CoF subjected to specified load and velocity.
In the full lubrication condition (phase III) and when the shear stress reached a certain level, lubricant tended to have a constant viscosity instead of reducing with an increase of shear stress. As shown in Fig. 2, the flow curve of lubricant shows a second Newtonian region with the second Newtonian viscosity. Katyal and Kumar [28]  modelled the shear-thinning phenomenon of EHL point contact by double Newtonian Carreau Viscosity equation, and the result showed the effect of the second Newtonian viscosity on the film thickness was determined by load, scale, and the relative speed at contacts.
The Newtonian model applied to the simplest rheological lubrication model where the film thickness between two surfaces never is broken. When oil film was broken, the asperities made direct contacts, and the behavior of the oil became non-Newtonian. For the contact surfaces with given roughness, the transition from the Newtonian to the non-Newtonian model might just be justified by the following ratio defined as [29] total initial surface roughness elastohydrodynamic oil film thickness D  For a value of D under 0.1, the surfaces were totally free of asperity contacts, and the Newtonian model was appropriate. The inverse of D was referred as lambda ratio () and used widely as a criterion to justify if the full-lubrication was formed.
The condition of  > 2 tended to a fully lubrication, and 1.5 <  < 2 corresponded to the mixed boundary lubrication [4].
Seabra et al. [30] studied the friction at roller-inner contacts subjected to the EHL condition. The behavior of the lubricant was represented by two rheological laws, i.e., the viscous and viscoelastic models of a Ree-Eyring fluid, and the Barus and Roeland expressions were used to describe piezov-iscosity and thermoviscosity laws of lubricants. Sander et al. [31] integrated the experimental and simulation methods to investigate the mixed lubrication in journal bearings, and they found that both of piezoviscous effect and the non-Newtonian behavior must be considered to determine the friction power losses under different dynamic loads and motion speeds.
Fang et al. [32] considered single-point contact under the EHL condition, the theory of thermal activation energy was used to describe the viscous behaviors, and the constitutive relations was expressed by the Eyrin-based rheological model as where  is the shear stress,  is the stress-strain rate, 0  is the Eyring representative stress before that the non-Newtonian behavior dominates the lubricant;  is the pressure-temperature dependency of the viscosity of: where S is called as the viscosity modification factor, and A and B are pressure and temperature indices, respectively. Bou-Chakra et al. [33] argued that the non-Newtonian effect of EHL was viscoelasticity, shearthinning, and plasticity, and the Ree-Eyring model was used to describe the rheology of a layer as where  is the shear rate,  is shear stress, t is the variable along with the thickness, and G is the shear modulus. h is central film thickness, ho is surface layer thickness, P is mean contact pressure,  is pressure coefficient of viscosity,  L is the viscosity of the adsorbed surface layer, 0L  and 0H  are the Newtonian limits of the lubricant and the surface layer, respectively. Anuradha and Kumar [34] investigated the rheological behaviors with the consideration of the shear-thinning effect, and the Carreau viscosity 210 Friction 9(2): 207-227 (2021) | https://mc03.manuscriptcentral.com/friction model was applied to represent the shearingthinning effect. The dimensionless rheological model for spur gears was: where     0 / is the low shear viscosity, U is the dimensionless speed parameter,      / u Y is the shear rate, n is the power-law index, and cr G is the critical stress corresponding to the Newtonian limit of lubrication.
For any of the three phases in Fig. 1, many factors contributed to CoF occulting to the contact surfaces. While it was infeasible to investigate all affecting factors, researchers have identified the following five variables as the main factors of CoF.
1) Normal Load. In a simplified friction model, CoF is constant, and the friction force is proportional to a normal load. While CoF depends on the contact condition in reality, and it varies with the change of asperity contacts, the deformation and stress distribution. An increase of normal load leads to a reduction of the gap of two rough surfaces, which involves more asperities in physical contacts and more plastic deformations over pre-contacting asperities. The shearing factor of the boundary layers at two contact surfaces causes the friction. A continuous decrease of the gap will increase the amount of friction force until all of the load is carried solely by asperities, and the friction force will not be increased further [35].
2) Viscosity. Viscosity is considered in determining the Hersey number in a Stribeck curve. When the area and gap of two contact surfaces are given, the resistance force depends on the viscosity and the shear rate. The higher the fluid viscosity is, the higher the resistance force is. The impact of the shear rate relates to the lubricant types, i.e., Newtonian or non-Newtonian. For a Newtonian fluid, the friction force is proportional to the velocity of the relative motion. For non-Newtonian fluid, the friction force is nonlinear, which depends not only on the velocity of the relative motion, but also the change of viscosity with respect to the shear rate.
3) Shapes. The conformance at the macro-scale of the surfaces affects the distribution of contact stress. The contacts can be made on points, lines, or areas. The visible features not only affect the distribution for contact stress, but also affect the interaction with the fluid flow. It brings an issue to define an EHL model since the boundary conditions have to be specified clearly; while an inlet or outlet of the fluid flow is open around the contacting points, lines or areas. Textures such as dimples and bumps over a surface can be treated as the features at a macro-scale in numerical simulation. A number of researches were reported on the impact of surface textures on tribological behaviors. For example, the Stribeck curves were investigated for journal bearings with dimpled bushings [36,37]. In their experiments, dimples were fabricated by chemical etching; the friction forces of the bearings were measured under different loads, lubrication types, the shapes and depths of dimples. Galda et al. [38] studied the impact of geometrical texture on the friction force in sliding, and their experiments indicated that the lubrication conditions for sliding were improved by appropriate shape, and dimensions of oil pockets. Menezes et al. [39] measured the friction forces for the contacts with four textures and surface roughness. It was found that CoF was insensitive to roughness when the texture remained the same; while CoF was significantly affected by the changes of surface textures even the surface roughness was the same. 4) Surface Finish. The surface finish affects the friction greatly under a boundary or mixed lubrication condition. When the surface roughness is significant in comparison with the thickness of the fluid film, two surfaces will have the contacts at asperities to share the normal load. The lubrication condition relates to a number of factors such as surface topography, the pressure of asperity contacts, fluid pressure, distribution of film thickness, and even the temperature at contacts [40]. It was found that the friction was affected by the surface texture in self-lubricated bearings [41]. The impact of surface texture and roughness on an EHL model was studied [42]; it turned into a semi-analytical model with Reynold's and elasticity equations. The surface pattern also affects the characteristics of the lubrication at contacts; Akbarzadeh and Khonsari [23,43] investigated the impact of surface patterns in the contacts of alumina alloy and steel plate in a dry condition and found that the friction depended on surface patterns rather than surface roughness [44]. 5) Temperature. The wasted energy by friction is transferred into the heat form. It will change the temperatures of objects. While the temperature affects material properties such as Young's modulus of solids and viscosity of fluids. If the temperature involves a significant change, it impacts on material properties and stress distribution. For example, shear-thinning can be introduced if the reduction of viscosity is due to an increase of fluid temperature. The characteristic of shear-thinning in gear transmissions was investigated [43]. This will reduce the thickness of the fluid film and thus the percentage of load share by fluid pressure. The dependence of temperature and pressure on density, pressure, viscosity, and localized shear stress were explored [45]. The thermal effect could be significant to change the friction-induced motion behavior in stick-slip actuators [46]. The mechanism is that the friction induces heat, which causes temperature rise and deforms the contact surface, and eventually increases the friction.

Friction in mixed lubrication conditions
Many researchers have contributed to model friction and wear in different lubrication conditions. Ludema [47], Siniawski et al. [48], and Yang et al. [49] found that over 300 equations were developed as friction and wear models more than 1,000 other papers discussed extensively on friction and wear. In an EHL model with mixed lubrication conditions, the friction consists of (1) the resistance at adhesive junctions and (2) the resistance caused by the viscosity of the fluid flow. Kragelsky [50] was the first author to discuss the relation of the friction and the elastic deformation and the shear strength at adhesion, where the asperities at the contacts were treated as a set of the rods with different lengths and fixed at their ends. The proposed statistical model showed that the friction force was determined by materials, pressure, and height distributions of asperities. Paouris et al. [51] proposed an empirical tribo-dynamic model to predict the friction under a mixed boundary condition in a gear transmission; the friction consisted of the shearing of tribo-film as asperity summits and the dry sliding frictions over the contacted asperities. Figure 3 gave the model of surface contacts at the mixed lubrication conditions [24]. When a normal load was applied, two surfaces made the contacts at some asperities, but the normal load was carried not only by the contact asperities but also the pressure of the lubricant layer. The load-sharing concept was introduced by Johnson et al. [52] where the normal load was shared by the deformed asperities and the fluid pressure at the contact regime. However, quantifying the load from asperities was not straightforward. Here, relevant literature were summarized based on different applications.

Resistance at asperity contacts
When two rough surfaces are forced to make contacts against each other, the contacts occur to peak asperities at the first. If the external load keeps increasing, the number of contacted asperities increases, so are the deformations of asperities. However, the percentage of the actual area of the contacted asperities over the total nominal area is extremely small, which typically varies from 0.00001% to 0.005% [21]. In addition, either plastic  | https://mc03.manuscriptcentral.com/friction or elastic deformations occur to the asperities; this leads to the resistance when two surfaces are involved in a relative motion. The resistance could be calculated by the magnitudes of the plastic strains over the asperities [53].
The statistical approaches were used to correspond to an average gap of two surfaces to the normal load applied at the interface. Greenwood and Williamson [54] proposed a contact model (i.e., the G−W model) statistically to evaluate the effective areas of plastic and elastic contacts, and the model was further used to estimate the contact stress in an ensemble of asperities. The G−W model was then extended [55] for multi-asperity contacts of two rough surfaces. The G−W model and its variants were used widely to analyze the contacts of two rough surfaces [56]. To make the G−W model more practical, Gupta et al. [57] and Gupta [58] evaluated the friction of the bearings analytically with the consideration of geometry, operating, and lubrication conditions. However, the G−W model and its variations were shown some limitations: (1) While the surface geometry is continual, but the G−W model treated it as discrete asperities; the materials beyond asperities were discarded even though it affected the deformation of asperities. (2) The G−W model was not able to consider the impacts of neighboring asperities. (3) It was critical to determine the number of asperities per unit; while it was difficult to characterize asperities and distribution adequately [59]. In the statistic model by Ogilvy [60], the surface roughness was assumed to be randomly Gaussian distribution; the Hertzian equation was applied to determine the number, locations, and sizes of contact asperities. The resistance was estimated as the multiplication of the shearing strength and actual contact area by contacted asperities. Molinari et al. [61] considered the resistance under dry sliding condition; they used the Lagrangian method to model large plastic deformations and nonlinear material properties, the purpose of using the Lagrangian method was to eliminate the mesh distortion by deformation continuously.
Deformation of an object is always associated with energy; alternatively, the law of energy conservation can be used to develop a contact model of surfaces. Tian and Bhshan [62] developed an energy-based model for elastic and elastic-plastic surface contacts.
With the specified load, the actual contact area was estimated for the minimized total complementary potential energy. In the simulation model by Hu et al. [63], the minimization of complementary potential energy was integrated with the Fast Fourier transform (FFT) to calculate the contact area.
Siniawski et al. [48] corresponded to the wear with the number of loading cycles directly based on a kinetic model of abrasion for sliding bearings.

Resistance of elastohydrodynamic lubrication (EHL)
The fluid flow has the second type of resistance to the relative motion of two surfaces. An EHL model is introduced to quantify the shearing resistance at the interface. The resistance from the fluid flow depends on pressure distribution, which is affected by the normal load, surface geometry, operating condition, and fluid properties. The Reynolds equation is widely used to model an EHL problem [64]. Figure 4 shows a thin film over a finite area of two surfaces in a relative motion. One surface is assumed flat in the X−Z plane, and the other surface is curved. The velocities of the fluid flow along X, Y, and Z are denoted as U, V, and W, respectively. At the given moment, h(x, z) is defined as the distance of the corresponding points at two surfaces with the same X and Z coordinates of (x, z). Accordingly, the following boundary conditions can be defined for this finite area: , , Since it is a thin film, it is assumed that p/y = 0, and the pressure changes   / p x and p/z along X and Z are independent of y; the Reynold's equation for the mass conservation of incompressible fluid is expressed as  While surface texture changes the film thickness then affects pressure distribution, early Reynold's equations were able to take into consideration of roughness stochastically in one-direction, i.e., the ridges were transversely or longitudinally oriented. Patir and Cheng [65,66] were able to extend the stochastic Reynold's equation for a generic roughness model. Due to randomness of height, the film thickness is also randomly distributed. It seems more practical to use mean pressure or film thickness. Therefore, the Reynold's equation was modified by averaging fluid flows, and roughness was modeled by the factors of pressure flow, which were determined from the numerical simulation. In such a way, and the EHL modelling problem was reduced to estimate the flow factors between two nominal boundaries. Tripp [67] adopted the perturbation expansion of pressure to estimate the flow factors in the Reynold's equation. The pressure in the fluid in the mixed lubricated area shares the normal load with contacted asperities. Figure 5 showed that the resistant force, i.e., friction, was the sum of the resistances at asperity contacts and the EHL layer [68]. Payvar and Salant [69] proposed the cavitation index to measure the shared load in mechanical seals. Harp and Salant [70] modeled the lubrication conditions in hydrostatic seals; the considered factors included fluid properties, solid contacts, and deformation of the seal. Both inter-asperity and macroscopic cavitation were considered to improve the Reynold's equation [71]. The effects of material properties were also investigated by some researchers. For example, Stoyanov et al. [72] investigated the friction force of WC materials over W in both dry and lubricated contacts. They found that the dry sliding of WC over W involved in the highest friction force; in addition, sliding in a dry condition increased the surface roughness, thus increased the friction for time being.

Fluids types and other factors
Resistance by shearing in the film depends on the viscosity of the fluid. The Newtonian fluid has a constant value of viscosity; while the viscosity of non-Newtonian fluid varies with the changes of temperature or shear rate. The effect on viscosity by shear rate is called as shear-thinning. It was found that using Newtonian fluid in the Reynold's equation was inappropriate to analyze tribological behaviors of fluid in mixed lubrication conditions [73−75] indicated that the resistance was determined by both film thickness and shear-thinning, and they further used the power-law rule in an empirical model to measure the impact of shear-thinning.
The non-Newtonian fluid has been widely investigated in the past 30 years. For example, Sharif et al. [76] proposed a thermal elastohydrodynamic lubrication model for non-Newtonian fluid to calculate film thickness and traction force for a type of rolling rigs. An EHL model used to treat the contact regime as isothermal since the heat from friction was neglected; an EHL model with consideration of thermal transfer was called as a thermal EHL (TEHL) model. Hacioğlu and Dursunkaya [77] studied the dependence of lubricant and relative motion for slider-crank pistons whose surface texture was modeled. When two surfaces have a fluid film in between, the shear stress in fluid resists the relative motion. Otero et al. [78] proposed a model to predict Mixed lubrication is certainly relevant to applications which involve numerous operating factors. Taking the bearings as an example, some major factors included the eccentricity ratio, surface finish, lubrication properties, pressure, and contacts of asperities [80]. Roshan et al. [81] studied the impact of lubricant ingredients and properties on the formation of boundary lubrication conditions. Fryza et al. [82] investigated the squeeze film action when EHL was subjected to an impact loading or abrupt change of motion experimentally; they found that entrapped film shape was mainly determined by the loading speed and the thickness of the central film by the approaching speed and viscosity.

Tribological behaviors
In the mixed lubrication, the friction is complex due to the coupling effect of contact mechanics and fluid dynamics. The deformation of asperities and the motion of fluid flow must be considered together to find the distributions of film thickness and pressure at the contacts. As shown in Table 1, Chang [83] argued that four governing equations may be assembled to model the tribo-contacts in the micro-EHL regime. For example, the Reynold's equation was developed as a simplified Navier− Stokes equation for the mass conservation in a narrow flow channel [35]. Shirzadegan et al. [84] developed the customized Reynold equation to simulate the EHL line contacts of a cam-roller follower; the corresponding rheological model considered rolling crowning, edge geometries, and the variants of pressure and film thicknesses. Table 1 Governing equations for mixed lubrication.

Equation Description
Reynold's equation Describes the relation of the hydrodynamic pressure with other variables in lubricated film.

Elasticity equation
Represents the deformation of the contact surface under fluid pressure and normal load.

Load equation
Represents the equilibrium condition of external load and the force of fluid pressure.

Energy equation
Governs energy conservation with consideration of heat generation and transfer at the interface.
It will be ideal that the equations in Table 1 are solved simultaneously to obtain a valid solution, which satisfies all of the constraints in four disciplines. It is unfortunate that up to date, no solver is able to deal with these equations in one simulation model; researchers made their efforts in solving these coupled equations iteratively in a sequence. For example, Ai et al. [1] decoupled the Reynold's equation and the model of contact asperities for journal bearings. Ghahbavieh et al. [85] established respective EHL model and solid contact model to predict film thickness and friction in gearing, the design variables in these two models included loads, roughness, hardness, and rolling speeds of gears. Olver and Spikers [86] modeled the heat transfer by considering the traction [56] to estimate friction in different traction regimes.
The mechanics at molecule-scale may give a new perspective to understand friction and wear; since the film thickness is at the same magnitude scale with molecules size. Fillot et al. [87] introduced molecular dynamics in continuum mechanics to quantify the slips at the nano-scale; a semi-analytical model was developed to model the slip over surfaces with the consideration of pressure, film thickness, and relative speeds. Ghaednia and Jackson [88] investigated the major factors that affecting nanoparticle friction; two models were developed to deal with roughness features and nano-particles respectively. They found that nano-particles were able to reduce friction due to the reduction of actual contact area; therefore, nano-particle sizes and distributions affected tribological behaviors when nano-lubricants were applied.

Numerical simulation
Due to the lacking of effective simulation techniques, experiments are still widely used in the study of tribology. For example, Bair [45] and Bair et al. [75] experimented to develop a shear rheological model for the lubrication with concentrated contacts. Lee and Polycarpou [89] used experiments to estimate the coefficient of friction (CoF) in different roughness, dwell time, displacement rate, lubrication trace, and wear debris at reciprocating interfaces. Note that the experiments should be the last resolution due to their disadvantages in time, cost, and the limitation of tested options.
Semi-experimental approaches may alleviate some issues of experimental approaches. A semi-experimental approach uses an empirical model which can be calibrated and refined by the data from the experiments. Zhao et al. [90] and Zheng and Liu [91] used the neuro-based methods to investigate the wear in a number of products. Kumar et al. [92] used an artificial neural network (ANN) to predict the friction and wear rate, and the experimental data was collected to train and test the ANN model. Nagaraj et al. [93] and Senatore et al. [94] used similar approaches to develop friction-prediction models; experiments were conducted to acquire data and train the ANN models. Okuyucu et al. [95] used an ANN model to investigate the friction stir welding of aluminum plates. Varade and Kharde [96] introduced the Taguchi method in their ANN model to find the wear rate of the PTFE composite. Similarly, Prater et al. [97] used the Taguchi method for statistical analysis to calculate the friction force in friction stir welding.
Note that either an experimental or semi-experimental method tends to be expensive and very time consuming [98]. Numerical simulation is desirable in predicting friction and wear subjected to different operating conditions. Numerous studies were reported towards this goal. Generally, the developed simulation models fall into two groups, i.e., statistic models or deterministic models.

Statistic models
Most of the friction models were statistical [99]. This is because the friction occurs to the micro-level surface structure while the surface finish has a random distribution in heights. A statistical model is suitable to represent the contacts of asperities and the tribological behavior at contacts. The majority of early friction models, such as the G−W model and its extensions, were statistic models. For example, the model by Francis [100] used the concept of probability in estimating the average deformation of the interface at two rough surfaces; they assumed that the actual contact area of two nominal planes was to that with the same load over a smooth plane. The equivalent area was estimated based on the Gaussian distribution.

Deterministic models
Statistic models have their limitations for microlevel statistical behaviors, which are inappropriate to assess the impact of features at the macro-level such as bumps or dimples at contacts. With a deterministic model, visible surface features are graphically modelled explicitly, and both solid contacts and the fluid flow at the contacts are simulated using finite element analysis (FEA) packages; for example, surface finishes were modeled as a set of individual asperities in the deterministic model [13]. By all means, computer programs are needed to solve deterministic friction models.

Specialized programs
Technically, the set of the governing equations in Table 1 can be numerically solved by specialized programs, and many researchers contributed to developing specialized programs in solving EHL problems. For example, Ai et al. [1] adopted the Garlerkin method to represent the governing equations and applied the Newtown−Raphson method to solve the assembled EHL model at the system-level. Andersson [2] formulated an ordinary dynamic model to represent the wearing at rolling and sliding contacts, and they suggested the technique to simplify differential terms to accelerate the simulation. To deal with the coupling of viscosity, pressure, and temperature, Bair [45] used the free volume theory to establish a simulation model where the film thickness was determined with the consideration of heat transfer. Beheshti and Khonsari [3] investigated the wear in the mixed EHL regime by estimating the stress at contacts by elasto-plastic asperities and analyzing the temperature with consideration of heat generation and transfer. Han et al. [101] developed a specialized program to predict the friction in helical gears. Bjorling et al. [19] attempted to develop a practical friction-prediction tool to avoid tuning of lubrication properties.
Friction and wear vary along with time. Kim et al. [102] developed a simulation model to calculate | https://mc03.manuscriptcentral.com/friction the wear rate of metal parts under oscillation. The model involved in adjustable meshes to compensate for the material loss by wears. The computation was reduced by discretizing wear propagation and adopting the extrapolation scheme.
Molecular dynamics has been widely used to model the interaction of solids and liquid at the molecular level. Mendonca et al. [103] developed a molecule-scale model to study ionic lubricant over metallic surfaces, and they found that the friction was mainly affected by flow velocity, load, surface texture, as well as modular structure such as length of molecular chains. The molecule-scale model by Stoyanov et al. [72] was used to study the friction between WC and W surfaces; they found that that the grains near the contact surfaces were refined due to friction; such refinements, in turn, affected the friction and caused the rotation of grains and the fatigue by multi-axial shears.

Generic FEA packages
Commercial FEA packages are generic and versatile, and they are widely used in multidisciplinary studies, such as developing simulation-based tools to predict friction [2]. Carden [25] discussed the applications of different FEA packages, such as FAST, PISDYN, ENGDYN, WAVE, and VALDYN, in predicting the friction forces of machine elements. The friction assessment simulation tool (FAST) was used to study the friction in engines, and the program was able to consider pressure, the masses of piston and rod, and the ring tension in the model to calculate the friction force on the piston.
Steen [104] surveyed the friction models developed in the Abaqus program. For example, Abaqus was used by Abu-Bakar and Ouyang [14] to determine wear rate and pressure under different surface finishes; experiments were conducted to measure surface profiles overtime at the validation of the simulation results. Ray [105] used Abaqus as the simulation tool to study the friction in self-lubricating bearings. McColl et al. [106] and Liu et al. [107] corresponded to the normal stress with the friction over contact surfaces which were subjected to a fretting load, a customized sub-routine was integrated with the Abaqus program to modify the meshes for loss of materials by wear. ANSYS Multiphysics was used by Belhocine et al. [108] to predict the friction force of braking discs; the impact of temperature changes was emphasized. Podra and Andersson [18] used ANSYS to calculate the sliding friction for the spherical pin-on-disc pair in dry lubrication; the simulation and experimental results gave a discrepancy of 40%−60%. LS−DYNA was used by Moghaddam [109] to model the surface roughness and predict the adhesion of surface contacts. Gustafsson [110] integrated multiple tools at different stages of simulation to predict the wear: ANSA was used for preprocessing, LS−DYNA was used for processing, and META was used for post-processing.
MSC. Marc MENTAT package was used by Békési [98] to investigate the friction of elastomers; they developed two models for sliding of a round asperity and a pin over rubber plate, respectively. Jiang et al. [111] modeled d the asperity contacts subjected to rolling, sliding, and the combined rolling and sliding. They found the massive computation was needed: even giving the asperity size of  0.02 m and the plate size of 0.2 m by 0.2 m; it took hours to finish the simulation. AutoForm was used by Wang et al. [112] to optimize sheet metal stamping with the friction reduction, and considered variables were lubrication conditions, pressure, and surface coating. Škurić et al. [113] proposed a lubricated contact model for metal forming processes, and the OpenFOAM was used as the solver to evaluate the hyperelastoplastic deformation of solid contacts. Both the deformation of asperities and the changes of rheology properties were taken into consideration.

Existing friction models
Understanding the tribological behavior of some machine elements such as seals, bearing, centrifugal pumps, mixers, turbines, engines, compressors, and propeller shafts is critical to improve product performance [114]. However, the level of difficulty to model and solve a friction model varies greatly with products, geometries, and operating conditions. In this section, existing friction models are classified by contact types.

Asperity contacts
Some simulation models were developed to study the friction at the contacts of asperities. García and Martini [115] represented asperities as partial spheres to predict the static friction on the plane, the surface roughness was modeled by small spheres at varying heights from the nominal plane; and the discrete convolution Fast Fourier Transform (DC−FFT) was applied to consider the effect of cold hardening. Békési [98] proposed the simulation models to estimate friction forces of rubber materials for (a) a round asperity sliding over a line, and (b) a pin with spherical end sliding over a plane. Bjorling et al. [19] simulated the contacts of asperities to a 3-D Stribeck surface by considering the variants of loads and sliding speeds. Ford [55] simplified the asperities as cylindrical pins to model the friction at multi-asperity contacts. Gustafsson [110] conducted the static analysis for the stress distributions of the asperity in simplified 2-D and 3-D contact models. Feng and Guo [116] represented a foil surface as a set of bumps, each bump was corresponded to a mass-spring system to calculate the friction in optimizing the structure of a foil bearing.

Line contacts
Journal bearings and gears are among the widely used machine elements; therefore, the lubrication mechanisms at line contacts have attracted a great deal of attention.
Gelinck and Schipper [117] developed a mixed lubrication model for line contacts to define Stribeck curves; it was used to predict the transition stage from one lubrication regime to another. Ai et al. [1] simplified the contact in a journal bearing as the line contact, and developed an FE model to evaluate the impact of film pressure and thickness on friction. A lubrication regime in gearing is either mixed or partial EHL. Andersson [2] considered both rolling and sliding in their friction model of gear transmission. Stahl and Jacobson [118] investigated the impact of the Hertzian pressure on film thickness at line contacts using a rheological model. Ma et al. [119] predicted the friction between the cylindrical liner and the piston ring with consideration of surface roughness, temperature, pressure, and lubrication. Akbarzaheh and Khonsari [43] modelled the meshing of spur gears to investigate the shear-thinning. Beheshti and Khonsari [3] modeled the wear at line contacts for both of the dry and lubrication conditions; the empirical formula was used to estimate contact pressure and a statistic model was used to calculate film thickness. Han et al. [101] studied the dependence of film parameters on surface roughness and load variation for helical gears. Gay [105] developed a friction model for self-lubricating composite bearings, and two main factors were the constitutive materials and contact pressure. Ghahnavieh et al. [85] simulated the meshing of straight bevel gears under the mixed lubrication to predict the film distribution and friction force.

Area contacts
The G−W model [54] was classic for the friction prediction of the surface contacts. In the G−W model, the friction was calculated statistically based on the estimation of deformed plastic and elastic areas; it was assumed that the gap of two nominal surfaces was given.
Technically, the geometries or textures of contact surfaces can be modeled in detail. For example, Harp and Salant [71] modeled the fluid flow between two rough surfaces in defining a modified Reynold's equation which included film thickness, pressure, and surface roughness for the flow factors. Berger [16] proposed different friction models for various surface contacts in dynamic simulation. Kim et al. [120] developed a friction model for the surface contacts in the small intestine in optimizing a self-propelled capsule endoscope to reduce friction. Abu-Bakar and Ouyang [14] developed a friction model for a disc brake with the area contact; they modelled the real surface topography in the simulation. Hacioğlu and Dursunkaya [77] modeled the surface roughness to consider its impact on fluid flow in developing the Reynold's equation for the pistoncylinder bearings. Ghaednia and Jackson [85] attempted to model surface contacts by nanoparticles to include both nano-level and micro-level features. Belhocine et al. [108] looked into the temperature effects in a disc brake with surface contacts.
The friction force is certainly relevant to the material properties of parts. For example, an elastomer part moving over a rough surface may have the friction | https://mc03.manuscriptcentral.com/friction by adhesion, abrasion, and hysteresis of materials [98]; the friction of elastomers was the synthetic effect of adhesion and internal friction from the hysteresis of deformation of counterpart. Pálfi et al. [121] discussed the friction of the rubber plate over a rough surface caused by hysteresis. Scaraggi and Persson [122] developed an analytical model to predict the friction of a rubber part in rolling.

Challenges in numerical simulation
Many researchers have contributed to develop friction and wear models for different parts under mixed or full lubrication conditions. Unfortunately, the authors have not found any software tool, which has been proven its practical value in supporting the product designs for better lubrication performances. While developing such a simulation tool sounds feasible technically, this section discusses some challenges to be addressed in the implementation.

Data preparation for simulation
A simulation model always needs some inputs such as part models, material properties, boundary, and loading conditions. Sufficient and valid inputs are critical to generate acceptable and reliable outputs from the simulation. In preparing for input data for an EHL model, Bayer [15] suggested to determine (1) major design factors such as part models (geometry, dimension, surface texture and roughness, and material) and operating factors such as lubrication types, ambient temperature, and motion conditions, (2) contact conditions at the interfaces including rolling or slipping types, motion types (unidirectional/reciprocating, low/ medium/ high speed, and fretting/gross sliding), lubricating (dry/wet/ mixed), loading (light/medium/ heavy, constant/ varying, steady/gradual/impact; contact nature (point/line/area, conforming/non-conforming); level of contact stress (elastic/plastic/mixed); temperature (constant/varying, low/medium/high), and environments (hostile/non-hostile, with or without abrasive particles).
By all means, it is challenging to quantify the aforementioned factors: for example, characterizing surface roughness is difficult seems the height is randomly distributed on the contact surface. Summing all resistance types as the friction force also involves in a high-level of uncertainty since different resistance types are affected by properties, coating, temperature, and in particular, micro-level randomness and impurities of materials. In addition, the basic tribological characteristic of materials, such as coefficient of friction (CoF) in the dry testing condition, is assumed in the simulation. There is the need for a minimal set of experiments to acquire the input data for simulation models [81].

Surface modeling
To model and analyze EHL behaviors, the features of the contact surfaces at both macro-and micro-levels have to be modelled. However, when the features at the micro-level are modeled deterministically, the mesh density must be fine and small enough in comparison to micro-level features, this brings the challenging to run the simulation model since the computation can be exponentially increased with the number of nodes and elements.
Since modeling surface roughness is critical to understand its impact on the formation of mixed lubrication conditions, how to represent the surface roughness in the simulation was investigated by a number of researchers. Akbarzadeh and Khonsari [23] proposed to use the ratio  of the film thickness and the standard deviation of asperity heights as the affecting measure of lubrication conditions: the mixed lubrication occurred when  was in the range of (1.0, 3.0) and the hydrodynamic lubrication occurred when  is larger than 3.0. The characterization of surface roughness was extensively discussed [123−125]. For example, Bhushan [124] classified the characteristic of contact surfaces: a contact surface was either of a non-homogeneous or homogeneous surface. The homogeneity of surfaces could be deterministic or random, and a random surface could be anisotropic non-Gaussian distribution or isotropic Gaussian distribution. Figures 6 and 7 gave two schematics to characterize surface roughness [3,124].
Early efforts in modelling contact surfaces used the parametric method, in which the profile |www.Springer.com/journal/40544 | Friction http://friction.tsinghuajournals.com  of a surface was modelled as the pattern of parametrized features such as cylinders, hemispheres, cones with spherical caps, ellipsoids, or sinusoidal surfaces [126]. Patir [127] used the linear transformation method over the random matrix based on statistical properties to model rough surfaces. The mating relations of the contact surfaces with such features must be defined with cautions to avoid interference [128]. Syafa'at et al. [129] treated a surface finish as a set of micro-scale wavy features with a height amplitude of 2%−5% of the wavelength; the surface roughness was modelled as spherical asperities in predicting sliding wears in Abaqus. Thompson [130] proposed four surface finish models and implemented them using ANSYS Parametric Design Language (APDL); the comparison was made to identify the advantages and disadvantages of each method.
Tavares [131] and McCormick and Duho [132] utilized the measurement in modelling surface roughness in simulation. Surfaces can be measured by one of optical-based, stylus-based, or microscopicbased techniques [40]. Lee and Cheng [133] tried to model contact surfaces based on real profiles. The contact surfaces were considered simultaneously to define the contacts with equivalent roughness; the model was used to evaluate the deformations under varying loads. Hu and Tonder [134] introduced a 2-D filter to model 3-D surfaces, the spectrum analysis was conducted to determine the coefficients of the filter. Talyor et al. [135] suggested characterizing surface by the roughness heights and asperity diameters. In characterizing asperities, Wu [99] proposed an auto-correction method to refine the height distribution of asperities; the curvatures of asperities were estimated as the root mean square curvature of profile. Pálfi et al. [121] defined rough surfaces used the wavelengths and height amplitudes to characterize rough surfaces, and further generated non-uniform rational basis splines (NURBS) to model surfaces. Figure 8 shows the example of surface characteristics with different sizes, i.e., micro-level roughness model by Bergström [136] and Pawlus et al. [137] shown in Fig. 8(a), and a bump feature at macro-level as shown in Fig. 8(b). To represent the features at two levels in one surface model, surface roughness can be mapped on the surface with the macro-level features, which is illustrated in Fig. 9. It will never be overemphasized that modelling of actual surface finish demands a huge amount of computation, memory, and storage.

Idealization
To obtain a valid simulation model, some trials and errors are needed to deal with uncertainties and incomplete data. To develop a friction-prediction model for EHL problems, the simulation model should be idealized on the following factors (1) exploring the feasibility of using deterministic surface models to represent statistic properties of surface roughness at contacts, (2) analyzing design factors and identifying the limited number of main factors, (3) simplifying boundary conditions with minimized uncertainties, and (4) defining the material properties of parts, coating, and lubrications reliably with the minimized assumptions.

Multidisciplinary coupling
The tribological behaviors at the contacts with lubrication are extremely complex due to the multi-   Example of the realistic surface model disciplinary coupling of solid mechanics, contact mechanics, fluid dynamics, and system dynamics modeling [16]. Each discipline brings a number of governing equations, which causes the difficulty to solve these equations simultaneously. To the authors' knowledge, no commercial simulation tool is available to integrate the models in contact mechanics and fluid dynamics as one system model in the simulation. Alternatively, some commercial FEA packages are able to decompose an EHL problem into two sub-models in contact mechanics and fluid dynamics, respectively, and further solve the sub-models as a system-level solution sequentially and iteratively. Due to the strong multidisciplinary couplings, there is no guarantee that the iterating process will be converged; even though it will, the solving process might take an extremely long time to obtain the results in an acceptable time duration.

Challenge of computation
In an EHL problem, contact surfaces interact at the asperity-scale, and this is accompanied by solid deformation, formation and fraction of fluid films, and even tribo-chemical reactions. Therefore, the meshes in solids and fluid volume with a very high are needed to represent the details and generate reasonable results [129, 138−140]. Chang [140] discussed the dependency of the solving steps and time for a converged solution on the numbers of nodes and elements; they found that several thousand-time steps were taken to obtain the solutions for the EHL model with the surface finish only; while over fifteen thousand time steps were taken when asperity interactions were modelled [83,139]. It is well-known that modelling of surface finish involves in the highly demand of computation. Karpenko and Akay [141] proposed the influence coefficients to represent surface finish in estimating the sliding friction from shearing. They concluded that the friction was linearly reduced with an increase of surface roughness.

Verification and validation
The simulation results must be validated against some trusted results that are either obtained from experiments or the proven data by others. In the verification and validation, one has to ensure that the results from experiments and simulations are comparable [142]: Both simulation and experiment results must be reproducible; the varying ranges of multiple measurements or simulations under the same conditions must be at an acceptable level of repeatability and accuracy. article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not in-cluded in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder.
Zhuming BI. He received the bachelor degree in manufacturing engineering in 1987 from Harbin University of Science and Technology, the master of Science in mechanical engineering and the Ph.D. degree in mechatronic control and automation from Harbin Institute of Technology in 1991 and 1994, respectively. He received the second Ph.D. degree in mechanical engineering from University of Saskatchewan in 2002. He is a professor of mechanical engineering in Purdue University Fort Wayne, and his research interests are modelling and simulation, manufacturing systems, robotics, and automation. Donald W. MUELLER. He received all of his bachelor degree, master of science, and Ph.D. degree in mechanical engineering from University of Missouri-Rolla. He is an associate professor of mechanical engineering in Purdue University Fort Wayne, and his research interests are thermal sciences, machine design, and numerical methods.