On an elastoplastic sliding model for a coated single asperity

(1)

lubricants

Article

On an Elastoplastic Sliding Model for a Coated

Single Asperity

Can Wang * and Dik J. Schipper

Laboratory for Surface Technology and Tribology, Department of Engineering Technology, University of Twente, P.O. Box 217, 7500 AE Enschede, The Netherlands; d.j.schipper@utwente.nl

* Correspondence: c.wang-1@utwente.nl; Tel.: +31-53-489-2663

Received: 11 September 2018; Accepted: 22 October 2018; Published: 1 November 2018

 Abstract: 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.

Keywords:boundary lubrication; friction model; tribofilm

1. Introduction

Sliding friction under boundary lubrication has been studied for almost a century [1–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 Figure1). 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?”

(2)

Lubricants 2018, 6, 96 2 of 18

Lubricants 2018, 6, x FOR PEER REVIEW 2 of 19

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 line-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–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–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–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–22]. However, properly choosing the size of single particles

(3)

Lubricants 2018, 6, 96 3 of 18

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.

2. 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.

2.1. Displacement-Driven Contact Model

First, a reference coordinate system is introduced, as shown in Figure2. 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;

epis the plastic strain in the elastoplastic object formed due to contact stress;

ue is the surface displacement due to contact pressure, a function of the contact pressure

distribution (p);

uris the surface displacement due to residual plastic strain, a function of plastic strains (ep);

h0is the initial surface separation relative to the initial undeformed geometry;

δzis the rigid body movement.

The initial geometry (h0) 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 (ep).

In the whole domain in which contact and noncontact areas both exist, the following equations and conditions hold:

Inside the contact area:

ue(p) +ur ep+h0=δz (1)

(4)

Lubricants 2018, 6, 96 4 of 18

Outside the contact area:

ue(p) +ur ep+h0≥δz (3)

p=0 (4)

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._{Lubricants 2018, 6, x FOR PEER REVIEW} _{4 of 19}

Figure 2. Schematic of the contact.

In the whole domain in which contact and noncontact areas both exist, the following equations and conditions hold:

Inside the contact area:

𝑢 (𝑝) + 𝑢 𝑒 + ℎ = 𝛿 (1)

𝑝 > 0 (2)

Outside the contact area:

𝑢 (𝑝) + 𝑢 𝑒 + ℎ ≥ 𝛿 (3)

𝑝 = 0 (4)

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.

2.1.1. 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:

𝑢( )= 𝐹𝑇[𝜑_,( )+ 𝑥𝜓( )_, + 𝑧 𝜓( )_, − (3 − 4𝜈 )𝜓( )] (5) 𝜎( )= 𝐹𝑇[𝜑,( )− 2𝜈 𝜓( ), + 𝜓( ), 𝛿 − (1 − 2𝜈 ) 𝜓( ), + 𝜓( ), + 𝑥𝜓( ), + 𝑧 𝜓( ), ] (6) 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]

Figure 2.Schematic of the contact. 2.1.1. 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:

e u(k)_i = 1 2Gk FThϕ(k)_,i +xψ_1,i(k)+zkψ(k)_3,i − (3−4νk)ψ_i(k) i (5) e σ_i(k)=FT h ϕ(k)_,ij −2νk ψ_1,1(k)+ψ_3,3(k) δij− (1−2νk) ψ_i,j(k)+ψ(k)_j,i +xψ_1,ij(k)+zkψ_3,ij(k) i (6) with: e ϕ(k)=A(k)e−αzk+A(k)eαzk e ψ(k)₁ =B(k)e−αzk+B(k)eαzk e ψ₃(k)=C(k)e−αzk+C(k)eαzk

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]

   uxe uye uz_e   = IFFT         e Cux qx Ceuqyx Ceupx e Cquxy Ce uy qy Ce uy p e Cuz qx Cequyz Ceupz   .∗    e qx e qy e p         (7) σe(k)=IFFT e Cσ(k) p .∗ep+Ce σ(k) qx .∗qex+Ce σ(k) qy .∗qey (8)

(5)

Lubricants 2018, 6, 96 5 of 18

where qxis the surface traction in x direction, qyis 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 qxand qyequal zero in Equation (7), the expressions of ue(p) can be acquired.

