On an Elastoplastic Sliding Model for a Coated Single

In this study, a sliding friction model for coated single asperity contacts is proposed. A displacement-driven layered contact algorithm is firstly introduced and verified by the finite element method. Then, this algorithm is applied to simulate the contact between two semispherical asperities. The full sliding contact process is discretized into a series of transient steps, and each of these steps are calculated by the displacement-driven contact algorithm. The effects of the interference depth and the properties of, respectively, the tribofilm (thickness, elastic modulus, and yield strength) and the nanocrystalline layer on the sliding coefficient of friction are investigated. The results suggest that when surface adhesion and asperity damage are ignored, the plastic deformation of the tribofilm is the main source of the sliding friction. Greater interference depth, tribofilm with greater thickness, higher elastic modulus or lower yield strength, and the presence of a nanocrystalline layer will lead to a higher coefficient of friction in single asperity sliding.


Introduction
Sliding friction under boundary lubrication has been studied for almost a century [1][2][3][4][5].In the boundary lubrication regime, the lubricant film is not sufficiently thick to separate the surfaces and solid-to-solid contact will occur.Under these conditions, primarily, the highest peaks (asperities) on the surfaces will come into contact.Therefore, the asperities are carrying the load locally, so on a microscale, the stresses exceed the yield of the material causing plastic flow even though the macroscopic pressure is well below the yield limit.This elastoplastic deformation of asperities has been regarded as one of the important sources of friction [2].Meanwhile, a solid layer is also generated in boundary lubrication and plays a crucial role in protecting the system against severe friction and wear.This layer consists of two separate layers, one produced by chemical reactions and one by the mechanical energy dissipated at the top layer of the bulk material (nanocrystalline (NC) layer) (see Figure 1).The properties and composition of these layers are dominated by the lubricants used (mainly the additives package) and the conditions under which the system is running.Using this representation of the system in previous publications, wear rates were successfully predicted [6,7].In these publications, mechanical properties available in the literature were used for both the NC layer and the tribofilm.However, the coefficient of friction was presumed known and constant.In this publication, the latter will be investigated in more detail and the main question to be answered is: "Can plastic deformation of the tribofilm and NC layer explain the coefficient of friction observed in boundary lubrication?" Figure 1.Layers present in a sliding asperity contact system.The top layer is the tribofilm, which is a mixture of oxides and chemical products of the lubricant.The second layer is a nanocrystalline layer formed at the top of the bulk material by severe plastic deformation under high hydrostatic pressure and high shear rates.
It has been known for many years that different additives influence the tribological behavior of a system (the level of friction and wear).This can be explained by the fact that different additives will produce different tribofilms [8] which tweak the contact status (contact pressure and deformation, real contact area, etc.) of asperities.The finite element method (FEM) is a commonly used approach in solving contact problems.Liu et al. [9] analyzed the elastic-plastic contact of rough surfaces in linecontact problems.This work was extended by Etsion [10] for point contacts.Furthermore, Jackson et al. [11] developed a finite element model predicting the normal and tangential resisting forces of two sliding elastic, perfectly plastic frictionless spherical asperities.Mulvihill et al. [12] improved this work by imposing greater interference depths and considering surface adhesion in the contact.However, it should be noted that performing FEM contact analysis is time consuming since a very fine mesh is needed in contact areas to ensure convergence and accurate results, while the domain needs to be large enough to avoid edge effects.In addition, the contact condition itself (no penetration of the two bodies) is not inherently incorporated to the FEM method and requires additional degrees of freedom (DOF) for a solution to be reached.Often, the penalty function is used, in which the magnitude of the contact stiffness is tuned using a Hertzian reference contact.While the obtained value may be true for the reference case, it is very possible that in different scenarios the values obtained do not hold.
In addition to FEM, which is governed by continuum mechanics, alternative options have been proposed to solve the contact process.Molecular dynamics (MD), for example, simulates the contact process at the atomic level.This method allows all atoms inside the calculation domain to move simultaneously, and interactions between these atoms are calculated by Newton's law.It is thus possible to track the position information of each atom at a given time.Therefore, MD is particularly suitable for solving contact mechanics at the nanoscale [13][14][15].However, with increasing size and time scale, the computational cost will become unacceptably high and limit the application of MD in contact simulation.In recent decades, the discrete element method (DEM), as another alternative, has gained popularity in numerical tribology.DEM considers that the contact areas are composed of heterogeneous particles interacting with each other, governed by conservation and dynamics equations.DEM can calculate the behavior of detached particles efficiently and therefore is widely applied in simulating the third body behavior [16,17] and wear process [18,19].DEM became even more versatile by combining with the cellular automata method: a series of dry friction problems between different materials were successfully simulated [20][21][22].However, properly choosing the size of single particles and defining the interaction force between them are a challenge [23].Besides, Figure 1.Layers present in a sliding asperity contact system.The top layer is the tribofilm, which is a mixture of oxides and chemical products of the lubricant.The second layer is a nanocrystalline layer formed at the top of the bulk material by severe plastic deformation under high hydrostatic pressure and high shear rates.
It has been known for many years that different additives influence the tribological behavior of a system (the level of friction and wear).This can be explained by the fact that different additives will produce different tribofilms [8] which tweak the contact status (contact pressure and deformation, real contact area, etc.) of asperities.The finite element method (FEM) is a commonly used approach in solving contact problems.Liu et al. [9] analyzed the elastic-plastic contact of rough surfaces in line-contact problems.This work was extended by Etsion [10] for point contacts.Furthermore, Jackson et al. [11] developed a finite element model predicting the normal and tangential resisting forces of two sliding elastic, perfectly plastic frictionless spherical asperities.Mulvihill et al. [12] improved this work by imposing greater interference depths and considering surface adhesion in the contact.However, it should be noted that performing FEM contact analysis is time consuming since a very fine mesh is needed in contact areas to ensure convergence and accurate results, while the domain needs to be large enough to avoid edge effects.In addition, the contact condition itself (no penetration of the two bodies) is not inherently incorporated to the FEM method and requires additional degrees of freedom (DOF) for a solution to be reached.Often, the penalty function is used, in which the magnitude of the contact stiffness is tuned using a Hertzian reference contact.While the obtained value may be true for the reference case, it is very possible that in different scenarios the values obtained do not hold.
In addition to FEM, which is governed by continuum mechanics, alternative options have been proposed to solve the contact process.Molecular dynamics (MD), for example, simulates the contact process at the atomic level.This method allows all atoms inside the calculation domain to move simultaneously, and interactions between these atoms are calculated by Newton's law.It is thus possible to track the position information of each atom at a given time.Therefore, MD is particularly suitable for solving contact mechanics at the nanoscale [13][14][15].However, with increasing size and time scale, the computational cost will become unacceptably high and limit the application of MD in contact simulation.In recent decades, the discrete element method (DEM), as another alternative, has gained popularity in numerical tribology.DEM considers that the contact areas are composed of heterogeneous particles interacting with each other, governed by conservation and dynamics equations.DEM can calculate the behavior of detached particles efficiently and therefore is widely applied in simulating the third body behavior [16,17] and wear process [18,19].DEM became even more versatile by combining with the cellular automata method: a series of dry friction problems between different materials were successfully simulated [20][21][22].However, properly choosing the size of single particles and defining the interaction force between them are a challenge [23].Besides, the application of DEM is also constrained by computational cost in large scale and a longer time period.
The fast Fourier transform (FFT) technique and semianalytical method (SAM) provide an alternative for solving elastic and elastoplastic contact more efficiently.Various problems including single and layered material contact have been solved by SAM [24].Also, techniques have been developed to accelerate the calculation and improve the accuracy.Polonsky and Keer [25] proposed an iterative scheme by applying the conjugate gradient method.Liu et al. [26] eliminated the periodic errors in regular FFT by introducing the discrete fast Fourier transform (DC-FFT) technique.However, the aforementioned works only handled purely elastic contacts.To solve elastoplastic contact problems, Jacq et al. [27] presented a semianalytical method by considering plastic deformation as eigenstrains, based on Chiu's results [28,29].Recently, Wang et al. [30] acquired stress and displacement fields of eigenstrains in two perfectly bonded layered halfspaces, making it possible to solve the elastoplastic contact of layered material.
In this paper, the sliding friction of a single asperity contact including a tribofilm and a nanocrystalline layer is solved by applying the DC-FFT technique and SAM based on the work in [22].Consequently, this work extends the model proposed by Boucly et al. [31] which involved only homogeneous material, to a layered system.The plasticity in the layered structure is considered to follow the J2 flow theory [32] and is solved by the method in reference [30].Green's assumption on steady state sliding [33] is applied: two asperities will not approach or separate over each other while sliding, e.g., a predefined separation of the two surfaces is given.To get an initial indication of the influence of the different aspects affecting friction, the effects of mechanical properties of the tribofilm (elastic modulus, Poisson's ratio, yield strength, thickness), interference depth, and nanocrystalline layer are studied separately.The adhesion between tribofilms and damage to asperities are ignored, as the main goal of the study is to investigate the contribution of plastic deformation to friction.

