Open access peer-reviewed chapter

Atomistic Simulation of Severely Adhesive Wear on a Rough Aluminum Substrate

Written By

Jun Zhong

Submitted: June 18th, 2020 Reviewed: September 14th, 2020 Published: October 27th, 2020

DOI: 10.5772/intechopen.94025

Chapter metrics overview

392 Chapter Downloads

View Full Metrics


In this Chapter, a severely adhesive wear on a rough aluminum (Al) substrate is simulated by molecular dynamics (MD) under a high velocity impact of a hard-asperity (a hard-tip) with the Al-asperity. Multiple simulations include effects of four factors: the inter-asperity bonding, the geometry overlap between two asperities, the impact velocity between two asperities and the starting temperature of the Al-substrate. It is observed that the deformation mechanism on the Al-substrate would involve a local melting (from 1200 to 2500 K) which forms liquid type layers (amorphous textures) in the contact area between two asperities. Also, temperature profiles on the hard-tip and the Al-substrate is depicted. Moreover, a method in the Design of Experiments (DOE) is employed to interpret above all simulations. The DOE results indicate that the inter-asperity bonding and the geometry overlap between two asperities would substantially increase the wear rate (for about 53.56% and 67.29% contributions), while the starting temperature of the Al-substrate and the impact velocity between two asperities would play less important roles (about 10.30% and 6.61%) in raising the wear rate.


  • the wear rate
  • molecular dynamics
  • the EAM potential
  • the Lennard-Jones potential
  • the design of experiment

1. Introduction

Surface wear of metals is very important to many industries such as automobile and aerospace etc. The wear rate often limits the lifetime and the durability of machinery parts, and thus leads to major economic losses. Wear phenomenon often involves the contacts between two or among more asperities on material surfaces, resulting in breaking old and forming new atomic bonds and plastic deformation in the contact areas. At the atomic scale, if only a few asperities come into contacts, the actual contact areas are very small when comparing with the macroscopic contact size. Thus local stresses in these areas can be exceptionally high, leading to high degrees of localized plastic deformation and heat generation, and even possibly local melting among asperities. In adhesive wear, aluminum (Al) wear is especially important because the Al is a relatively soft metal which is highly reactive with oxygen. Generally, a fresh (newly-formed) aluminum surface has little or no protection from oxides, and is less stable than the alumina (Al2O3). Therefore, excessive stresses and temperatures in the Al contact areas can provide an activation energy to initiate extremely exothermic reactions of fresh Al surfaces with oxides. For example, during the rolling of an aluminum sheet by a steel roller, if excessive stresses are applied onto the sheet surface without lubricant coverings, the newly-formed (fresh) Al surface may easily react with oxides on the steel roller, resulting in the local melting of aluminum to bond onto the roller surface, a so-called severely adhesive wear which causes the catastrophic breakage of the steel roller.

During past six decades, studies on the Al wear under a dry-sliding (no lubricants) constraint have revealed that the large plastic strain would occur near sub-surfaces on the Al-substrate when the wear process took place. Some experimental observations suggested that the Al wear rate be inversely proportional to the Al hardness because the higher Al hardness usually led to less plastic deformation on the Al-substrate [1]. Therefore, the understandings of the wear process may provide some valuable information for mechanisms of friction, lubrication and adhesion at the nano-scale [2, 3, 4, 5, 6, 7]. Even if, during asperity contacts at the atomic-scale, it was still difficult to observe wear mechanisms in nano-seconds.

In theoretical study, molecular dynamics (MD) may simulate the nano-scale phenomena in a very short time. Thus during past four decades, advances in the MD simulations have helped researchers understand atomic mechanisms which brought two kinds of materials into contact. For examples, Landman et alused the MD simulations to observe the hard-tip (Ni) jump-to-contact, the plastic yielding, the adhesion to induce the atomic flow and the slip generation in the Au-substrate [8, 9]. Plimpton et alused the MD simulations to observe the nucleation of partial dislocation loops occurring within the contact areas where a displacement-controlled hard sphere indented into a gold (Au) substrate surface. They found that the dislocation loops would grow rapidly into the substrate, but emerge at the surface edges, and then dislocation slips may produce complicated structures in the substrate [10, 11, 12]. Tanaka et alran the MD simulations to observe that, during the two−/three-body sliding contacts, an amorphous phase transformation would take place on the silicon substrate, i.e., the deformation on the silicon substrate would fall into adhering and plowing, but no wear regimes [13, 14]. Mendelev et alemployed the MD simulations to observe that, during a flat ruthenium (Ru) slab downward into a gold (Au) substrate with a single asperity, the Au was very ductile at 150 and 300 K, while the Ru showed the considerably less plasticity at 300 and 600 K [15]. In our former works, we have ran several MD simulations on the Al deformation at the nano-scale, including the wear [16] and the nano-indentation on the Al-substrates [17]. However, our these studies adopted a very low strength of inter-asperity bonding between the Al-substrate and a hard-tip, in which we found that, even if there was a large plastic deformation on the Al-substrate, no Al atoms were removed from the Al-substrate if the inter-asperity bonding was below a critical value.