2.1.2. 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 Figure3[30]. A full space is composed of two halfspaces I and II. These two halfspaces are perfectly bonded, while eigenstrains exist in halfspace I.

𝑢 𝑢 𝑢 = 𝐼𝐹𝐹𝑇 ⎩ ⎨ ⎧ ⎣ ⎢ ⎢ ⎡𝐶 𝐶 𝐶 𝐶 𝐶 𝐶 𝐶 𝐶 _{𝐶 ⎦}⎥ ⎥ ⎤ .∗ 𝑞 𝑞 𝑝 ⎭ ⎬ ⎫ (7) 𝜎( )= 𝐼𝐹𝐹𝑇(𝐶 ( ).∗ 𝑝 + 𝐶 ( ).∗ 𝑞 + 𝐶 ( ).∗ 𝑞 ) (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.

2.1.2. 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.

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

The stress and displacement field in both halfspaces are solved. For example, the displacement field in halfspace I can be expressed as

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

The stress and displacement field in both halfspaces are solved. For example, the displacement field in halfspace I can be expressed as

2µu_ijI,out(x1α, x2β, x3γ)/C = ∑L ν=1 M ∑ η=1 N ∑ ξ=1 U_ijkl(0)(x1α−x1ξ, x2β−x2η, x3γ −x3ν)ekl(x1ξ, x2η, x3ν) + ∑L ν=1 M ∑ η=1 N ∑ ξ=1 U_ijkl(1)(x1α−x1ξ, x2β−x2η, x3γ +x3ν)ekl(x1ξ, x2η, x3ν) +x3ν L ∑ ν=1 M ∑ η=1 N ∑ ξ=1 U_ijkl(5)(x1α−x1ξ, x2β−x2η, x3γ +x3ν)ekl(x1ξ, x2η, x3ν) +x2_3ν ∑L ν=1 M ∑ η=1 N ∑ ξ=1 U_ijkl(5)(x1α−x1ξ, x2β−x2η, x3γ +x3ν)ekl(x1ξ, x2η, x3ν) (1≤α≤M, 1≤β≤N, 1≤γ≤L1). (9)

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).

(6)

Lubricants 2018, 6, 96 6 of 18

2𝜇𝑢, 𝑥_𝛼, 𝑥 , 𝑥 /𝐶 = 𝑈( ) 𝑥 𝛼− 𝑥 , 𝑥 − 𝑥 , 𝑥 − 𝑥 𝑒 𝑥 , 𝑥 , 𝑥 + 𝑈( ) 𝑥𝛼− 𝑥 , 𝑥 − 𝑥 , 𝑥 + 𝑥 𝑒 𝑥 , 𝑥 , 𝑥 + 𝑥 𝑈( ) 𝑥𝛼− 𝑥 , 𝑥 − 𝑥 , 𝑥 + 𝑥 𝑒 𝑥 , 𝑥 , 𝑥 + 𝑥 𝑈( ) 𝑥𝛼− 𝑥 , 𝑥 − 𝑥 , 𝑥 + 𝑥 𝑒 𝑥 , 𝑥 , 𝑥 (1 ≤

𝛼

≤ 𝑀, 1 ≤ 𝛽 ≤ 𝑁, 1 ≤ 𝛾 ≤ 𝐿 ). (9)

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).

(a) Schematic of layered halfspace and eigenstrains

(b) Layered contact problem subproblems I and II

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

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.

2.1.3. Model Route

With ue(p) and ur(ep) being expressed, the pressure distribution pand 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.

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

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 FIand 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 Figure4b.

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.

2.1.3. 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 epare 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

(7)

Lubricants 2018, 6, 96 7 of 18

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 Figure5). 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 Figure6.

Figure 5. Flowchart of iteration loop.

Figure 6. Inner plastic loop. Figure 5.Flowchart of iteration loop.

Figure 5. Flowchart of iteration loop.

Figure 6. Inner plastic loop. Figure 6.Inner plastic loop.

(8)

Lubricants 2018, 6, 96 8 of 18

2.2. 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 Section2.1.

2.2. 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.

Figure 7. Sliding contact between two asperities [20].

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

𝛿, = 2 ∗ 𝑅 − 𝐶 𝐶 ,.