Theory and Methodology
In this section, the general model is discussed first.The main backbone is a displacement-driven contact model for a layered elastoplastic system.This general model is used to investigate a single asperity contact with different mechanical properties.

Displacement-Driven Contact Model
First, a reference coordinate system is introduced, as shown in Figure 2.For simplicity but without loss of generality, a contact of an elastoplastic object loading on a rigid flat is considered, where: p is the contact pressure distribution built within the contact area; e p is the plastic strain in the elastoplastic object formed due to contact stress; u e is the surface displacement due to contact pressure, a function of the contact pressure distribution (p); u r is the surface displacement due to residual plastic strain, a function of plastic strains (e p ); h 0 is the initial surface separation relative to the initial undeformed geometry; δ z is the rigid body movement.The initial geometry (h 0 ) of contact objects (and their mechanical properties), together with a rigid body movement (δ z ), serves as input for the model.The model should be able to predict the contact pressure (p) and the subsurface plastic strains (e p ).
In the whole domain in which contact and noncontact areas both exist, the following equations and conditions hold: Inside the contact area: u e (p) + u r e p + h 0 = δ z (1) Outside the contact area: Equations ( 1) and ( 2) indicate that in contact areas, surface separation must be zero (no penetration).In noncontact areas, however, surface separation should be greater than zero, which is shown in Equations ( 3) and (4).
Since the ultimate goal of the model is to find the pressure distribution and plastic strains that make Equations ( 1)-( 4) hold, the exact expressions of u e (p) and u r (e p ) should be built first.This will be described in the following subsections.In the whole domain in which contact and noncontact areas both exist, the following equations and conditions hold: Inside the contact area: Outside the contact area: Equations ( 1) and ( 2) indicate that in contact areas, surface separation must be zero (no penetration).In noncontact areas, however, surface separation should be greater than zero, which is shown in Equations ( 3) and ( 4).
Since the ultimate goal of the model is to find the pressure distribution and plastic strains that make Equations ( 1)-( 4) hold, the exact expressions of ue(p) and ur(ep) should be built first.This will be described in the following subsections.

