A Boundary Lubrication Model and Experimental Study Considering ZDDP Tribofilms on Reciprocating Friction Pairs

Boundary lubrication state may dominate the friction pairs operating under severe conditions, yet its mechanism is not clearly understood and related numerical models are still lacking. A boundary lubrication model considering zinc dialkyldithiophosphate (ZDDP) tribofilms, which impact friction and wear performances, was developed in this study. A series of reciprocating experiments were conducted to verify this model and also to investigate the effects of the tribofilm on friction and wear under various temperatures and loads. Moreover, the experimental data were employed to modify the tribofilm removal equation, which enabled the present boundary lubrication model to be applied under a wide range of loads. The results showed that the friction force and wear depth both decline with the increasing lubricant temperature due to thicker tribofilms formed. As the load becomes heavier, the wear depth keeps increasing, while the tribofilm thickness first increases then decreases.


Introduction
Nowadays the structures of friction pairs are designed to be more and more compact and their operating conditions are becoming more and more severe to improve power efficiency. As a result, the friction pairs are inevitable to experience or even maintain boundary lubrication state under the harsh operating conditions [1]. In this situation, the contacting surfaces are not well lubricated by a fluid film [2] but separated by a protective surface-bonded tribofilm in the presence of lubricant additives. The tribofilm will impact the states of fluid lubrication, asperity contact, friction, and wear during surface rubbing [3]. However, a numerical model to describe these effects is still lacking, especially for engineering friction pairs under dynamic working conditions. Therefore, a boundary lubrication model considering the tribofilms, which is calibrated and verified by experiments, is developed in this study.
One of the most efficient tribofilms to reduce surface wear are ZDDP tribofilms and many studies have been conducted on its formation mechanism, friction and wear performances. ZDDP tribofilms tightly bonds to the rubbing surfaces with a thickness as high as 200 nm [4]. They consist of rough, patchy, and pad-like features that are composed of pyro-or orthophosphate glasses in the bulk, with an outer nanoscale layer of zinc polyphosphates and a sulfur-rich layer near the substrate surface [3]. It was demonstrated that a high contact pressure and lubricant temperature could promote the formation of ZDDP tribofilms [5], but it would be inhibited if the pressure was too high [6]. Gosvami et al. revealed that the tribofilm growth was stress-dependent by singleasperity sliding contact experiments and its growth rate fitted a stress-activated Arrhenius model [7]. They assumed this driving stress was normal pressure [7]. While Zhang and Spikes argued that it was primarily the shear stress [4], and they showed that tribofilms could grow even under full film lubrication where there was no asperity contact if only the shear stress was high enough [8]. Zhang et al. also studied the boundary friction of ZDDP tribofilms under both rolling/sliding and pure sliding contacts [9]. They found that the coefficient of friction (CoF) of the tribofilms depended on the alkyl structure of ZDDPs. It is incontrovertible that ZDDP tribofilms can efficiently mitigate surface wear especially at extreme pressure [10], yet studies also found that oils with higher ZDDP concentrations produced more micro pitting [11,12]. The stress-activated model for the tribofilm growth is widely accepted. However, the process of tribofilm removal is still under debating and there is no reliable equation that can model it.
Most researches on ZDDP tribofilms have been conducted experimentally, while a practical boundary lubrication model for engineering friction pairs has not yet been developed. Nevertheless, there are still several studies striving to model the tribofilm formation and its effects using numerical methods. Ghanbarzadeh et al. proposed a semi-deterministic wear model considering the effects of ZDDP tribofilms [13]. In their model, contact pressure was calculated by the elastic-perfectly plastic theory and surface wear was evaluated by a modified Archard's wear model [14] whose wear coefficient was varied with tribofilm growth. Later Akchurin et al. also developed a model for the simulation of tribofilm formation by integrating a boundary element method-based contact model and the stress-activated Arrhenius tribofilm growth equation [15]. Surface wear in the presence of tribofilms was regarded as mild wear resulting from atom loss of substrate material by tribofilm removal. Considering both fluid lubrication and tribochemistry reaction, Azam proposed a deterministic mixed lubrication model where the tribochemistry, lubrication, and asperity contact were coupled [16]. Similarly, a modified Archard's wear model was applied in their work and asperity contact pressure was used to estimate surface wear. Recently, Chen et al. modeled the formation and removal of ZDDP tribofilms on rough surfaces, and applied this model to the cylinder liner surface of an internal combustion engine (ICE) [17]. They simulated the tribofilm growth on both top dead center (TDC) and bottom dead center (BDC), which induced a distinguished wear topography at the end of the simulation. Yet these models/simulations are all based on deterministic asperity contacts, which are quite time-consuming when coupled with hydrodynamic lubrication, as the formation of steady tribofilms may require several-hours sliding. Besides friction force in the presence of tribofilms is not predicted in these simulations.
Therefore, a boundary lubrication model, in which ZDDP tribofilms growth and removal were considered, was presented in this study to capture their effects on both friction and wear performances. In this model, asperity contact pressure was obtained using a statistical contact model, which is more practical and efficient for engineering friction pairs. Then experiments for reciprocating friction pairs under various conditions were conducted to calibrate and verify this model. Moreover, a modified tribofilm removal equation with an impact factor reflecting the state of asperity contact was proposed based on the experimental results. This study could help to understand the effects of tribofilm formation on the friction and wear under boundary lubrication.