With 𝛿, calculated, the contact force in each step Fcontact,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 Fcontact,i:

𝐹

_{_}

= 𝐹

_{, _}

=

𝑝 ∗ 𝑑

Ω

Ω ∈ contact area

∗ 𝑐𝑜𝑠𝜃

𝐹

_

= 𝐹

, _

=

𝑝 ∗ 𝑑

Ω ∗

𝑠𝑖𝑛𝜃

Ω ∈ contact area

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.

3. 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.

Figure 7.Sliding contact between two asperities [20].

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

δz,i =2 ∗R − C1C2,i.

With δz,i calculated, the contact force in each step Fcontact,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 Fcontact,i:

Ff riction_i =Fcontact, i_ x= x Ω ∈ contact area p∗dΩ∗cosθ Fnormal_i=Fcontact,i_y = x Ω ∈ contact area p∗dΩ∗sinθ where θ is the angle between dΩ and the horizontal plane.

For the whole sliding process, the coefficient of friction (COF) is defined as COF= ∑ Ff riction_i

∑ Fnormal_i

= ∑ Fcontact, i_ x ∑ Fcontact,i_y

.

(9)

Lubricants 2018, 6, 96 9 of 18

3. 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 Table1.

Table 1.Properties of ZDDP tribofilm and steel asperities.

Parameter Values

Young’s modulus, EZDDP(GPa) 35 [35]

Poisson’s ratio, νZDDP 0.3

Yield stress, σZDDP(MPa) 500 [35]

Young’s modulus, Easperity(GPa) 205

Poisson’s ratio, νasperity 0.3

Yield stress, σasperity(MPa) 1200

Radius of asperity (µm) 10 [36]

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.

3.1. Verification of Displacement-Driven Contact Model

The displacement-driven contact model proposed in Section2.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 Table2. 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.

Table 2.Parameters of the indentation model.

Parameter Values

Young’s modulus, Ecoating(GPa) 35

Poisson’s ratio, νcoating 0.3

Yield stress, σcoating(MPa) 500, 800

Young’s modulus, Esubstrate(GPa) 205

Poisson’s ratio, νsubstrate 0.3

Yield stress, σsubstrate(MPa) 1200

Radius of asperity (µm) 10 Indentation depth (nm) 10

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 Figure8. 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.

(10)

Lubricants 2018, 6, 96Lubricants 2018, 6, x FOR PEER REVIEW 10 of 1811 of 19

(a) (b)

Figure 8. Comparison of finite element method (FEM) and semianalytical method (SAM) results along the z axis (x = 0, y = 0): (a) von Mises stress; (b) equivalent plastic strain.

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.

(a) _(b)

Figure 9. Contour of von Mises stress: (a) σtribofilm = 500 MPa; (b) σtribofilm = 800 MPa.

(a) (b)

Figure 10. Contour of equivalent plastic strain: (a) σtribofilm = 500 MPa; (b) σtribofilm = 800 MPa.

Figure 8.Comparison of finite element method (FEM) and semianalytical method (SAM) results along the z axis (x = 0, y = 0): (a) von Mises stress; (b) equivalent plastic strain.

Contour plots of the von Mises stresses and equivalent plastic strain in the x-z plane are shown in Figures9and10. 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.

(a) (b)

(a) _(b)

(a) (b)

Figure 9.Contour of von Mises stress: (a) σtribofilm= 500 MPa; (b) σtribofilm= 800 MPa.

(a) (b)

(a) _(b)

(a) (b)

(11)

Lubricants 2018, 6, 96 11 of 18

3.2. 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 Figure11 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.

3.2. 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.

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

3.2.1. 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.

Figure 11.Tangential force curve for pure elastic and elastoplastic sliding asperity contacts. 3.2.1. Effects of Interference Depth

A series of interference depths were applied (20, 35, 50, and 65 nm), while other parameters were held constant (see Table1). The results are shown in Figure12, 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.

(12)

Lubricants 2018, 6, 96 12 of 18

(a) _(b)

(c) (d)

Figure 12. Results for different interference depth: (a) frictional force; (b) normal force; (c) coefficient of friction; (d) max. equivalent plastic strain.

3.2.2. 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.

(a) _(b)