Elastic Displacement and Stress Fields
By employing Papkovich-Neuber potential and Fourier transformation, O'Sullivan and King [34] derived the displacement and stress fields in a layered material: with: Schematic of the contact.

Elastic Displacement and Stress Fields
By employing Papkovich-Neuber potential and Fourier transformation, O'Sullivan and King [34] derived the displacement and stress fields in a layered material: with: where k = 1 indicates coating material and k = 2 indicates substrate material.A, B, and C are coefficients to be determined by boundary conditions (see [34] for more details).
Once the frequency responses are acquired, the displacement and subsurface stress can then be linked to surface tractions by DC-FFT [26] Lubricants 2018, 6, 96 where q x is the surface traction in x direction, q y is the surface traction in y direction, and p is the normal pressure.The tilde sign on top of variables stands for 2D Fourier transformation.By letting q x and q y equal zero in Equation ( 7), the expressions of u e (p) can be acquired.

Residual Displacement and Stress Fields
The SAM was used in this work to evaluate the residual displacement and stress fields.In SAM, the plastic strains are considered as "eigenstrains", while their influence on the domain as a whole is solved by a combination of analytical and numerical approaches.The method proposed recently by Wang et al. [30] was used in this work and is described briefly below.
Firstly, a simple case is introduced, see Figure 3 [30].A full space is composed of two halfspaces I and II.These two halfspaces are perfectly bonded, while eigenstrains exist in halfspace I.

𝜎 ( ) = 𝐼𝐹𝐹𝑇(𝐶 ( ) . * 𝑝 + 𝐶 ( ) . * 𝑞 + 𝐶 ( ) . * 𝑞 ) (8)
where qx is the surface traction in x direction, qy is the surface traction in y direction, and p is the normal pressure.The tilde sign on top of variables stands for 2D Fourier transformation.By letting qx and qy equal zero in Equation ( 7), the expressions of ue(p) can be acquired.