To summarize, although many interesting MD simulations for deformation and wear on metal surfaces have been discussed at the atomic scale, there have not yet any MD simulations to focus on investigating a severely adhesive wear between two sliding surfaces. In the actual manufacturing, this kind of wear would occur during the rolling of aluminum sheet and many other forming processes. So in this chapter, such the wear will be discussed by the MD simulations to find out what may occur when a soft Al-asperity on an Al slab is contacted by another asperity on a hard tool surface when these two surfaces are sliding relative to another. In our MD investigation, multiple MD simulations were conducted by varying some constraint factors: the inter-asperity bonding, the geometry overlap between two asperities, the relative impact velocity between two asperities and the starting temperature of the Al-substrate etc. In details, the Al-substrate was simulated by an EAM potential [18, 19, 20], while the hard-asperity was simulated by a strong Lennard-Jones (L-J) potential which may serve as a model for hard tool surfaces (with their oxides). Moreover, different L-J potentials from weak to strong values may describe the bonding strength between the Al- and the hard-asperities, which may reveal an effect of the inter-asperity interactions on the severely adhesive wear. In Section 2, we describe our methodology. In Section 3, we list the simulation procedures. In Section 4, we assess the results of our MD simulations. Finally, in Section 5, we summarize the results.


2. Methodology

2.1 Molecular dynamics

Molecular dynamics (MD) is a methodology which depicts motions of a many-particle system using classical Newton’s equations. It has been over 40 years since the first MD’s application to a hard sphere system by Alder and Wainwright [21, 22]. The MD method is particularly useful for studying dynamical properties of materials, and help researchers extract physical-insights from the modeling works.

In the MD simulations, for a finite size system, classical trajectories of particles in real space are traced by solving the Lagrange equations numerically: when choosing the Cartesian coordinates, these equations would become Newton’s equations


where Vtotis the total potential energy of the system, miis the mass of particle i, Piis the momentum of particle i, and Fiis the total inter-particle force acting on particle i, and the “” denotes the first order of the time derivative. In practice, the convenience of using the MD simulation would mainly depend upon some particularly arithmetic algorithms to solve Eq. (1).

2.2 Algorithms and statistical ensembles

A standard way to solve Eq. (1) was the finite-difference method: given a configuration of a system (positions and velocities of all particles) and other dynamical information at time t, the numerical integration would determine the new configuration of the system at a later time t + Δt(Δtis the time step). Commonly used methods in the MD simulations were the Verlet algorithm[23], the Leapfrog algorithm[24], and the Gear predictor–corrector(GPC) algorithm[25]. An ideal algorithm should be simple, run fast, require little memory and permit using a long time step Δtto hold the whole system trajectories as true as possible, and preserve conservation laws of momentum and energy. According to these, the Verlet algorithmcan be regarded as the most widely used algorithm. The leapfrog algorithmwas essentially identical to the Verlet algorithm. The GPCmethod was usually more accurate as well as more complicated than others.

Commonly used statistical ensembles in the MD simulations were the microcanonical ensemble(NVE), the canonical ensemble(NVT), the isothermal-isobaric ensemble(NPT) and the grand canonical ensemble(μVT). Please note, thermodynamic variables in parenthesis for each ensemble were fixed during the MD simulations. Usually, the (NVE) ensemble was most convenient to realize since all equations of motion conserve the total energy Eand the particle number N, and the constant volume Vwas fixed by using periodic boundary conditions. The (NVT) ensemble was very commonly used for practical calculations. However, it needed to control the temperature: T. So, there were a number of ways to realize it: (1) the simplest one was the velocity scaling method, which was simple but crude in the MD simulations [26]. (2) the Nosé-Hoover’s method was more conceptually fundamental, but had an adjustable parameter of thermal inertia in practice [27]. For the (NPT) ensemble, Anderson [28] proposed an extended system method for a constant pressure MD. In this method, a term of “external piston” was added into the Lagrange functional for adjusting the value of external pressure. In another way, Parrinello and Rahman [29] proposed “a variable cell constant MD scheme” which allowed the simulation box to change its shape as well as its size. This method, so far, is the most general constant pressure method in the MD simulations. The (μVT) is very useful for an open system: in this ensemble, the total number of particles: Nwill not be conserved. Regarding this, see Ref. [30] for details.