Figure 12.Results for different interference depth: (a) frictional force; (b) normal force; (c) coefficient of friction; (d) max. equivalent plastic strain.

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 Figure13.

(a) _(b)

(c) (d)

Figure 12. Results for different interference depth: (a) frictional force; (b) normal force; (c) coefficient of friction; (d) max. equivalent plastic strain.

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.

(a) _(b)

(13)

Lubricants 2018, 6, 96 13 of 18

(c) (d)

Figure 13. Results for different tribofilm thickness: (a) frictional force; (b) normal force; (c) coefficient of friction; (d) max. equivalent plastic strain.

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.

3.2.3. 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.

(a) _(b)

Figure 13.Results for different tribofilm thickness: (a) frictional force; (b) normal force; (c) coefficient of friction; (d) max. equivalent plastic strain.

The results shown in Figure14show 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.

(c) (d)

Figure 13. Results for different tribofilm thickness: (a) frictional force; (b) normal force; (c) coefficient of friction; (d) max. equivalent plastic strain.

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.

(a) _(b)

(14)

Lubricants 2018, 6, 96Lubricants 2018, 6, x FOR PEER REVIEW 15 of 19 14 of 18

(c) (d)

Figure 14. Results for different tribofilm elastic modulus: (a) frictional force; (b) normal force; (c) coefficient of friction; (d) max. equivalent plastic strain.

3.2.4. 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.

(a) (b)

(c)

(d)

Figure 14. Results for different tribofilm elastic modulus: (a) frictional force; (b) normal force; (c) coefficient of friction; (d) max. equivalent plastic strain.

The results are shown in Figure15.

(c) (d)

Figure 14. Results for different tribofilm elastic modulus: (a) frictional force; (b) normal force; (c) coefficient of friction; (d) max. equivalent plastic strain.

(a) (b)

(c)

(d)

Figure 15. Results for different tribofilm yield strengths (a) frictional force; (b) normal force; (c) coefficient of friction; (d) max. equivalent plastic strain.

(15)

Lubricants 2018, 6, 96 15 of 18

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. 3.2.5. 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 Figure16.

Figure 15. Results for different tribofilm yield strengths (a) frictional force; (b) normal force; (c) coefficient of friction; (d) max. equivalent plastic strain.

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.

Figure 16. Hardness of nanocrystalline (NC) layer at different depths, reproduced from Figure 4 in [38].

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.

(a) (b)

Figure 16.Hardness of nanocrystalline (NC) layer at different depths, reproduced from Figure4in [38]. The results are shown in Figure17.

Figure 15. Results for different tribofilm yield strengths (a) frictional force; (b) normal force; (c) coefficient of friction; (d) max. equivalent plastic strain.

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.

Figure 16. Hardness of nanocrystalline (NC) layer at different depths, reproduced from Figure 4 in [38].

(a) (b)

(16)

Lubricants 2018, 6, 96Lubricants 2018, 6, x FOR PEER REVIEW 16 of 1817 of 19

(c) _(d)

Figure 17. Results for presence of NC layer: (a) frictional force; (b) normal force; (c) coefficient of friction; (d) max. equivalent plastic strain.

The coefficient of friction decreases with a higher yield strength of the NC layer. The relatively softer NC layer undergoes higher plastic strain and therefore dissipates more energy, which is consumed as friction. The maximum plastic strains in the tribofilm decrease with the presence of an NC layer because a softer NC layer makes it more possible for the plastic strain to extend downwards into the substrate. This makes the plastic strain accumulate less inside the tribofilm, which reduces the risk of tribofilm damage.

4. Conclusions

As in a previous study [7] it was shown that mild wear in the boundary lubrication regime may very well be governed by the plastic deformation of the topmost material. The question is whether friction is also dominated by the mechanism. To answer this question, an elastoplastic coated single asperity sliding model was developed, with surface adhesion being ignored. The sliding process was solved in an incremental manner based on a displacement-driven contact model, accelerated by applying the conjugate gradient method (GCM) and DC-FFT technique. Analyses were carried out at various interference depths and tribofilm parameters, namely, thickness, elastic modulus, and yield strength. The main conclusion is that with adhesion being ignored, to obtain realistic values (~0.1) regarding the coefficient of friction, an unrealistically great interference depth would be needed and cause very high plastic strains. This would, however, in practice result in unrealistic values for wear. Thus, in layered boundary lubrication, other sources, e.g., surface adhesion, ploughing, and removal of asperities [39,40] may play a role in sliding friction. It is reported that with asperity failure included, the coefficient of friction is slightly higher than when the asperity failure ignored [12]. These effects, however, are beyond the scope of this paper. Nevertheless, the following insights can be obtained from the results:

(1) The frictional force in a single asperity sliding is closely related to the plastic deformation of the asperities: greater plastic strain, higher sliding friction.

(2) Interference depth, tribofilm thickness, and the mechanical properties of the tribofilm will all affect the frictional force. Higher interference depth, thicker tribofilm, lower yield strength, and higher elastic modulus will contribute to a higher frictional force.

(3) The presence of a nanocrystalline layer will alter the behavior of the tribosystem. The presence of an NC layer will increase the coefficient of friction and ease the plastic strain accumulating in the tribofilm.

Author Contributions: Conceptualization, C.W. and D.J.S.; Methodology, C.W. and D.J.S.; Software, C.W.; Validation, C.W.; Formal Analysis, C.W.; Investigation, C.W.; Resources, C.W.; Data Curation, C.W.; Writing-Original Draft Preparation, C.W.; Writing-Review & Editing, C.W. and D.J.S.; Visualization, C.W.; Supervision, D.J.S.; Project Administration, D.J.S.; Funding Acquisition, D.J.S.

Figure 17. Results for presence of NC layer: (a) frictional force; (b) normal force; (c) coefficient of friction; (d) max. equivalent plastic strain.

The coefficient of friction decreases with a higher yield strength of the NC layer. The relatively softer NC layer undergoes higher plastic strain and therefore dissipates more energy, which is consumed as friction. The maximum plastic strains in the tribofilm decrease with the presence of an NC layer because a softer NC layer makes it more possible for the plastic strain to extend downwards into the substrate. This makes the plastic strain accumulate less inside the tribofilm, which reduces the risk of tribofilm damage.

4. Conclusions

As in a previous study [7] it was shown that mild wear in the boundary lubrication regime may very well be governed by the plastic deformation of the topmost material. The question is whether friction is also dominated by the mechanism. To answer this question, an elastoplastic coated single asperity sliding model was developed, with surface adhesion being ignored. The sliding process was solved in an incremental manner based on a displacement-driven contact model, accelerated by applying the conjugate gradient method (GCM) and DC-FFT technique. Analyses were carried out at various interference depths and tribofilm parameters, namely, thickness, elastic modulus, and yield strength. The main conclusion is that with adhesion being ignored, to obtain realistic values (~0.1) regarding the coefficient of friction, an unrealistically great interference depth would be needed and cause very high plastic strains. This would, however, in practice result in unrealistic values for wear. Thus, in layered boundary lubrication, other sources, e.g., surface adhesion, ploughing, and removal of asperities [39,40] may play a role in sliding friction. It is reported that with asperity failure included, the coefficient of friction is slightly higher than when the asperity failure ignored [12]. These effects, however, are beyond the scope of this paper. Nevertheless, the following insights can be obtained from the results:

(1) The frictional force in a single asperity sliding is closely related to the plastic deformation of the asperities: greater plastic strain, higher sliding friction.

(2) Interference depth, tribofilm thickness, and the mechanical properties of the tribofilm will all affect the frictional force. Higher interference depth, thicker tribofilm, lower yield strength, and higher elastic modulus will contribute to a higher frictional force.

(3) The presence of a nanocrystalline layer will alter the behavior of the tribosystem. The presence of an NC layer will increase the coefficient of friction and ease the plastic strain accumulating in the tribofilm.

Author Contributions: Conceptualization, C.W. and D.J.S.; Methodology, C.W. and D.J.S.; Software, C.W.; Validation, C.W.; Formal Analysis, C.W.; Investigation, C.W.; Resources, C.W.; Data Curation, C.W.; Writing-Original Draft Preparation, C.W.; Writing-Review & Editing, C.W. and D.J.S.; Visualization, C.W.; Supervision, D.J.S.; Project Administration, D.J.S.; Funding Acquisition, D.J.S.

(17)