Residual Displacement and Stress Fields
The SAM was used in this work to evaluate the residual displacement and stress fields.In SAM, the plastic strains are considered as "eigenstrains", while their influence on the domain as a whole is solved by a combination of analytical and numerical approaches.The method proposed recently by Wang et al. [30] was used in this work and is described briefly below.
Firstly, a simple case is introduced, see Figure 3 [30].A full space is composed of two halfspaces I and II.These two halfspaces are perfectly bonded, while eigenstrains exist in halfspace I.The stress and displacement field in both halfspaces are solved.For example, the displacement field in halfspace I can be expressed as The stress and displacement field in both halfspaces are solved.For example, the displacement field in halfspace I can be expressed as Readers are referred to [30] for more details on the influence matrices/coefficients.However, this simple case solution is not instantly applicable to the model because in layered structure contact problems, a halfspace is needed instead of a full space, as shown in Figure 4a.This halfspace is composed of a fully extended halfspace (the substrate) and a finite halfspace (the tribofilm, due to its finite thickness).
Readers are referred to [30] for more details on the influence matrices/coefficients.However, this simple case solution is not instantly applicable to the model because in layered structure contact problems, a halfspace is needed instead of a full space, as shown in Figure 4a.This halfspace is composed of a fully extended halfspace (the substrate) and a finite halfspace (the tribofilm, due to its finite thickness).The principle of superposition states that for a given small deformation problem domain, if the state I is a solution to the fundamental elasticity equations with prescribed body forces FI and surface tractions TI, and state II is a solution to the fundamental equations with prescribed body forces FII and surface tractions TII, then state I + II will be a solution to the problem with body forces FI + FII and surface tractions TI + TII.Therefore, the layered structure contact problem can be split into two subproblems: an eigenstrain problem in a full space plus a set of imaginary surface tractions to make the aggregate solution satisfy the halfspace boundary conditions (surface traction free), as seen in Figure 4b.
Note that the first subproblem is the eigenstrain problem in a full space, which can be solved by the method mentioned in Equation ( 9).The second subproblem, however, is a surface traction problem, which can be solved by Equations ( 7) and (8), by letting qx, qy, and p equal the imaginary surface tractions.By summing up the displacements caused by the eigenstrains and the imaginary stresses, the expressions of ur(ep) are acquired.

Model Route
With ue(p) and ur(ep) being expressed, the pressure distribution p and plastic strains ep are now explicitly linked to the inputs (rigid body movement δz and initial surface geometry h0) by Equations ( 1)-( 4).This is in essence an M*N linear system.However, it should be noted that the subsurface stress generated by the surface contact pressure causes plastic strains.These plastic strains cause residual displacement and stresses, which in turn alter the geometry and change the surface contact pressure.In other words, p and ep are coupled.Thus, a direct solution is difficult to acquire.
Therefore, an iteration loop is proposed to address the problem.First, an initial status is set.The pressure and subsurface stress are solved first, followed by a radial return algorithm to determine plastic deformation for the current loop.Then, the residual displacement is calculated and surface geometry is updated.Afterwards, the pressure distribution is recalculated based on the new surface geometry.This loop will keep running till a convergence in pressure distribution is achieved (see Figure 5).Note that in each iteration, an inner iteration (return mapping) is needed to map the stress back to yield surface and update plastic deformation.This process is shown in Figure 6.The principle of superposition states that for a given small deformation problem domain, if the state I is a solution to the fundamental elasticity equations with prescribed body forces F I and surface tractions T I , and state II is a solution to the fundamental equations with prescribed body forces F II and surface tractions T II , then state I + II will be a solution to the problem with body forces F I + F II and surface tractions T I + T II .Therefore, the layered structure contact problem can be split into two subproblems: an eigenstrain problem in a full space plus a set of imaginary surface tractions to make the aggregate solution satisfy the halfspace boundary conditions (surface traction free), as seen in Figure 4b.
Note that the first subproblem is the eigenstrain problem in a full space, which can be solved by the method mentioned in Equation ( 9).The second subproblem, however, is a surface traction problem, which can be solved by Equations ( 7) and (8), by letting q x , q y , and p equal the imaginary surface tractions.By summing up the displacements caused by the eigenstrains and the imaginary stresses, the expressions of u r (e p ) are acquired.

Model Route
With u e (p) and u r (e p ) being expressed, the pressure distribution p and plastic strains e p are now explicitly linked to the inputs (rigid body movement δ z and initial surface geometry h 0 ) by Equations ( 1)-( 4).This is in essence an M*N linear system.However, it should be noted that the subsurface stress generated by the surface contact pressure causes plastic strains.These plastic strains cause residual displacement and stresses, which in turn alter the geometry and change the surface contact pressure.In other words, p and e p are coupled.Thus, a direct solution is difficult to acquire.
Therefore, an iteration loop is proposed to address the problem.First, an initial status is set.The pressure and subsurface stress are solved first, followed by a radial return algorithm to determine plastic deformation for the current loop.Then, the residual displacement is calculated and surface geometry is updated.Afterwards, the pressure distribution is recalculated based on the new surface geometry.This loop will keep running till a convergence in pressure distribution is achieved (see Figure 5).Note that in each iteration, an inner iteration (return mapping) is needed to map the stress back to yield surface and update plastic deformation.This process is shown in Figure 6.Lubricants 2018, 6, 96 8 of 18