Boundary Lubrication Model
The research on boundary lubrication has commenced for a complete century ever since it was first systematically studied by Hardy [18] in 1922. It is believed that from a macroscopic perspective, the rubbing surfaces of friction pairs working under boundary lubrication are separated by thin films, which are called boundary films. Boundary films can be generally classified into absorption films and reaction films [19]. Adsorption films by physisorption or chemisorption are not considered in the present study because they are easily desorbed under loaded conditions [1]. As for the main topic in the present study, ZDDP tribofilms growth on the rubbing surfaces by tribochemistry reactions between lubricant additives and substrate materials belong to the latter. ZDDP tribofilms partly cover substrate surface with an evolving thickness and they directly impact the contact state during rubbing. In the boundary lubrication model, the tribofilm as well as the substrate contact pressure are obtained by a layered statistical asperity contact model. The yielded results are then used to evaluate the boundary friction force and surface wear. At the same time, the tribofilm thickness is determined by the stress-activated growth model and the removal equation.

Layered Statistical Asperity Contact Model
The tribofilm formed by anti-wear lubricant additives firmly bonds on the substrate surface. As the tribofilm grows, it will bear more and more applied load. Therefore, the evolution of the tribofilm thickness determines the asperity contact state. A layered contact model is proposed in this paper, where there are two layers, the tribofilm and the substrate from the outside to the inside surface. The contact of two rough surfaces can be equivalent to one equivalent single rough surface in contact with a rigid flat surface, and it shows as: The tribofilm thickness is represented by h tri in Fig. 1. The reference plane is located at the mean height of the substrate surface asperities. Considering the rough, patchy, and pad-like features of ZDDP tribofilms, the substrate is partly and unevenly covered by the tribofilm. Thus the substrate-to-substrate contact and tribofilm-to-tribofilm contact exist simultaneously. Analogically, in the layered contact model, the rigid surface will contact with the tribofilm and also penetrate it to contact with the substrate. Assume the separation between the rigid flat surface and the reference plane is h and the probability density function of the asperity height is Φ(h) , the substrate contact pressure ( p sub ) according to the statistical model proposed by Greenwood and Tripp (GT) [20] is given as: where E ′ sub is the equivalent elastic modulus of the substrate material. σ is the composite surface roughness, is the asperity density, and is the asperity radius. These parameters for a given friction pair are shown in Sect. 3.1.
The tribofilm contact pressure ( p tri ) is closely related to its thickness and is given as: The integration from h − h tri to h represents the contact proportion of the tribofilm and it is actually the difference by subtracting the contact proportion of the substrate (from h to ∞ ) from the whole contact area (from h − h tri to ∞ ). E ′ tri is the equivalent elastic modulus of the tribofilm, which ranges from 40 to 90 GPa [21] and will be discussed in Sect. 2.3.3.
The total asperity contact pressure is the sum of the substrate and the tribofilm contact pressure: The contact model above is simple yet logical and it depicts the process that the tribofilm will bear more applied load as it grows thicker. The model will be validated by comparing the simulation and experimental results in Sect. 4.
The hydrodynamic pressure ( p lub ) generated between rubbing surfaces is obtained by the average Reynolds equation [22]. The lubricant viscosity has a significant effect on the performance of the hydrodynamic pressure and varies with the pressure and temperature. It is characterized by Roelands's model [23]. Moreover, the load is borne by both the fluid and asperities. The distributions of hydrodynamic pressure and asperity contact pressure are supposed to satisfy the load balance equation. It can be expressed as: where w is the applied load and A is the contact area.