Lubricants 2018, 6, 96 17 of 18

Funding:The authors thank Bosch Transmission technology b.v. for the financial support. Conflicts of Interest:The authors declare no conflict of interest.

References

1. Hardy, W.B.; Doubleday, I. Boundary lubrication—The paraffin series. Proc. R. Soc. Lond. A 1922, 100, 550–574. [CrossRef]

2. Komvopoulos, K.N.N.P.; Saka, N.; Suh, N.P. The mechanism of friction in boundary lubrication. J. Tribol. 1985, 107, 452–462. [CrossRef]

3. Persson, B.N.J. Theory of friction and boundary lubrication. Phys. Rev. B 1993, 48, 18140. [CrossRef] 4. Dupont, P.E.; Dunlap, E.P. Friction modeling and control in boundary lubrication. In Proceedings of the

American Control Conference, San Francisco, CA, USA, 2–4 June 1993; pp. 1910–1915. 5. Spikes, H. The history and mechanisms of ZDDP. Tribol. Lett. 2004, 17, 469–489. [CrossRef]

6. Bosman, R.; Schipper, D.J. Mild wear prediction of boundary-lubricated contacts. Tribol. Lett. 2011, 42, 169–178. [CrossRef]

7. Bosman, R.; Schipper, D.J. Mild wear maps for boundary lubricated contacts. Wear 2012, 280, 54–62. [CrossRef]

8. Stachowiak, G.W.; Batchelor, A.W. Physical Properties of Lubricants. In Engineering Tribology; Butterworth-Heinemann: Burlington, UK, 2006.

9. Liu, G.; Zhu, J.; Yu, L.; Wang, Q.J. Elasto-plastic contact of rough surfaces. Tribol. Trans. 2001, 44, 437–443. [CrossRef]

10. Etsion, I. Elastic–plastic contact analysis of a sphere and a rigid flat. ASME J. Appl. Mech. 2002, 69, 657–662. 11. Jackson, R.L.; Duvvuru, R.S.; Meghani, H.; Mahajan, M. An analysis of elasto-plastic sliding spherical

asperity interaction. Wear 2007, 262, 210–219. [CrossRef]

12. Mulvihill, D.M.; Kartal, M.E.; Nowell, D.; Hills, D.A. An elastic–plastic asperity interaction model for sliding friction. Tribol. Int. 2011, 44, 1679–1694. [CrossRef]

13. Cheng, S.; Luan, B.; Robbins, M.O. Contact and friction of nanoasperities: Effects of adsorbed monolayers. Phys. Rev. E 2010, 81, 016102. [CrossRef] [PubMed]

14. Yang, C.; Persson, B.N.J. Molecular dynamics study of contact mechanics: Contact area and interfacial separation from small to full contact. Phys. Rev. Lett. 2008, 100, 024303. [CrossRef] [PubMed]

15. Anciaux, G.; Molinari, J.F. Contact mechanics at the nanoscale, a 3D multiscale approach. Int. J. Numer. Methods Eng. 2009, 79, 1041–1067. [CrossRef]

16. Dubujet, P.; Ghaoudi, A.; Chaze, M.; Sidoroff, F. Particulate and Granular Simulation of Third Body Behavior. In Proceedings of the 22nd Leeds-Lyon Symposium on Tribology, Lyon, France, 5–8 September 1995; Elsevier: Amsterdam, The Netherlands, 1995; Volume 31, pp. 355–365.

17. Iordanoff, I.; Seve, B.; Berthier, Y. Solid third body analysis using a discrete approach: Influence of adhesion and particle size on macroscopic properties. J. Tribol. 2002, 124, 530–538. [CrossRef]

18. Tan, Y.; Yang, D.; Sheng, Y. Discrete element method (DEM) modeling of fracture and damage in the machining process of polycrystalline SiC. J. Eur. Ceram. Soc. 2009, 29, 1029–1037. [CrossRef]

19. Fillot, N.; Iordanoff, I.; Berthier, Y. Modelling third body flows with a discrete element method—A tool for understanding wear with adhesive particles. Tribol. Int. 2007, 40, 973–981. [CrossRef]

20. Müller, M.; Ostermeyer, G.P.; Bubser, F. A contribution to the modeling of tribological processes under starved lubrication. Tribol. Int. 2013, 64, 135–147. [CrossRef]