Single Asperity Sliding
In this work, the sliding asperity contact was idealized as a normal contact between two semispheres, following the work of Boucly [31] (see Figure 7).Using a semisphere to represent the asperity is a commonly adopted approach when studying a single asperity contact, as it allows a relatively simple calculation of the contact parameters.The two semispheres are aligned with the distance between two basis planes being held constant as h, while the interference depth is d 0 .Elastoplastic behavior of both tribofilm and substrate materials are considered, while adhesion between the surfaces is ignored.
The whole contact process is calculated in an incremental manner.In each incremental step, the corresponding rigid body movement is determined first, then the tangential and normal forces generated by the sliding contact, which are regarded as frictional and normal force, are calculated by the process proposed in Section 2.1.

Single Asperity Sliding
In this work, the sliding asperity contact was idealized as a normal contact between two semispheres, following the work of Boucly [31] (see Figure 7).Using a semisphere to represent the asperity is a commonly adopted approach when studying a single asperity contact, as it allows a relatively simple calculation of the contact parameters.The two semispheres are aligned with the distance between two basis planes being held constant as h, while the interference depth is d0.Elastoplastic behavior of both tribofilm and substrate materials are considered, while adhesion between the surfaces is ignored.
The whole contact process is calculated in an incremental manner.In each incremental step, the corresponding rigid body movement is determined first, then the tangential and normal forces generated by the sliding contact, which are regarded as frictional and normal force, are calculated by the process proposed in Section 2.1.First, the whole sliding contact process is discretized into n equally distributed incremental steps.In each of these incremental steps, the rigid body movement in the contact direction can be derived.In the ith step, note that in the triangle C2,original-C1-C2, the distance between the two centers of the semispheres is C1C2, Therefore, the rigid body movement is With  , calculated, the contact force in each step Fcontact,i can be calculated by the displacementdriven contact model.The frictional force and normal force in each step then can be acquired as tangential and vertical components of Fcontact,i: where θ is the angle between dΩ and the horizontal plane.
For the whole sliding process, the coefficient of friction (COF) is defined as .
Readers are referred to Boucly's work [31] for detailed expressions.

Results and Discussion
In this section, a tribosystem of a pair of steel asperities (25MnCr5) with a ZDDP tribofilm is set as a baseline, which is commonly seen as a soft-tribofilm and hard-substrate combination in industrial applications.The properties of the ZDDP film, steel substrate, and asperity radius were collected from the literature and are listed in Table 1.First, the whole sliding contact process is discretized into n equally distributed incremental steps.In each of these incremental steps, the rigid body movement in the contact direction can be derived.In the ith step, note that in the triangle C 2,original -C 1 -C 2 , the distance between the two centers of the semispheres is C 1 C 2 , Therefore, the rigid body movement is With δ z,i calculated, the contact force in each step F contact,i can be calculated by the displacement-driven contact model.The frictional force and normal force in each step then can be acquired as tangential and vertical components of F contact,i : where θ is the angle between dΩ and the horizontal plane.
For the whole sliding process, the coefficient of friction (COF) is defined as Readers are referred to Boucly's work [31] for detailed expressions.

Results and Discussion
In this section, a tribosystem of a pair of steel asperities (25MnCr5) with a ZDDP tribofilm is set as a baseline, which is commonly seen as a soft-tribofilm and hard-substrate combination in industrial applications.The properties of the ZDDP film, steel substrate, and asperity radius were collected from the literature and are listed in Table 1.Note: the Poission's ratio and yield strength of ZDDP film have not yet been reported in the literature.The yield strength above is estimated as one-third of the hardness in [35].Poission's ratio is assumed to be 0.3.
Next, a verification of displacement-driven contact model is presented, then the effects of different tribofilm properties on sliding friction are quantified.