2.3 Inter-particle potentials

Traditionally, the total potential energy of a system Vtotin the MD simulations was calculated using empirical potentials. There have been many kinds of models like pair-wise potentials: Lennard-Jones or Morse potentials; multi-body interactions [26], etc. In early 1980s, M.S. Daw and M.I. Baskes [18, 19] proposed a model for semi-empirical potential: the embedded-atom method(EAM) which reflected the many-body effect in materials. Generally, this model was written as two following terms


where rijis the distance between atoms iand j. In Eq. (2), on the right-hand side, the first term reflects a many-body effect in nature, meaning to embed ionic-cores into electron gas, so the Fi(ρi) was termed the embedding function; the second term represents the screening Coulombic repulsion among ionic-cores, which is a sum of pair-wise interactions. The total atomic-charge density ρiat the position rican be expressed as


where the ρiwas approximate to a sum of individual atomic-charge density fjaround site ineighbors, which was opposed to true charge distributions that resulted from the self-consistent ab-initiocalculations. During past four decades, there have been several fitting techniques to successfully constitute the EAM model, among which a particularly reliable method created by Adams et alwas developed by fitting to a large database of density functional calculated forces and experimental data, a so-called “the force-matching method,” see Ref. [20].

In this work, since we were primarily interested in the wear deformation on an Al-substrate, we chose a reliable aluminum EAM potential developed by Liu, Ercolessi and Adams through the force-matching method [31]. This potential has proven reliable for many types of bulk and surface simulations [20, 31, 32]. Also, since little deformation was assumed to occur in the contact tool surface (steel roller surface), we used a simple L-J potential with a very high bond strength to describe the hard-tip (the LJ-tip, hard-asperity). Similarly, another simple L-J potential was applied to the interaction between an Al-asperity and the hard-tip, in which the strength of inter-asperity bonding would vary from weak to strong in order, so as to simulate a wide range of generic systems. The L-J potential was adopted in a simple pair-wise formula as shown below


where rwas the distance between two atoms, εand σwere parameters corresponding to the bond strength and the bond length between two atoms, respectively. Table 1 provides relevant fitting parameters of those L-J potentials.

Bond strengthε(Å)σ(Å)
an Al-saperity with a LJ-tip0.20 (weak)2.85
0.65 (medium)2.85
1.00 (strong)2.85

Table 1.

Parameters of the Lennard-Jones potentials for interactions between an Al-asperity and a hard-tip (LJ-tip).

2.4 An implement simulation package: LAMMPS