Boundary Friction and Wear Evaluation
Bowden and Tabor in the 1940s [24,25] proposed a widely accepted model for boundary lubrication. In their model, the friction force of the rubbing surface comprises of adhesive shear at metallic contact areas and viscous shear at oil lubricating areas. This model in understanding the process of boundary lubrication is simple yet efficient. In the present model, the total friction force F fri comprises of the tribofilm and substrate contact friction forces F tri and F sub as well as the lubricant shear force F lub : The tribofilm and substrate contact friction forces F tri , F sub are obtained as: where tri and sub are the CoFs of the tribofilm and substrate respectively. The CoF of ZDDP tribofilms varies slightly within the temperature range of 60-100 °C [26] and it is 0.1 used in [15]. In the present study, it is obtained by experimental calibration. The CoF of the substrate is obtained by dry sliding tribotest for the given friction pair.
For lubricant shear force F lub , it calculates as: where the terms f , fs and fp are friction-induced flow factors [22,27]. u is the sliding speed, is the lubricant viscosity, and x is the coordinate parallel to the sliding direction.
In the boundary lubrication, the substrate is patchily covered by the tribofilm, and substrate-to-substrate contact and tribofilm-to-tribofilm contact both exist. Therefore, the loss of the substrate origins from two parts, the adhesive wear by substrate-to-substrate contact and the mild wear by tribofilm removal [28]. The adhesive wear ( h mech w ) is derived from the direct substrate contact. It can be evaluated by Archard's wear model [14]: While the latter form of wear is caused by diffusion of the substrate atoms into the tribofilm, which are subsequently loss if the tribofilm is removed. The loss of the substrate material due to atom diffusion is calculated considering both the tribofilm removal rate and the concentration of the substrate material atoms in the tribofilm. The concentration in the tribofilm is higher close to the substrate-tribofilm interface and reduces away from the interface. The loss of the substrate material (mild wear, h mild w ) can then be calculated [15]: where the term e −ch tri is the concentration of substrate material and c is an unknown constant which needs be determined experimentally. In this paper, it is obtained from [15]. ( is the tribofilm removal rate and is determined in the following section.