21. Müller, M.; Ostermeyer, G.P. A Cellular Automaton model to describe the three-dimensional friction and wear mechanism of brake systems. Wear 2007, 263, 1175–1188. [CrossRef]

22. Dmitriev, A.I.; Österle, W.; Wetzel, B.; Zhang, G. Mesoscale modeling of the mechanical and tribological behavior of a polymer matrix composite based on epoxy and 6 vol.% silica nanoparticles. Comput. Mater. Sci. 2015, 110, 204–214. [CrossRef]

23. Renouf, M.; Massi, F.; Fillot, N.; Saulot, A. Numerical tribology of a dry contact. Tribol. Int. 2011, 44, 834–844. [CrossRef]

24. Ju, Y.; Farris, T.N. Spectral analysis of two-dimensional contact problems. J. Tribol. 1996, 118, 320–328. [CrossRef]

(18)

Lubricants 2018, 6, 96 18 of 18

25. Polonsky, I.A.; Keer, L.M. A numerical method for solving rough contact problems based on the multi-level multi-summation and conjugate gradient techniques. Wear 1999, 231, 206–219. [CrossRef]

26. Liu, S.; Wang, Q.; Liu, G. A versatile method of discrete convolution and FFT (DC-FFT) for contact analyses. Wear 2000, 243, 101–111. [CrossRef]

27. Jacq, C.; Nelias, D.; Lormand, G.; Girodin, D. Development of a three-dimensional semi-analytical elastic-plastic contact code. J. Tribol. 2002, 124, 653–667. [CrossRef]

28. Chiu, Y.P. On the stress field due to initial strains in a cuboid surrounded by an infinite elastic space. ASME J. Appl. Mech. 1977, 44, 587–590. [CrossRef]

29. Chiu, Y.P. On the stress field and surface deformation in a half space with a cuboidal zone in which initial strains are uniform. ASME J. Appl. Mech. 1978, 45, 302–306. [CrossRef]

30. Wang, Z.; Yu, H.; Wang, Q. Analytical solutions for elastic fields caused by eigenstrains in two joined and perfectly bonded half-spaces and related problems. Int. J. Plast. 2016, 76, 1–28. [CrossRef]

31. Boucly, V.; Nelias, D.; Green, I. Modeling of the rolling and sliding contact between two asperities. J. Tribol. 2007, 129, 235–245. [CrossRef]

32. Simo, J.C.; Taylor, R.L. Consistent tangent operators for rate-independent elastoplasticity. Comput. Methods Appl. Mech. Eng. 1985, 48, 101–118. [CrossRef]

33. Green, A.P. Friction between unlubricated metals: A theoretical analysis of the junction model. Proc. R. Soc. Lond. A 1955, 228, 191–204. [CrossRef]

34. King, R.B.; O’sullivan, T.C. Sliding contact stresses in a two-dimensional layered elastic half-space. Int. J. Solids Struct. 1987, 23, 581–597. [CrossRef]

35. Aktary, M.; McDermott, M.T.; McAlpine, G.A. Morphology and nanomechanical properties of ZDDP antiwear films as a function of tribological contact time. Tribol. Lett. 2002, 12, 155–162. [CrossRef]

36. Van der Heide, E.; Schipper, D.J. On the frictional heating in single summit contacts: Towards failure at asperity level in lubricated systems. J. Tribol. 2004, 126, 275–280. [CrossRef]

37. Bosman, R. Mild Wear Modeling in the Boundary Lubrication Regime. Ph.D. Thesis, University of Twente, Enschede, The Netherlands, 2011.

38. Shakhvorostov, D.; Pöhlmann, K.; Scherge, M. Structure and mechanical properties of tribologically induced nanolayers. Wear 2006, 260, 433–437. [CrossRef]

39. Yang, L.; Wang, D.; Guo, Y.; Liu, S. Tribological behaviors of quartz sand particles for hydraulic fracturing. Tribol. Int. 2016, 102, 485–496. [CrossRef]

40. Sandeep, C.S.; Senetakis, K. Grain-scale mechanics of quartz sand under normal and tangential loading. Tribol. Int. 2018, 117, 261–271. [CrossRef]

© 2018 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).