The Large-scale Atomic/Molecular Massively Parallel Simulator(The LAMMPS) is a molecular dynamics package from the Sandia National Laboratories (see, to make use of the MPI for parallel communication. The parallel algorithm for the LAMMPS was a force-decomposition method in which a subset of pair-wise force computations was assigned to each processor. To improve the computational efficiency, the LAMMPS used neighbor lists to keep tracking particles nearby, which were optimized for systems with particles that were repulsive at short distances, so that the local density of particles never became too dense. The LAMMPS is a free software with open-source, distributed under the terms of the GNU General Public License. Till now, the LAMMPS is particularly efficient (in a parallel computing sense) to systems whose particles fill a 3D rectangular box with approximately uniform density.


3. Simulation procedures

3.1 Geometry constructions

Figure 1 shows schematic models for the MD simulations. They consisted of a hemispherical hard-tip (the LJ-tip) with 4440 L-J atoms, plus one slab with a sinusoidal asperity totally with 45,299 Al atoms (the Al-substrate). Table 2 provides the geometric information for the Al-substrate and the hard-tip at 0 K, both of which were perfect FCC single crystals.

Figure 1.

The MD simulation models for the Al-substrate and the LJ-tip.

ComponentShapeOrientationLattice (0 K, Å)Atom numbers
the LJ-tipHemisphereFCC4.0094440
the Al-asperitySinusoidFCC4.0324299
the Al-slabCubicFCC4.03241,000

Table 2.

Geometric information in the MD simulations.

3.2 Constraints in the MD simulations

These MD simulations were carried out through the LAMMPS which may apply the EAM and the L-J potentials to pure and alloy metal systems [10]. In each of these simulations, the lattice constant (a0 )of the Al bulk was determined at each desired temperature by the Hoover’s thermostat method [30]. And then, the simulation system was constructed by an appropriate a0 equilibrated for approximately 20,000 time steps (a time step = 0.001 ps) at one desired temperature. Bottom edge layers of the Al-substrate were fixed to prevent motions of the substrate during the whole MD simulations. Uppermost layer of the LJ-tip was constrained to move at a constant velocity along the x-direction (translational motion parallel to the Al slab surface), with no forces acting on this layer. All the MD simulations of asperity-asperity shear were performed at the constant energy (except for the fixed atoms), which allowed natural and realistic heat generation and diffusion due to formation of adhesive bonding and mechanical deformation in the contact area between two asperities.

3.3 Factors selected for design of experiments

Design of experiments (DOE) is a method which may be applied to analyses for experimental data, so that the relations of different selected factors (such as the strength of asperity–asperity bond) with response variables (such as the wear rate) can be determined. This method is widely applied in different fields to help researchers understand the factors which impact experimental processes most, so as to improve the desired experimental results [33].

During the MD simulations, the wear rate of aluminum was defined as the number of aluminum atoms removed from the Al-substrate, and was dependent upon many possibly significant factors. In order to get a better understanding of the wear process, some initial screening simulations were fulfilled to determine the most important factors and the range of these factors for investigation. Based upon the screening results, four factors were believed the most important ones: the inter-asperity bonding (DI), the impact velocity (VI) of the LJ-tip toward the Al-asperity along x-direction, the starting temperature of the MD system (TI), and the geometry overlap between two asperities (OI), see Table 3. The selected design technique was a so-called “The 24 full factorial design” with one central point for each factor (the mid-point of the range) in the Design Of Experiment (DOE). In this design, we did not consider any other replicates or blocking because they would not offer additionally valuable information due to the nature of computer simulations. So the chosen factors, their levels and their central points are listed in Table 3. In this Table, TI = 100 K, 400 K and 700 K; VI = 1.00 Å/ps, 5.50 Å/ps and 10.00 Å/ps; OI = 1.5 a0, 2.5 a0 and 3.5 a0; and DI = 0.20 eV, 0.65 eV and 1.10 eV. The velocities (VI) were corresponding to 100 ∼ 1000 m/s, which were relevant to the surface edge of steel roller and some other forming processes; these high speeds were also chosen because they were more computationally efficient to the MD simulations. Factors that were held fixed in the DOE analyses: crystal defects (single crystals with no defects), crystal orientations, the LJ-tip moving direction, the MD model size and geometry, the LJ-tip geometry, and the contact load, etc. So totally, 17 MD simulations were carried out for different combinations of four selected factors, and with combinations chosen through the software of DesignExpertTM Version 6 [33].

FactorVariableData typeLow (−)Central (0)High (+)
ATI (K)numerical100400700
BVI (Å/ps)numerical1.005.5010.00
COI (Å)numerical1.5 a02.5 a03.5 a0
DDI (eV)numerical0.200.601.00

Table 3.

Four selected factors and their levels in the DOE analyses.

3.4 Wear rate as a response variable

Since the wear rate is defined as the loss of materials from a rough substrate surface due to an interaction of the surface with its contacting environment, according to this, we chose the total number of Al atoms removed from the Al-asperity during asperity-asperity shear as the response variable.


4. Results and discussions

4.1 Wear rate

During a typical MD simulation of asperity-asperity shear, three different simulation stages took place: contact of two asperities, plowing of the LJ-tip through the Al-asperity, and final necking/fracture of the Al-asperity, as depicted in Figure 2. At the first stage, top of the LJ-tip impacted the Al-asperity, leading to the local heating, adhesion and plastic deformation at the contact area between two asperities, which was generally the region of maximum heat generation. At the second stage, the LJ-tip plowed through the Al-substrate, forming an Al neck to glide along a favorable slip system: . Finally, at the third stage, the Al neck was stretched and fractured, resulting in some Al atoms being removed from the Al-substrate. The extent of the Al necking was very sensitive to the VI and the TI. However, at the third stage, if the DI reached below a critical value (approximately 0.04 eV), the removed Al atoms would slip off the LJ-tip surface. And then they would return to the Al residual substrate, resulting in no net removal of Al atoms [16].

Figure 2.

Three stages of MD simulation for asperity-asperity shear: (a) the LJ-tip contact; (b) the LJ-tip plowing; (c) the Al necking and fracture.

Furthermore, occurrence of brittle or ductile fracture of the Al substrate during the MD simulations was observed to depend primarily upon the VI and the TI. That was, if the VI was high but the TI was low, there was very little or no elongation during the Al necking, and the final fracture of the Al neck was brittle. Conversely, if the VI was low while the TI was high, the Al neck would be deformed and elongated more, and the final fracture of the Al neck was more ductile.

4.2 Effects of four selected factors

Outputs of the MD simulations indicated that the number of Al atoms removed from the Al-substrate (wear rate) varied from 262 to 3144, depending upon different factors in the MD simulations. Analyses for these data indicated that, increasing the DI and the OI would increase the wear rate mostly, increasing the TI would just slightly increase the wear rate, whereas raising the VI exhibited a weak inverse correlation with the wear rate. In details, (1) if the DI became stronger, two asperities would adhere together via solid–liquid welding in the contact area, and more bonding states would occur along with more heat generation to increase temperature profiles in the contact area. So, a stronger DI may result in more removal of Al atoms. Due to thermal diffusion, temperature profiles in other areas of two asperities were also increased; (2) if enlarging the OI which determined the contact area and hence controlled the volume of Al atoms being plowed by the LJ-tip, the deformation of the Al-substrate would increase. And then, it resulted in a higher wear rate. Komanduri et alalso reported similar results in their studies [9]; (3) if increasing the TI, dislocations in the Al-substrate can move and glide more easily through the elevated temperature region, resulting in large plastic deformation, i.e., the transformation from the brittle fracture (a small Al necking) at lower temperatures to the ductile fracture at higher temperatures can be easily observed; (4) it has found that the VI would play a minor inverse correlation role in the wear rate. This could be due to the increase of strain rate along with the rise of the VI, which decreased the ductility of materials so that the wear process on the Al-substrate may undergo less plastic deformation before fracture. As discussed earlier, more plastic deformation at the stage of the LJ-tip plowing may enlarge the area of adhesion. Therefore, when less plastic deformation took place, fewer Al atoms were removed by the LJ-tip. Similar experiments on dry-sliding wear of Al-Si alloys indicated that the coefficient of wear in these alloy systems was highly dependent upon the disc speed: at a faster disc speed (0.356 m/s), the wear rate was found to be low or moderate; While at a slower speed (0.089 m/s), the wear rate increased dramatically [34]. So, it should be noticed that, in experiments when the VI increased, more asperity-asperity collisions could take place, whereas in our simulation, only two single asperities were simulated. Thus when comparing with experimental observations, one would have to consider both the damage per asperity and the increased number of two interacting asperities.

4.3 Plastic deformation during asperity-asperity shear

It was well known that, during the wear process, severe temperature fluctuations could result in major changes of mechanical properties. That was, during a higher temperature region, moduli and yield strength of metals could usually decrease, so the localized plasticity would increase and the welding could even occur if a critical temperature was reached [35]. It was found that, near the contact area of two asperities, the Al neck was deformed amorphously under very high stresses and temperatures. In this area, linear motion of debris along the sliding direction was transformed into thermally random atomic motions. Therefore, the higher the VI, the much higher the temperature profiles in the contact area became, so a thin and soft/liquid like layer may form. As a result, the deformed Al neck would behave like a viscous liquid under very high temperatures and stresses [36]. When at a very high VI (= 10.0 Å/ps), amorphous plasticity became the major deformation mechanism. However, at a very low VI (= 1.00 Å/ps) and other lower TI, DI and OI, in addition to amorphous deformation at the contact area, dislocation cores were identified by using the evaluation of the centrosymmetry parameter [10], so they were found to emit into the Al-substrate from the high temperature region (amorphous deformation region) along a favorable system 111101¯, see Figure 3.

Figure 3.

Amorphous deformation of the Al neck and emission of dislocation cores in the Al-substrate.

4.4 Thermal analysis

Thermal distributions from four different MD simulations are shown in Figures 4 and 5, which give insight into heat generation and heat transfer during asperity-asperity shear. Temperature profiles of the Al-substrate were calculated by using the following equation:

Figure 4.

(a). side view of the LJ-tip at the contact stage; (b) slide view of the LJ-tip at the plowing stage; (c) temperature profiles for a slice of the LJ-tip at the final stretching stage. (TI = 700 K,OI = 3.5a0,DI = 0.20 eV andVI = 10.0 Å/ps).

Figure 5.

Thermal analyses for two different MD simulations underTI = 700 K,OI = 3.5a0,DI = 1.00 eV, and (a)VI = 1.00 Å/ps and (b)VI = 10.0 Å/ps.


where mk, Tiand vkwere mass, temperature, and velocity for atom iand k, respectively. kBwas the Boltzmann constant, and nwas the number of atoms for the Al substrate and the LJ-tip within a sphere of about 27 Å in diamond at a point iconsidered [16]. For the moving LJ-tip and the Al atoms removed by the LJ-tip, the temperature profiles were calculated by


where VIwas the impact velocity of the LJ-tip along the x-direction.

For examples, Figure 4 shows three selected views of narrow slice for the LJ-tip during the MD simulations along with relevant temperature profiles. In Figure 4(a) and (b), atoms in the first layer on the Al-asperity were found to follow almost same stacking sites as those L-J atoms being located at the top of the LJ-tip when two asperities came into bond together. And then in Figure 4(c), the order of subsequent layers in the Al-asperity became more random in high temperature region (about 1200 ∼ 1550 K): atoms in these layers moved randomly around the top of the Al-asperity (local melting) until they bonded to and transferred heat to the cooler LJ-tip.

Figure 5 shows the temperature profiles for two selected systems: TI = 700 K, OI = 3.5 a0, DI = 1.00 eV, and VI = 1.00 Å/ps, 10.0 Å/ps, respectively. Figure 6 shows the temperature profiles for two selected systems: TI = 100 K, OI = 3.5 a0, DI = 1.00 eV, and VI = 10.0 Å/ps, 1.00 Å/ps, respectively. In these two Figures, (1) at the contact stage of the LJ-tip, the local heating occurred as new atomic bonds formed in the contact area between two asperities. Thus for the LJ-tip, the temperature gradient was positive from the contact area to its top layers; while for the Al-substrate, this positive gradient was found from the deformed Al neck to its far ends; (2) at the plowing stage of the LJ-tip, the heat diffused into the LJ-tip bulk from the contact area; while for the Al-asperity, the Al necking also produced more heat in the remnant surface where the Al necking root glided, so the heat diffused into the Al-substrate from the remnant surface, and hence increased the temperature profiles there.

Figure 6.

Thermal analyses for two different MD simulations underTI = 100 K,OI = 3.5a0,DI = 1.00 eV, and (a)VI = 1.00 Å/ps and (b)VI = 10.0 Å/ps.

In these two Figures, the maximum local temperature and the temperature gradient were quite different: the heat generation was much faster than the heat diffusion under a very high VI, so the VI seemed to play a major role in generating the maximum local temperature. In addition, a higher local temperature may lead to a larger local softening of the Al-substrate, so the deformation was mostly limited to a narrow region, and hence materials were removed. It was noticed that the maximum local temperature was much higher than the starting temperature TI, so the TI seemed not to play a major role in removing the Al atoms. It should also be noticed that the maximum local temperature can briefly exceed the boiling point of the Al, and yet the local liquid was not boiling. The reason was that, the increase of temperature was very brief (about 10–100 ps), while the activation energy for atoms to jump into the gas phase was much higher (approximately the cohesive energy), so the kinetics prevented any significant amount of evaporation from occurring.

4.5 Statistical analysis

Seventeen outputs for the MD simulations are listed in Table 4. The goal of the statistical analysis by the DOE is to determine the “best fit” equation to describe the wear rate as a function of four simulation variables (inter-asperity bonding, geometry overlap, impact velocity and starting temperature). It was assumed that the effect of each variable is additive, and there is no interaction between each of two variables (no cross term). So, this simplest model may work very well for the DOE analysis [16]. However, because the wear rate varied over a large range (the ratio of wear rate from its maximum to minimum was about 12), the wear rate must by transformed by using a natural log. Therefore, the equation for the wear rate in terms of the selected four factors was expressed as follows

TI (K), A-termVI (Å/ps), B-termOI (Å), C-termDI (eV), D-termWear rate (atom number)
1001.001.50 a00.20409
7001.001.50 a00.20453
10010.01.50 a00.20262
70010.01.50 a00.20292
1001.003.50 a00.20947
7001.003.50 a00.201231
10010.03.50 a00.20842
70010.03.50 a00.20982
1001.001.50 a01.00711
7001.001.50 a01.001433
10010.01.50 a01.00694
70010.01.50 a01.00933
1001.003.50 a01.001346
7001.003.50 a01.003144
10010.03.50 a01.001402
70010.03.50 a01.001855
4005.002.50 a00.501244

Table 4.

Seventeen combinations of the MD simulations for the DOE analyses.


where the Wear(wear rate) was the number of Al atoms removed by the LJ-tip, Awas the starting temperature (K), Bwas the impact velocity (Å/ps), Cwas the geometry overlap (Å), and Dwas the inter-asperity bonding (eV). Please note, Eq. (7) applied an exponential formula to describing the wear rate because of the large variation in the wear rate (see the detailed analyses of variance for the selected factors in Table 5).

SourceSum of squareDFMean squareF values
Prob > F
Model6.3341.5834.44< 0.0001
C3.0913.0967.29< 0.0001
D2.4612.4653.56< 0.0001
Core Total6.8816

Table 5.

ANOVA for those selected factorial models provided by the design expert™ software.

Please note: ATI; BVI; COI; DDI.

It is more useful to describe the wear rate using the normalized variables (variables are normalized to a scale from −1 to +1, the low and high levels for each variable) because the magnitude of A, B, Cand Dcoefficients allows one to easily determine the relative importance of each variable.

Through using the normalized variables, Eq. (7) can be rewritten as follows


In Eq. (8), coefficient magnitudes of variables indicated that A(= +0.19) term had a small effect, B(= −0.12) term had a minor inverse correlation, while C(= +0.42) and D(= +0.37) terms had major effects on the wear rate. That was, the inter-asperity bonding (D) and the geometry overlap (C) had much more effects on the wear rate than the starting temperature (A) and the impact velocity (B).


5. Conclusions

In this work, a severely adhesive wear was investigated by simulating asperity-asperity shear between a fast moving rigid LJ-tip toward an Al-asperity. Molecular dynamics simulations were conducted for 17 different combinations of four variables: the starting temperature, the relative velocity, the geometry overlap and the inter-asperity bonding between two asperities. It was found that the wear process occurred in three stages: the contact, the plowing, and the necking/fracture on the aluminum substrate. Thermal analyses indicated that the heat generated during the MD simulations stemmed from the adhesive reaction in the contact area between two asperities, and then in the remnant surface on the Al residual substrate where the Al necking root glided. Bond formation and mechanical deformation during asperity-asperity shear may result in large increases of local temperature in the contact area (1200 ∼ 2500 K), so the primary mechanism of deformation on the Al-substrate was due to amorphous plasticity and local melting. Generations and motions of dislocation cores were observed under milder conditions where little melting occurred. A method: The 24 full factorial designin the Design Of Experiment was adopted in analyzing effects of those four variables (factors) on the wear process. Analysis results indicated that, the inter-asperity bonding and the geometry overlap between two asperities would play much more important roles in the wear process than the starting temperature and the impact velocity.



This work was sponsored by the National Science Foundation (Grant no.: DMR9619353), USA; the General Program of The National Nature Science Foundation of China (Grant no.: 11972348), the Joint Fund Program of The National Nature Science Foundation of China (Grant no.: U1738108), and the Hebei Provincial Key Laboratory of Thermal Protection Materials (SZX2020038, in preparation), China.


  1. 1. Wert J.J, Singerman S.A, Galdwell S.G, Quarles R.A. The role of stacking fault energy and induced residual stresses on the sliding wear of aluminum bronze. Wear. 1983;91: 253–267. DOI:
  2. 2. Syed S.A, Wahl K.J, Colton R.J. Nanoindentation and contact stiffness measurement using force modulation with a capacitive load-displacement transducer. Review of Scientific Instruments. 1999;70: 2408–2413. DOI:
  3. 3. Engelder J.T, Scholz C.H. The role of asperity indentation and ploughing in rock friction—II: Influence of relative hardness and normal load. International Journal of Rock Mechanics and Mining Sciences & Geomechanics Abstracts. 1976;13: 155–163. DOI:
  4. 4. Putman C and Kaneko R. Experimental observation of single-asperity friction at the atomic scale. Thin Solid Film. 1996;273: 317–321. DOI:
  5. 5. Maw W, Stevens F, Langford S.C, Dickinson J.T. Single asperity tribochemical wear of silicon nitride studied by atomic force microscopy. Journal of Applied Physics. 2002;22: 5103–5109. DOI:
  6. 6. Xue X, Polycarpou A.A, Phinney LM. Measurement and Modeling of Adhesion Energy Between Two Rough Microelectromechanical System (MEMS) Surfaces. Journal of Adhesion Science and Technology. 2008;22: 429–455. DOI:
  7. 7. Mate C.M, McClelland G.M, Erlandsson R, Chiang S. Atomic-scale friction of a tungsten tip on a graphite surface. Physical Review Letters. 1987;59: 1942–1945. DOI:
  8. 8. Landman U, Luedtke W.D, Nancy A, Colton R.J. Atomistic Mechanisms and Dynamics of Adhesion, Nanoindentation, and Fracture. Science. 1990;248: 454–461. DOI: 10.1126/science.248.4954.454
  9. 9. Komanduri R, Chandrasekaran N, Raff L.M. MD simulation of indentation and scratching of single crystal aluminum. Wear. 2000;240: 113–143.
  10. 10. Kelchner C.L, Plimpton S.J, Hamilton J.C. Dislocation nucleation and defect structure during surface indentation. Physical Review B. 1998;58: 11085.
  11. 11. Saraev D, Miller R.E. Atomistic simulation of nanoindentation into copper multilayers. Modelling and Simulation in Materials Science and Engineering. 2005;13: 1089–1100. DOI:
  12. 12. Zimmerman J.A, Kelchner C.L, Hamilton J.C, Foiles S.M. Surface Step Effects on Nanoindentation. Physical Review Letters. 2001;87: 165507. DOI:
  13. 13. Zhang L, Tanaka H. Towards a deeper understanding of wear and friction on the atomic scale—a molecular dynamics analysis. Wear. 1997;211: 44–53. DOI:
  14. 14. Zhang L, Tanaka H. Atomic scale deformation in silicon monocrystals induced by two-body and three-body contact sliding. Tribology International. 1998;31: 425–433. DOI:
  15. 15. Fortini A, Mendelev M.I, Buldyrev S, Srolovitz D. Asperity contacts at the nanoscale: Comparison of Ru and Au. Journal of Applied Physics. 2008;104: 074320. DOI:
  16. 16. Zhong J, Adams J.B, Hector L.G. Jr. Molecular dynamics simulations of asperity shear in aluminum. Journal of Applied Physics. 2003;43: 4306–4314. DOI:
  17. 17. Yu H, Adams J.B, Hector L.G. Jr. Molecular dynamics simulation of high-speed nanoindentation. Modelling and Simulation in Materials Science and Engineering. 2002;10: 319–330. DOI:
  18. 18. Daw M.S, Baskes M.I. Semiempirical, Quantum Mechanical Calculation of Hydrogen Embrittlement in Metals. Physical Review Letters. 1983;50: 1285–1288. DOI:
  19. 19. Daw M.S, Baskes M.I. Embedded-atom method: Derivation and application to impurities, surfaces, and other defects in metals. Physical Review B. 1984;29: 6443–6453. DOI:
  20. 20. Ercolessi F, Adams J.B. Interatomic Potentials from First-Principles Calculations: The Force-Matching Method. Europhysics Letters. 1994;26: 583–588. DOI:
  21. 21. Alder B.J, Wainwright T.E. Studies in Molecular Dynamics. I. General Method. The Journal of Chemical Physics. 1959;31: 459–466. DOI:
  22. 22. Alder B.J, Wainwright T.E. Studies in Molecular Dynamics. II. Behavior of a Small Number of Elastic Spheres. The Journal of Chemical Physics. 1960;33: 1439–1451. DOI:
  23. 23. Verlet L. Computer “Experiments” on Classical Fluids. I. Thermodynamical Properties of Lennard-Jones Molecules. Physical Review. 1967;159: 98–103. DOI:
  24. 24. Beeman D. Some multistep methods for use in molecular dynamics calculations. Journal of Computational Physics. 1976;20: 130–139. DOI:
  25. 25. Gear G.W.Numerical Initial Value Problems in Ordinary Differential Equations(Prentice-Hall Inc., Englewood Cliffs, New Jersey, 1971). ISBN-13: 978–0136266068
  26. 26. Heermann D.W.Computer Simulation Method in Theoretical Physics(Springer Verlag, Berlin, 1990). ISBN: 978–3–540-52210-2
  27. 27. Nosé S. A Molecular Dynamics Method for Simulations in the Canonical Ensemble. Molecular Physics. 1984;52: 255–268. DOI: 10.1080/00268978400101201
  28. 28. Anderson H.C. Molecular dynamics simulations at constant pressure and/or temperature. The Journal of Chemical Physics. 1980;72: 2384–2393. DOI:
  29. 29. Parrinello M, Rahman A. Polymorphic transitions in single crystals: A new molecular dynamics method. Journal of Applied Physics. 1981;52: 7182–7190. DOI:
  30. 30. Hoover W.G.Computational Statistical Mechanics(Elsevier Science, 1991). ISBN: ISBN-13: 978–0444881922
  31. 31. Liu X-Y, Ercolessi F, Adams J.B. Aluminium interatomic potential from density functional theory calculations with improved stacking fault energy. Modelling and Simulation in Materials Science and Engineering. 2004;12: 665–670. DOI:
  32. 32. Liu X-Y, Adams J.B, Ercolessi F, Moriarty J.A. EAM potential for magnesium from quantum mechanical forces. Modelling and Simulation in Materials Science and Engineering. 1996;4: 293–304. DOI:
  33. 33. Montgomery D.C. Design and Analysis of Experiment (8th ed., Wiley, 2012). ISBN-13: 978–1118146927
  34. 34. Lasa L, Rodriguez-Ibabe J.M. Effect of composition and processing route on the wear behaviour of Al–Si alloys. Scripta Materialia. 2002;46: 477–481. DOI:
  35. 35. Markov D, Kelly D. Mechanisms of adhesion-initiated catastrophic wear: pure sliding. Wear. 2000;239: 189–210. DOI:
  36. 36. Ling F.F, Edward Saibel E. Thermal aspects of galling of dry metallic surfaces in sliding contact. Wear. 1957;1: 80–91. DOI:

Written By

Jun Zhong

Submitted: June 18th, 2020 Reviewed: September 14th, 2020 Published: October 27th, 2020