Verification of Displacement-Driven Contact Model
The displacement-driven contact model proposed in Section 2.1 was verified by an asperity indented on a coated flat plane.For simplicity without losing generality, the asperity was assumed rigid, while both materials of the coated flat plane showed elastic-perfectly plastic behavior.The material properties and contact conditions are listed below in Table 2. Two different yield strengths of the coating material were used for comparison.The solution domain was meshed into an 80 × 80 × 80 grid (20 grid over the thickness direction of tribofilm).With a single element of size 6 × 6 × 6 nm, this resulted a calculation domain sized 480 × 480 × 480 nm, which was large enough to allow the plastic deformation to fully extend.The results of the SAM model were verified by comparing the solutions with those from commercial FEM software (ABAQUS, from Dassault Systèmes, Vélizy-Villacoublay, France).The FEM model used four-node bilinear axisymmetric quadrilateral with full integration (type CAX4) elements.The element sizes varied from near the contact area (indicating the volume which has this fine mesh, 6 × 6 nm) to 100 × 100 nm in faraway areas.To limit the effect of boundary conditions in the FEM simulation, the modelling domain had to be at least 10 times larger than the contact zone.This gives a total of 55,214 elements in the FEM model.
A comparison of the results is presented in Figure 8.The von Mises stress and equivalent plastic strain as a function of depth at the center of the contact are shown.Despite some minor differences, good agreement is observed.Contour plots of the von Mises stresses and equivalent plastic strain in the x-z plane are shown in Figures 9 and 10.The results show that plastic strains solely occur in the tribofilm due to its lower yield strength than steel.For the 500 MPa case, the plastic strain tends to concentrate near the bonding interface, while in the 800 MPa case, the plastic strain is more homogeneously distributed.Contour plots of the von Mises stresses and equivalent plastic strain in the x-z plane are shown in Figures 9 and 10.The results show that plastic strains solely occur in the tribofilm due to its lower yield strength than steel.For the 500 MPa case, the plastic strain tends to concentrate near the bonding interface, while in the 800 MPa case, the plastic strain is more homogeneously distributed.Contour plots of the von Mises stresses and equivalent plastic strain in the x-z plane are shown in Figures 9 and 10.The results show that plastic strains solely occur in the tribofilm due to its lower yield strength than steel.For the 500 MPa case, the plastic strain tends to concentrate near the bonding interface, while in the 800 MPa case, the plastic strain is more homogeneously distributed.Contour plots of the von Mises stresses and equivalent plastic strain in the x-z plane are shown in Figures 9 and 10.The results show that plastic strains solely occur in the tribofilm due to its lower yield strength than steel.For the 500 MPa case, the plastic strain tends to concentrate near the bonding interface, while in the 800 MPa case, the plastic strain is more homogeneously distributed.

Single Asperity Sliding
Without considering adhesion effects, the sole mechanism behind friction in single asperity sliding is plastic deformation.The tangential force during the simulated sliding process is shown in Figure 11 for both an elastic and an elastoplastic contact.
In the pure elastic deformation case, the friction force curve acts in a centrosymmetric manner, and thus the net frictional force is zero.However, when plastic deformation is involved, asymmetry of the tangential force can be observed.The reason for this change is that residual displacements modify the surface geometry and break the symmetry of the system.In other words, a certain amount of energy is dissipated during plastic deformation.
Therefore, anything that affects the plastic behavior of the system affects friction.In the following parts, effects of interference depth, tribofilm thickness, elastic/plastic properties of the tribofilm, and the presence of an NC layer on sliding friction are presented and discussed.
Also, the maximum calculated plastic strain in the tribofilm is introduced to measure the risk of tribofilm wear.As indicated in [7] in single asperity sliding, most plastic strain occurs in the tribofilm.The higher the maximum plastic strain is, the more the tribofilm is worn.

Single Asperity Sliding
Without considering adhesion effects, the sole mechanism behind friction in single asperity sliding is plastic deformation.The tangential force during the simulated sliding process is shown in Figure 11 for both an elastic and an elastoplastic contact.
In the pure elastic deformation case, the friction force curve acts in a centrosymmetric manner, and thus the net frictional force is zero.However, when plastic deformation is involved, asymmetry of the tangential force can be observed.The reason for this change is that residual displacements modify the surface geometry and break the symmetry of the system.In other words, a certain amount of energy is dissipated during plastic deformation.
Therefore, anything that affects the plastic behavior of the system affects friction.In the following parts, effects of interference depth, tribofilm thickness, elastic/plastic properties of the tribofilm, and the presence of an NC layer on sliding friction are presented and discussed.
Also, the maximum calculated plastic strain in the tribofilm is introduced to measure the risk of tribofilm wear.As indicated in [7] in single asperity sliding, most plastic strain occurs in the tribofilm.The higher the maximum plastic strain is, the more the tribofilm is worn.

Effects of Interference Depth
A series of interference depths were applied (20,35,50, and 65 nm), while other parameters were held constant (see Table 1).The results are shown in Figure 12, where it becomes clear that the tangential force increases with interference fit, which makes sense as more material is plastically deformed.This is also reflected by the increase in calculated coefficient of friction.It should be noted that the maximum plastic strain in the SAM calculations should not exceed 10% to ensure stable results [37].The increase in plastic deformation of the tribofilm also indicates that an increase in wear could be expected for the system under consideration, suggesting a relation between friction and wear.