Tribochemistry Reaction
The tribofilm thickness is an essential variable in the layered contact model (Sect. 2.1) to determine the contact portion of the tribofilm. The tribofilm is continuously worn and replenished during surface rubbing and the evolution of the tribofilm thickness stems from the balance of its growth and removal rates. The change rate ( h tri t ) in the tribofilm thickness is calculated by subtracting the removal rate from the growth rate: where ( are the tribofilm growth and removal rate respectively. Equation 11 also denotes the net growth rate of the tribofilm.

Tribofilm Growth
The growth rate ( ( h tri t ) gro ) of the ZDDP tribofilm is calculated based on Eq. 12, which is known as the stress-activated Arrhenius equation [7].
In this equation, Γ 0 is the pre-factor, k B is Boltzmann's constant, T is the absolute temperature, ΔU act is the internal activation energy (i.e., the energy barrier in the absence of stress), and ΔV act is the activation volume. It is worth noting that in this equation is shear stress [4] and is the CoF. A formulated oil (10W40) was used in the presents study, which is a mixture of base oil with other additives, including detergents, dispersants and friction modifiers [29,30]. These additives will compete with ZDDP to adsorb on the sliding surfaces [31], slowing down the formation of ZDDP tribofilm. Therefore, the parameters in Eq. 12 ( ΔU act and ΔV act ), which governs the tribofilm growth, need to be calibrated by fitting the experimental results.
The stress-activated Arrhenius equation was initially derived through a single asperity (an AFM tip) sliding contact experiment and subsequently applied in deterministic models in which the contact pressure of identical asperity was calculated [15,17]. In this study, a statistical model is used and in each mesh grid (10 × 10 μm 2 ), the asperity contact pressure is averaged based on the statistical characteristics of the surface topography. The tribofilm thickness in each mesh grid is represented by a uniform value based on the contact pressure and the tribofilm formation is indeed determined by pressure/ shear stress. It yields different tribofilm thickness under various contact states, which may result from different rubbing time, lubricant temperature or applied load. This is practical and more efficient as the size of an engineering friction pair (like piston ring and cylinder liner pair in ICEs) is in macroscopic scale and the computation cost of a deterministic model is quite high.

Tribofilm Removal
It was observed in experiments that the removal rate would increase linearly with the growth of the tribofilm [32], which is attributed to the difference in the wear resistance between the bulk of the substrate (higher wear resistance) and the free surface (lower wear resistance). Based on Archard's wear model with a variable wear coefficient, the tribofilm removal rate ( ( ) is calculated as follows: where K is the removal rate (or wear rate) coefficient, K sub is the wear rate of the substrate, and is a constant describing the durability change rate inside the tribofilm. They can be calibrated by experiments. Equations 14 and 15 indicate that the tribofilm removal rate varies linearly with its thickness, and similar models were also used in [17,24].
However, experiments also showed that the tribofilm would collapse under extremely high contact pressure and a threshold pressure for the tribofilm to maintain a stable thickness exists [7,33]. Above the threshold pressure, the tribofilm formation is not stable and the removal process would be dominant over the growth process. In this situation, the results using the linear tribofilm removal rate equation (Eqs. 14 and 15) may be far from the realistic. Based on Eqs. 14 and 15, the authors propose a modified tribofilm removal equation and introduce an impact factor reflecting the state of asperity contact based on the experimental results. The details will be discussed in Sect. 4.2.

Mechanical Properties of the Tribofilm
ZDDP tribofilms have different mechanical properties compared to the substrate, and several studies have focused on the assessment of these properties [29,34,35]. The mechanical properties of the tribofilm determine the contact pressure and then the growth rate. It is generally agreed that the hardness of the tribofilm varies with its height basically in a linear way. Therefore, as the tribofilm grows and removes, its hardness evolves along with its height. To model this variable hardness, the approach presented by Andersson et al. is used in the present study [36].
Regarding as the effect of temperature, it was found that the hardness and Young's modulus of ZDDP tribofilms closely depended on the temperature, especially under a high temperature [21,37] (above 220 °C, ZDDP tribofilms start to degrade [38]). It was observed by experiments that the hardness of ZDDP tribofilms would be lower at a higher temperature [37], and based on this, Bosman and Schipper [39] proposed a compensation factor for the hardness as a function of temperature. According to Ref. [15], the hardness of the tribofilm at a given temperature with a certain thickness can then be obtained by multiplying a compensation factor from [39]. The mechanical model is generally reasonable as it follows the trend of the tribofilm that its hardness decays with film thickness and temperature. The Young's modulus of the tribofilm is related to its hardness, and after a threshold of hardness, it shows a linear increase with a slope of 35 [39,40].
Asperity interaction and fluid viscous shear could produce considerable heat and the temperature in the contact area will rise to some extent under boundary lubrication. Studies showed that flash temperature might play a role in promoting the growth of ZDDP tribofilms at high sliding speed, but not at low sliding speed [5,41]. In the present study, the sliding speed is extremely low (an average speed of 0.04 m/s in a stroke), and thus the temperature variation due to frictional heating is minor. Nevertheless, the temperature increase was evaluated using the method by Bosman and De Rooij [42] and the maximum was found to be less than 1 °C. As a result, the flash temperature at the sliding interface has a marginal effect on the formation of ZDDP tribofilms as well as the macroscopic friction and wear. Thus the temperature increase by friction is neglected in the current simulations.

Computation Flowchart
The flowchart of the boundary lubrication model is shown in Fig. 2. The computation starts with inputting the initial parameters, including the mechanical properties, macroscale geometry, roughness, load, temperature and sliding speed. Then the oil film distribution is determined by the ring profile and an assumed minimum film thickness. The oil film pressure and asperity contact pressure could be obtained by solving the average Reynolds equation and the layered statistical contact model respectively. Initially, there is no tribofilm ( h tri = 0 ). After the calculation of contact pressure and shear stress, the tribofilm growth and removal rate are determined using Eqs. 12-15. The tribofilm thickness is then updated accordingly. Successively, it is used in the contact model as the thickness of the tribofilm layer. The mechanical properties are updated based on the new tribofilm thickness. Next, the convergence on load balance is checked using Eq. 4. If the convergence criterion is not satisfied, adjust the minimum film thickness and the calculation repeats. Otherwise, the friction force and wear depth are obtained from Eqs. 5-10. The simulation time is then increased by a certain time step and the sliding speed also updates. There are 180 time steps in a stroke, and therefore, a time step is 0.00139 s . The whole computation will stop when reaches the set time (2 h).  Fig. 2 Flowchart of the computation procedures for the boundary lubrication model

Boundary Lubrication Experiment
The parameters in the boundary lubrication model need to be calibrated based on the experimental results. Successively the tribofilm thickness, CoF, and wear depth predicted by the model were compared with the ones measured in the reciprocating experiments under various lubricant temperatures and applied loads.

Test Materials
The piston ring and cylinder liner segments were selected to be the test samples. They simulated a piston ring and cylinder liner (PRCL) contact, which usually experiences boundary lubrication state near the TDC and BTC in ICEs [43]. In this study, the piston ring segments were cut off from a top compression ring with a diameter of 110 cm. The ring working surface had diamond-like carbon (DLC) coating. Spheroidal graphite cast cylinder liner in an ICE, which was mating with the ring, was sectioned into liner samples every 1.1 • in the circumferential direction (corresponding to about 10 mm arc) and 15 mm in the axial direction. The thickness of the liner sample is 5 mm Note that the liner sample has one curved surface rubbing against the piston ring. Each pair formed a conformal contact and the contact geometry was cylindrical with an area of 1.2 × 1.1mm 2 . Surface roughness measurements were performed by Mitutoyo CV-4500S4 surface profiler over a data length of 4 mm and lc filter cut off length of 0.8 mm.
The mechanical properties and surface roughnesses of the piston ring and liner samples used in this study are listed in Table 1.
To simulate the lubrication environment in engineering friction pairs more realistically, an industrial lubricant, SAE 10W40, which is a ZDDP-fortified mineral engine oil for classic cars, was used in the present study. The properties of the lubricant are shown in Table 2.

Experimental Rig
The experiments were conducted using a reciprocating tribometer (MFT-R4000). The tribometer is shown in Fig. 3a and b, and a schematic diagram of the ring-on-liner contact is shown in Fig. 3c. The ring segment was oscillating against the liner segment with a stroke length of 10 mm while the latter was fixed in the oil bath. A conformal ring holder was used to prevent any inclination to the counterbody. The reciprocating frequency was as low as 2 Hz to ensure the boundary lubrication appears, which resulting in an average speed of 0.04 m/s in a stroke. Note that the sliding speed was sinusoidal in a stroke and the oil film thickness was calculated based on the sliding speed at a certain position, as described in Sect. 2.4. Before each test, the oil bath was filled with lubricant, which completely covered the piston ring, i.e. the piston ring was under fully flooded lubrication. The lubricant temperature was set by the heating system under the oil bath from ambient temperature to 200 °C. The maximum applied load was up to 300 N. In this study, the duration time of each test was 2 h. The CoF, the test liner temperature (measured by a thermocouple at the bottom of the liner segment), the applied load, and the stroke were continuously monitored and recorded by the tribometer software.
Before each test starts, 30-min running-in was conducted under the load of 20 N and the temperature of 30 °C. Continually, the load and temperature were increased to the set values. Each test was repeated three times to ensure consistency of results. Considering the surface roughness changes rapidly during running-in and slowly thereafter [44], the surface roughnesses of the test ring and liner samples after running-in were measured and used as input values for the simulation, which were 0.408 and 0.571 µm, respectively. Note that the piston ring is resilient and it is compressed and then assembled into an engine, which means the ring segment curvature could be different to the liner without compression. Yet the running-in process could make the two surfaces of the ring and liner segments conformal.

Surface Characterization and Tribofilm Measurement
The wear depth of the liner segment was measured using the surface profiler. Measurements were taken from multiple locations across the wear track to determine standard deviation and the final wear depth was the average of these measurements. The surfaces after the experiments were observed without any conductive coating by the scanning electron microscope (SEM) operating under high vacuum. The chemical composition analyses were carried out by energy dispersive X-ray spectroscopy (EDX) and quantitative analysis of the EDX spectra was performed using the Aztec software.
In previous studies, the tribofilm thickness could be mapped in situ by spacer layer interferometric method cooperated with mini traction machine (MTM) ball-on-disc rig [30,45] or scanned ex situ by atomic force microscopy (AFM) [46,47]. The surface roughness of the tested samples in these studies was as low as dozens of nanometers or even lower. Yet in the present study, the surface roughness of the tested liner segment was much higher (hundreds of nanometers). It would be very difficult to measure the tribofilm thickness using these two methods aforementioned.
The resolution in the height direction of a 3D laser-scanning confocal microscope (Keyence VK-X100/X200) is sub nanometer and its measurement area can be expanded to millimeters. Therefore, it is suitable for surface topology measurement and tribofilm thickness evaluation in the present study. It has been reported that ethylenediaminetetraacetic acid (EDTA) could dissolve ZDDP tribofilms without damaging the substrate material [48]. After each test, the liner segment was carefully cleaned using anhydrous ethanol. A texture dimple was then marked within the macroscopic wear scar on the liner surface by a highly accurate picosecond laser machine (Trumpf, TruMicro 5050, German). The depth of the texture dimple was micrometer scale and it completely penetrated the tribofilm reaching the substrate material. Subsequently, the same region including the texture dimple was measured using the confocal microscope twice: (i) before removing the ZDDP tribofilm; (ii) after removing the ZDDP tribofilm with 0.05 mol∕L EDTA disodium salt in The tribofilm distribution and thickness were then obtained by subtracting the latter measurement results from the former (see Fig. 4). Similar tribofilm measurement method was used in [46,49].

Results and Discussion
In the boundary lubrication model, there are some key parameters that require calibration before predictions of surface wear and friction force can be made. They were determined using experimental results for a given operating condition and then used for all other simulations under different conditions. The calibration experiment was conducted under the lubricant temperature of 90 °C and load of 160 N . The CoF of the substrate was obtained by dry sliding tribotest and here it was 0.15 for the present test samples. As for the CoF of the ZDDP tribofilm, it was determined by fitting the simulated CoF with the recorded CoF by the calibration experiment, and it was 0.11. In Eq. 12, ΔU act and ΔV act have a significant impact on the growth rate of the tribofilm and they were calibrated by fitting the experimental results. These two parameters were 0.711 eV and 89 Å 3 , respectively. In the tribofilm removal and substrate wear model, K sub (the wear rate of the substrate) and were calibrated based on the wear depth and tribofilm thickness measurements. In the present study, equaled to 0.21 nm −1 and K sub was 1.43 × 10 −16 Pa −1 .
The simulation results below were obtained using the calibrated boundary lubrication model. The sliding speed is sinusoidal in a stroke (Fig. 5a, red line), leading to a varying oil film thickness (Fig. 5a, blue line), which makes the friction force also varying (Fig. 5a, green line)-highest at the two ends of a stroke and lowest in the middle of a stroke. The oil film thickness ratio, which equals to h∕ , represents the lubrication state and was also predicted in one stroke (blue dash line) under the lubricant temperature of 100 °C and load of 150 N . Its maximum value is 0.18 and is much lower than the threshold of boundary lubrications. Figure 5b The CoFs recorded by the tribometer and predicted by the boundary lubrication model exhibit a consistent tendency with the sliding time: they both decline rapidly initially and then remain stable. This is different from the results of previous ball-on-disk experiments, like [50,51]. It is explained as follows. The surface roughnesses of the test specimens (ball and disk) were only tens of nanometers, which were lower than the tribofilm thickness generated on the specimen surface during rubbing. It could substantially reshape the topography of the initial rubbing surface and increase the roughness significantly. Successively, the oil fluid entrainment was inhibited and asperity contact friction increased [51]. However, in the present study, the surface roughnesses of the test ring and liner samples were as high as hundreds of nanometers, which were one order of magnitude higher than the tribofilm thickness (tens of nanometers). The formation of ZDDP tribofilms had a minor effect on the surface roughness. On the contrary, the CoF of the tribofilm (0.11) Fig. 4 The schematic principle of tribofilm measurement was lower than the one of the substrate material (0.15). Thus the total friction would decrease (Fig. 5b, red and blue dash lines) as the tribofilm grows (Fig. 5b, green dash line). At the same time, the applied load would be more born by the tribofilm with a higher thickness according to the layered statistical contact model (Sect. 2.1). The part born by the oil film became lighter, which resulted in the oil film thickness being lifted slightly (Fig. 5b, orange dash line). Azam et al. [52] simulated the effects of tribofilms on lubrication. They found that the presence of tribofilm increased the lubricant film thickness and also improved the lubrication regime by helping entrain more lubricant inside the contact area.
The results from SEM imaging and EDX spectra are shown in Fig. 6. The SEM image clearly shows wear tracks parallel to the sliding direction (S.D.) on the surface of the liner sample. As typical elements of ZDDP molecular, the enrichment of zinc and phosphorus (in EDX chemical maps) on the wear track indicates the formation of corresponding poly-phosphate tribofilm. Note that higher concentrations of zinc and phosphorus appearing in the wear tracks means a thicker tribofilm, which was also observed in previous research using Raman spectroscopy [46]. It demonstrates that the tribofilm mainly forms in the wear track area, where the contact pressure and shear stress are much higher. According to the stress-activated tribochemistry theory, the tribofilm growth will be promoted in these areas. Similar element distributions where zinc and phosphorus are concentrated in the wear tracks were obtained in the ball-ondisk reciprocating experiments [53]. High magnification SEM image and EDX spectra (the bottom panels of Fig. 6) confirm that the enrichment of zinc and phosphorus emerges in the wear track.
As mentioned in Sect. 3.3, the tribofilm measurement was using the 3D laser-scanning confocal microscope and calibrated by a texture dimple. The dimple position was on the middle of the stroke, i.e., the center of the macroscopic wear track on the liner sample for each test, considering there might be wear debris or other contaminations accumulating at the ends of the stroke. The 3D surface topographies pre-EDTA and post-EDTA scanned by the confocal microscope were set on the same baseline under the same magnification scale. Then their gap in the height direction could be obtained by subtracting the latter measurement results from the former, representing the average ZDDP tribofilm thickness (Fig. 4). The baseline was defined by the bottom of the texture dimple with a diameter of 150μm marked within the wear scar. To avoid the complex metal oxidation and sublimation area by lase ablation, the tribofilm thickness evaluation areas ( 150 × 150μm 2 ) were chosen away from the texture dimple in its four directions (seen Fig. 7a). Figure 7b displays the cross-sectional profile of the texture dimple by the surface profiler. On the one hand, the depth of the texture dimple is about 1.8μm , and it completely penetrates the tribofilm layer reaching the substrate material. Therefore the EDTA treatment doesn't change the topography of the dimple bottom. On the other hand, the bottom profile is relative smooth and quite flat owing to the high-accuracy picosecond laser machine. It guarantees the reliability of the dimple bottom to be the baseline. The measurement method described above is practical as the ZDDP tribofilm is uneven Fig. 6 Representative SEM images and EDX maps of P and Zn demonstrating a tribofilm generated on the liner segment from the reciprocating tribometer and not continuous. The measurement results of the tribofilm thickness under various temperatures and loads will be shown and discussed below. Note that the CoFs presented in Figs. 10 and 13 under various lubricant temperatures and applied loads are the final/steady ones at the end of the test, considering relatively steady CoF and tribofilm thickness are formed after 2-h sliding.

Effects of Lubricant Temperature
The tests were firstly conducted over a wide range of lubricant temperatures under the load of 150 N . Figure 8 summarizes the tribofilm measurement areas (marked by cyan dash squares in Fig. 7) before and after EDTA treatment under the lubricant temperatures from 40 to 140 °C. These focal microscope images clearly demonstrate that ZDDP tribofilms have been entirely removed from the surface by EDTA treatment and the substrate material is exposed. The quantitative tribofilm thickness values are plotted as a function of temperature in Fig. 9. As the lubricant temperature rises, the tribofilm thickness increases rapidly and basically in an exponential way. The predicted tribofilm thickness by the boundary lubrication model agrees well with the experimental data, which again confirms the validity of the stressaugmented thermal activation model. This trend is similar to the previous works observed in ball-on-disk experiments [4,8] and simulated by a deterministic tribofilm growth model [15]. It is noticeable that the topographies of the tribofilm shown in Fig. 9 display a strip-like distribution parallel to the sliding direction. It is reasonable as the tribofilm mainly grows in the wear track area and it also matches the EDX map shown in Fig. 6. The tribofilm thickness results demonstrate that the temperature can evidently promote the tribochemistry reaction forming ZDDP tribofilms.
To further study the effect of lubricant temperature on the friction and wear under boundary lubrication, the CoFs and wear depths measured in experiments and predicted by simulations with and without tribofilm are compared in Fig. 10. Here simulations without tribofilm means that the stress-augmented thermal activation model is not included and the tribofilm formation is not considered in the simulations. Therefore, there are only lubricant fluid pressure and substrate contact pressure to support the applied load and the surface wear merely stems from mechanical interaction. It is in fact a mixed lubrication model. As shown in Fig. 10 (green line chart), the CoFs predicted by the model without tribofilm tends to rise with the increasing temperature. This is easy to understand as the lubricant viscosity will decline when the temperature rises, which leads to a thinner oil film. Then more asperities come into contact. As a result, the friction force increases. However, the CoFs recorded by the tribometer (red line chart) in the same temperature range shows an opposite tendency: it declines with the increasing temperature. When the growth of tribofilm is included in the lubrication model, the simulated CoFs (blue line chart) matches the experimental data well. It is reasonable as the tribofilm thickness increases with the rising temperature ( Fig. 9), and according to the layered statistical contact model (Sect. 2.1), the tribofilm will share more applied load. Due to its low friction characteristics, the total friction force calculated by Eq. 5 will decrease. In the ball-on-disk experiments by Xu et al. [46], the tribofilm formation rate was also promoted under a higher lubricant temperature followed by a decreased CoF. The inverse trends presented by the simulations with and without tribofilm emphasize the difference between the boundary and mixed lubrication.
Analogously, the wear depths (Fig. 10, histogram charts) predicted by the model excluded tribofilm becomes larger under a higher lubricant temperature owing to the severer contact state. When the tribofilm growth model is integrated in the simulations, the wear depths decrease with the rising temperature, which is also consistent with the experiment data. As aforementioned in Sect. 2.5, here substrate wear in the present model comprises two parts: the mild wear and the adhesive wear. The mild wear Fig. 7 The texture mark in the wear scar on the liner surface: a the optical microscope image of the texture dimple and tribofilm measurement areas; b the texture cross-sectional profile originates from the diffusion of the substrate material atoms into the tribofilm and subsequently lost by tribofilm removal [39]. The atom concentration of substrate material substantially determines this kind of wear (Eq. 10) and it is reducing with the increasing tribofilm thickness.
The mild wear will therefore decrease with the increasing temperature as a thicker tribofilm is generated. In the simulations by Akchurin et al. [15] as well as by Ghanbarzadeh et al. [13], the predicted mild wear was also decreasing under a higher lubricant temperature. For all that, the authors would like to point out that in the boundary lubrication experiments [13], only mild wear was considered, yet substrate contact and adhesive wear would inevitably occur noting the tribofilm is highly uneven and patchy. Apart from the mild wear, the adhesive wear also displays a negative relation to the temperature shown in Fig. 10. The reduction in the adhesive wear is attributed to the formation of the protective tribofilm film, which avoids the direct contact and effectively inhibits the adhesive wear of substrate materials. The wear resistance of ZDDP tribofilms is more significant under a higher lubricant temperature, which is also consistent with the simulation results by Chen et al. [17].
The simulation analysis above can properly explain the experimental results and it demonstrates that the model proposed in the present paper is able to capture the friction and wear behaviors of the boundary lubrication.

Effects of Applied Load
Contact pressure can enhance the growth of ZDDP tribofilms, while studies have also shown that the tribofilm wear/ removal would be dominant if the pressure is too high [7]. The final tribofilm thickness depends on the balance of the tribofilm growth and removal [54]. To study the effect of the contact pressure and shear stress, tests were then conducted over a wide range of applied loads under the lubricant temperature of 100 °C. The confocal microscope images of the tribofilm evaluation areas before and after EDTA treatment from 50 to 300 N are summarized in Fig. 11. The quantitative tribofilm thickness values are plotted as a function of load in Fig. 12 and similar strip-like distributions of the tribofilm topography are also displayed. The experimental data (red line chart) shows that as the applied load goes up, the tribofilm thickness will be increasing first but declines under the maximum load. It is supposed that a moderate applied load boosts the tribofilm growth but a much higher load could depress its formation by augmenting the tribofilm removal rate. In the single AFM tip sliding experiments by Gosvami et al. [7], they observed an unstable growth rate when the contact pressure exceeded 5.2 GPa . Bayat et al. [33] pointed out that there might be a threshold of pressure for the tribofilm to maintain a stable thickness. Above this threshold, the tribofilm formation is depressed and the removal is dominant. Yet the threshold pressure is sensitive to the surface roughness and thus hard to determine.
In the present study, the contact area was 1.2 × 1.1mm 2 , and thus the contact pressures varied from 37.88 to 227. 27 MPa , corresponding to the applied loads ranging from 50 to 300 N . Note that the contact pressures in the present study were generally lower than the ones in previous ball-on-disk experiments, in which point contacts formed [4,8]. This was because the contact between the ring and liner samples was conformal and the contact area was relatively large, compared with Hertz contact, i.e., point or line contact. Nevertheless, considering the surface roughness of the test ring and liner samples were one order of magnitude higher than the ones of the test ball and disk, there would be more local extreme-pressure/shear areas, where the tribofilm removal was dominant, resulting in thinner tribofilms. Besides, the authors would like to point out that under the highest applied load, the tribofilm growth rate using Eq. 12 would continually increase, and it is the tribofilm removal that increases even faster and inhibits a thicker tribofilm forming. The authors initially used a linear tribofilm removal equation (Eqs. 14 and 15) based on the Archard's  wear model whose removal rate coefficient varies with the tribofilm thickness. The predicted tribofilm thickness using this model agrees well with the experimental results under the moderate applied loads. However when the applied load exceeds 250 N , the gap between the simulation and the experiment is unbearable. Considering the double-edge effect of the contact pressure, the authors proposed a modified tribofilm removal equation by introducing an impact factor reflecting the state of asperity contact to the removal rate coefficient K . Specifically, the present model introduced Here β p p 0 + 1 is the multiplication factor and is highly nonlinear with the contact pressure. β and are constants and can be calibrated according to experiments. In the present study, they are 24.86 and 11.96, respectively. p 0 is the contact pressure when the nominal distance between the two rough surfaces is zero, that is, h = 0 . Note that the minimum values of β p p 0 + 1 equals to 1 when p is 0 and its maximum is β + 1 . Under a light load condition where the contact state is mild, the tribofilm removal rate yielded by the modified model is quite close to the one by the linear model. While the former will be significantly higher than the latter by dozens of times under a severe contact state. As show in Fig. 12, the predicted tribofilm thickness by the modified model is in line with the experimental data in the whole range of the applied load.
To further examine the modified tribofilm removal equation, the CoFs and wear depths predicted by the linear and modified model are both shown. They are compared with the measured results in the experiments and discussed below. The measured CoFs (Fig. 13, red line chart) declines first following the increasing applied load and reaches its minimum values about 0.106. Then it is increasing as the load becomes higher. The simulation results using the linear removal equation (blue line chart) keep decreasing along with the increasing load and are getting deviated from the experiment when the load exceeds 250 N, just like the tribofilm thickness results in Fig. 12. Using the modified model, the simulation results (orange line chart) can catch the variation of the experimental results. The authors propose that the rise of CoF under the high applied load may result from the decrease of tribofilm thickness, comparing the results of the linear and modified model. The wear depths no matter by simulations or experiments are increasing with the applied load. For the simulation results, the mild wear only counts a small portion and it is the adhesive wear that mainly contributes to the total wear increase. It can be concluded from Eq. 10 that although the removal rate is increased under a higher load, the concentration of substrate material is also lower as the tribofilm increases (Fig. 12). Therefore, the mild wear changes slightly. It is noteworthy that when the load increases from 250 N to 300 N , the experimental wear depth rises significantly, which is also matches the variation of CoF. The simulation using the modified model shows a similar tendency with the experiment result.
The results and analysis above demonstrate that the linear tribofilm removal equation may not suit for severe contact state, while the modified equation shows a favorable adaptation under a wide range of applied load.

Conclusions
A boundary lubrication model considering the ZDDP tribofilm growth and removal was developed, in which the asperity contact pressure is calculated using a layered statistical contact model. The simulated tribofilm thickness, friction and surface wear were compared with the reciprocating experiments under various temperatures and loads. The results can be summarized as follows:  1. The ZDDP tribofilm growth will be promoted under a high temperature and its thickness increases basically in an exponential way. As the applied load becomes heavier, the tribofilm increases first, yet it will decrease if the load is too high due to the dominant tribofilm removal process. To catch the variation of the tribofilm thickness under the severe contact state, a modified removal equation is proposed and it fits the experimental results well. 2. In the present study, as the tribofilm grows, the friction force will decline attributing to the lower CoF of the tribofilm, comparing with the CoF of the substrate. The friction force shows an decreasing tendency with the rising temperature from 40 to 140 °C, which are opposite to the simulation results without tribofilms. 3. The wear depth keeps decreasing with the rising temperature, and the effect on wear resistance of the ZDDP tribofilm reaches its maximum under the highest lubricant temperature. On the contrary, the wear depth will be higher under a heavier applied load, which is mainly attributed to the adhesive wear as the mild wear only count a small portion of the total wear. 4. The simulation results of CoF and wear depth using the modified tribofilm removal equation are well agreed with the experiments, which demonstrates that the present boundary lubrication model is validate.