Effects of Interference Depth
A series of interference depths were applied (20,35,50, and 65 nm), while other parameters were held constant (see Table 1).The results are shown in Figure 12, where it becomes clear that the tangential force increases with interference fit, which makes sense as more material is plastically deformed.This is also reflected by the increase in calculated coefficient of friction.It should be noted that the maximum plastic strain in the SAM calculations should not exceed 10% to ensure stable results [37].The increase in plastic deformation of the tribofilm also indicates that an increase in wear could be expected for the system under consideration, suggesting a relation between friction and wear.

Effects of Tribofilm Thickness
To investigate the effect of the thickness of the tribofilm, four thicknesses were applied in the model: 10, 60, 80, and 100 nm, respectively.The results are shown in Figure 13.

Effects of Tribofilm Thickness
To investigate the effect of the thickness of the tribofilm, four thicknesses were applied in the model: 10, 60, 80, and 100 nm, respectively.The results are shown in Figure 13.

Effects of Tribofilm Thickness
To investigate the effect of the thickness of the tribofilm, four thicknesses were applied in the model: 10, 60, 80, and 100 nm, respectively.The results are shown in Figure 13.The normal contact force decreases with thicker tribofilm, while softer tribofilms act like a cushion.The thicker the cushion is, the greater the effect will be.The friction increases with thicker tribofilm, which shows an almost linear relationship.This can be understood in terms of plastic deformation and energy dissipation.However, the relationship between the amount of plasticity and friction is highly nonlinear, showing a greater increase in plasticity than in friction.

Effects of Elastic Modulus
The elastic modulus of the tribofilms were chosen from 20 to 65 GPa, which refer to tribofilms softer and stiffer than ZDDP (e.g., the reference case, i.e., 35 GPa).
The results shown in Figure 14 show that the frictional and normal forces increase with a stiffer tribofilm than expected.Also, a tribofilm with higher elastic modulus leads to a higher coefficient of friction.This is because with a higher elastic modulus, the tribofilm reaches its elastic limit at a lower strain level, which leads to a greater volume of plasticity.Additionally, a higher elastic modulus means more energy is needed for the same amount of plastic deformation.So for these two reasons, more energy is absorbed by the system, thus generating greater friction.The normal contact force decreases with thicker tribofilm, while softer tribofilms act like a cushion.The thicker the cushion is, the greater the effect will be.The friction increases with thicker tribofilm, which shows an almost linear relationship.This can be understood in terms of plastic deformation and energy dissipation.However, the relationship between the amount of plasticity and friction is highly nonlinear, showing a greater increase in plasticity than in friction.

Effects of Elastic Modulus
The elastic modulus of the tribofilms were chosen from 20 to 65 GPa, which refer to tribofilms softer and stiffer than ZDDP (e.g., the reference case, i.e., 35 GPa).
The results shown in Figure 14 show that the frictional and normal forces increase with a stiffer tribofilm than expected.Also, a tribofilm with higher elastic modulus leads to a higher coefficient of friction.This is because with a higher elastic modulus, the tribofilm reaches its elastic limit at a lower strain level, which leads to a greater volume of plasticity.Additionally, a higher elastic modulus means more energy is needed for the same amount of plastic deformation.So for these two reasons, more energy is absorbed by the system, thus generating greater friction.The normal contact force decreases with thicker tribofilm, while softer tribofilms act like a cushion.The thicker the cushion is, the greater the effect will be.The friction increases with thicker tribofilm, which shows an almost linear relationship.This can be understood in terms of plastic deformation and energy dissipation.However, the relationship between the amount of plasticity and friction is highly nonlinear, showing a greater increase in plasticity than in friction.

Effects of Elastic Modulus
The elastic modulus of the tribofilms were chosen from 20 to 65 GPa, which refer to tribofilms softer and stiffer than ZDDP (e.g., the reference case, i.e., 35 GPa).
The results shown in Figure 14 show that the frictional and normal forces increase with a stiffer tribofilm than expected.Also, a tribofilm with higher elastic modulus leads to a higher coefficient of friction.This is because with a higher elastic modulus, the tribofilm reaches its elastic limit at a lower strain level, which leads to a greater volume of plasticity.Additionally, a higher elastic modulus means more energy is needed for the same amount of plastic deformation.So for these two reasons, more energy is absorbed by the system, thus generating greater friction.

Effects of Tribofilm Yield Strength
In this part, the effect of the yield stress of the tribofilm on sliding friction is quantified.The yield stresses of tribofilms were set in a range from 350 to 800 MPa, which covers the tribofilms with higher and lower yield strengths than the ZDDP tribofilm.
The results are shown in Figure 15.

Effects of Tribofilm Yield Strength
In this part, the effect of the yield stress of the tribofilm on sliding friction is quantified.The yield stresses of tribofilms were set in a range from 350 to 800 MPa, which covers tribofilms with higher and lower yield strengths than the ZDDP tribofilm.
The results are shown in Figure 15.

Effects of Tribofilm Yield Strength
In this part, the effect of the yield stress of the tribofilm on sliding friction is quantified.The yield stresses of tribofilms were set in a range from 350 to 800 MPa, which covers the tribofilms with higher and lower yield strengths than the ZDDP tribofilm.
The results are shown in Figure 15.From the results, it can be concluded that a tribofilm with a higher yield strength will cause a higher normal force, which is due to less plastic strain/displacement and thus a smaller contact area.Also, the coefficient of friction will decrease due to a higher yield strength of the tribofilm resulting from less plastic deformation.Therefore, the system will be affected less and show lower friction.

Effects of Presence of the NC Layer
In this subsection, effect of the yield stress of the NC layer to sliding friction is quantified.In order to do so, the presence of a nanocrystalline layer is considered.The hardness of the nanocrystalline layer was taken from sample B in [38], see Figure 16.From the results, it can be concluded that a tribofilm with a higher yield strength will cause a higher normal force, which is due to less plastic strain/displacement and thus a smaller contact area.Also, the coefficient of friction will decrease due to a higher yield strength of the tribofilm resulting from less plastic deformation.Therefore, the system will be affected less and show lower friction.

Effects of Presence of the NC Layer
In this subsection, the effect of the yield stress of the NC layer to sliding friction is quantified.In order to do so, the presence of a nanocrystalline layer is considered.The hardness of the nanocrystalline layer was taken from sample B in [38], see Figure 16.The results are shown in Figure 17.With the presence of an NC layer, the frictional and normal contact forces increase.The reason for this effect is similar with a higher yield strength of the tribofilm leading to a higher normal contact force: the less strong NC layer allows higher plastic strain and displacement.Therefore, the geometry can better accommodate the contact deformation.The results are shown in Figure 17.
With the presence of an NC layer, the frictional and normal contact forces increase.The reason for this effect is similar with a higher yield strength of the tribofilm leading to a higher normal contact force: the less strong NC layer allows higher plastic strain and displacement.Therefore, the geometry can better accommodate the contact deformation.From the results, it can be concluded that a tribofilm with a higher yield strength will cause a higher normal force, which is due to less plastic strain/displacement and thus a smaller contact area.Also, the coefficient of friction will decrease due to a higher yield strength of the tribofilm resulting from less plastic deformation.Therefore, the system will be affected less and show lower friction.

Effects of Presence of the NC Layer
In this subsection, the effect of the yield stress of the NC layer to sliding friction is quantified.In order to do so, the presence of a nanocrystalline layer is considered.The hardness of the nanocrystalline layer was taken from sample B in [38], see Figure 16.The results are shown in Figure 17.With the presence of an NC layer, the frictional and normal contact forces increase.The reason for this effect is similar with a higher yield strength of the tribofilm leading to a higher normal contact force: the less strong NC layer allows higher plastic strain and displacement.Therefore, the geometry can better accommodate the contact deformation.

Figure 2 .
Figure 2. Schematic of the contact.

Figure 3 .
Figure 3. Schematic of two joined elastic halfspaces inclusion problem, where two halfspaces are perfectly bonded together with inclusions in halfspace I [30].

Figure 3 .
Figure 3. Schematic of two joined elastic halfspaces inclusion problem, where two halfspaces are perfectly bonded together with inclusions in halfspace I [30].
(a) Schematic of layered halfspace and eigenstrains Lubricants 2018, 6, x FOR PEER REVIEW 7 of 19 (b) Layered contact problem subproblems I and II

Figure 4 .
Figure 4. Schematic of the superposition method for solving a layered structure [30].

Figure 4 .
Figure 4. Schematic of the superposition method for solving a layered structure [30].

Figure 11 .
Figure 11.Tangential force curve for pure elastic and elastoplastic sliding asperity contacts.

Figure 11 .
Figure 11.Tangential force curve for pure elastic and elastoplastic sliding asperity contacts.

Table 2 .
Parameters of the indentation model.