- Email: [email protected]
An analysis method for atomistic abrasion simulations featuring rough surfaces and multiple abrasive particles
Accepted Manuscript An analysis method for atomistic abrasion simulations featuring rough surfaces and multiple abrasive particles S.J. Eder, D. Bianchi, U. Cihak-Bayr, A. Vernes, G. Betz PII: DOI: Reference:
S0010-4655(14)00183-0 http://dx.doi.org/10.1016/j.cpc.2014.05.018 COMPHY 5345
To appear in:
Computer Physics Communications
Received date: 7 February 2014 Revised date: 16 April 2014 Accepted date: 23 May 2014 Please cite this article as: S.J. Eder, D. Bianchi, U. Cihak-Bayr, A. Vernes, G. Betz, An analysis method for atomistic abrasion simulations featuring rough surfaces and multiple abrasive particles, Computer Physics Communications (2014), http://dx.doi.org/10.1016/j.cpc.2014.05.018 This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65
An analysis method for atomistic abrasion simulations featuring rough surfaces and multiple abrasive particles S. J. Edera,∗, D. Bianchia , U. Cihak-Bayra , A. Vernesa,b , G. Betzb a Austrian b Institute
Center of Competence for Tribology, Viktor-Kaplan-Straße 2, 2700 Wiener Neustadt, Austria of Applied Physics, Vienna University of Technology, Wiedner Hauptstraße 8–10 /134, 1040 Vienna, Austria
Abstract We present a molecular dynamics (MD) model system to quantitatively study nanoscopic wear of rough surfaces under two-body and three-body contact conditions with multiple abrasive particles. We describe how to generate a surface with a pseudo-random Gaussian topography which is periodically replicable, and we discuss the constraints on the abrasive particles that lead to certain wear conditions. We propose a postprocessing scheme which, based on advection velocity, dynamically identifies the atoms in the simulation as either part of a wear particle, the substrate, or the sheared zone in-between. This scheme is then justified from a crystallographic order point of view. We apply a distance-based contact zone identification scheme and outline a clustering algorithm which can associate each contact atom with the abrasive particle causing the respective contact zone. Finally, we show how the knowledge of each atom’s zone affiliation and a timeresolved evaluation of the substrate topography leads to a break-down of the asperity volume reduction into its components: the pit fill-up volume, the individual wear particles, the shear zone, and the sub-surface substrate compression. As an example, we analyze the time and pressure dependence of the wear volume contributions for two-body and three-body wear processes of a rough iron surface with rigid spherical and cubic abrasive particles. Keywords: atomistic abrasive wear, molecular dynamics, multiple particle abrasion PACS numbers: 46.55.+d, 31.15.xv, 68.35.Ct
1. Introduction The increasing miniaturization of electronic and mechanical devices requires certain surfaces to be finished to a high degree of smoothness. Even some macroscopic surfaces, such as those of molding parts for the optical industry, must be virtually scratch-free in order to produce high-quality optical products. As scratch dimensions enter the (sub-)nanometric length scale, experimental treatment and analysis of the workpieces becomes complicated, expensive, or even impossible. Similarly, with manufacturing tolerances becoming ever smaller, the wear mechanisms leading to geometrical changes in components have to be understood on the atomistic scale. Here molecular dynamics (MD) simulation can reveal important insights into the processes of interest [1]. In literature, MD simulations of two-body and three-body contact sliding are discussed either in terms of tribological phenomena such as abrasive wear, or as surface finishing processes such as grinding and polishing. An extensive overview of the simulation of the grinding process at several length scales can be found in Ref. [2], where it is noted that compared to nano-cutting, little MD literature is available on grinding. The first earnest attempts to explicitly simulate two-body and three-body wear at an atomistic level were undertaken in the mid-1990s [3, 4, 5]. Most MD simulations treat grinding as a nanometric ∗ Corresponding
author Email address: eder at ac2t.at (S. J. Eder)
Preprint submitted to Elsevier
April 16, 2014
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65
cutting process and focus on the interaction between one single abrasive particle and a usually flat substrate [6, 7, 8]. The atomistic wear of even a single roughness feature due to asperity-asperity interaction has only very recently received some attention, e.g., in Refs. [9, 10]. While these approaches can lead to an understanding of a process under controlled conditions, e.g., at constant penetration depth and orientation of the abrasive particle, both wear and finishing processes may involve circumstances such as changing rough topographies, random rake angle distributions or wear debris build-up between neighboring abrasives, which cannot be studied by considering only a single abrasive interacting with a flat surface. For these reasons, the implementation of a statistically based surface roughness, similar to what was recently used in sliding friction simulations [11], combined with multiple abrasive particles acting on the surface, is highly desirable. Literature on the atomistic modeling and simulation of three-body wear is scarce if one is not primarily interested in the chemo-mechanical planarization of silicon [8, 12, 13, 14]. Especially the extended degrees of freedom (rolling motion) of the particles have been considered in very few publications [15, 16]. An MD approach to the polishing of a silicon surface with several randomly distributed asperities and abrasive particles was proposed in Ref. [17], leading to a phenomenologically based material removal rate. In this work, we present an easily scalable and extendable MD model system which can be used to quantify the wear volume and the contact zone occurring during wear processes on a rough surface with multiple hard abrasive particles. The model systems were intentionally kept simple and do not include a medium or phenomena such as the wear of the abrasive particles. However, our approach is general enough to allow its application to the analysis of material removal rates in processes such as polishing and grinding as well. Note that a comprehensive discussion, from a materials science point of view, of the results obtained for our chosen model systems is beyond the scope of this publication. In Sec. 2, we discuss the model setup, the procedure for pseudo-random Gaussian topography generation as well as the constraints which lead to two-body or three-body contact sliding. Sec. 3 deals with the identification of wear particles based on their advection velocities. In Sec. 4 we define and quantify the contact zone between the abrasives plus the wear particles and the remaining substrate. Here we also discuss a clustering algorithm which allows the affiliation of every contact zone with an abrasive particle. Finally, in Sec. 5, we perform a mesh-based volumetric substrate topography analysis on which we base the wear volume estimation. The knowledge of the wear particles allows a break-down of the total asperity reduction volume into several contributions. By applying a modified version of the clustering algorithm, the wear particle volume is then broken down further into contributions due to single abrasive particles. 2. Model setup and simulation procedure The substrate consists of a block of bcc Fe with lateral dimensions of 28.5×28.5 nm2 and sufficient initial height to accommodate the pseudo-random Gaussian roughness distribution, the generation procedure of which will be discussed in more detail later. At its base, three monolayers of atoms form a rigid support, covered by two monolayers of Fe thermostatted to 300 K with a Langevin thermostat with a time constant of 0.5 ps, acting only in y direction in order not to overly interfere with the load acting in z direction and the movement of the abrasive particles, which is predominantly in x direction. Periodic boundary conditions are applied in the x and y directions, which means that if an atom leaves the simulation box through one boundary, it re-enters at the opposing one. No matter is lost during a simulation run. The abrasive particles are modeled as rigid bodies arranged in a staggered 4×4 grid, rotated randomly between 0 and 90 degrees about all three Cartesian axes, and are initially placed several nanometers above the surface. Two abrasive particle geometries are considered: spherical with a radius of 3.25 nm and cubic with a side length of 4.2 nm. While these extreme geometries are unlikely to occur in real processes, they were chosen for illustrative purposes due to their likelihood of producing clear differences in the results. For similar reasons, the areal density of the abrasive particles was substantially increased compared to realistic systems in order to speed up the generation rate of wear debris and thus reduce computational expense. Snapshots of the initial configurations are shown in Fig. 1. A Finnis-Sinclair potential is used for the Fe–Fe interactions [18]. The interactions between abrasive particles and the Fe surface is modeled using a Lennard-Jones potential with ε = 0.125 eV and σ = 2
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65
0.2203 nm, a parameter set comparable to Refs. [4, 5, 7, 19, 20]. As the abrasive spheres in the threebody contact simulations have the possibility to touch each other, we assumed a very weakly attractive Lennard-Jones potential for this interaction with the parameters ε = 0.03 eV and σ = 0.2203 nm. We built a periodically replicated surface topography with Gaussian height distribution. The parameters which influence the final topography are the roughness parameter Sq , which is the standard deviation from (Asp) the mean surface height zsubst [21] as well as the nominal asperity width dx/y . First, these nominal values are replaced with ones of similar size so that the simulation box may be filled in x and y direction with an integer number of the resulting cells. The cell centers constitute the control points of the topography and are assigned a random height (z) value such that the total distribution is as close to Gaussian as can be achieved with the relatively small number of control points. The resulting 2D-matrix of z values is then padded at each boundary with two rows (or columns) from the opposite side to account for periodicity. The respective x and y values of the support points are offset by random values between 0 and a maximum (offset) (Asp) offset dx/y < dx/y , and then also padded to produce a periodic boundary. For the initial topography used as an example throughout this work, we assigned the parameter values zsubst = 4 nm, Sq = 0.8 nm, (Asp) (offset) dx/y = 2 nm and dx/y = 1 nm. On a realistic surface, the typical height of a roughness feature would be approximately 2–3 orders of magnitude smaller than its lateral dimensions. As we are primarily interested in (Asp) the illustration of our analysis method, in this example Sq and dx/y were chosen to be of similar proportion, leading to considerably steeper asperity flanks. The scattered point cloud generated in this fashion constitutes a very coarse representation of the topography. In order to filter an initial set of substrate atoms according to this surface, it is necessary to know the z value at any given lateral position. To achieve this, the scattered data was smoothed and mapped to a fine mesh (of the order of the lattice parameter of Fe) using a modified ridge estimator algorithm [22]. The application of the resulting highly resolved topography as a criterion for surface generation yields a perfectly periodic atomistic topography with smooth transitions between roughness features. It should be noted that the originally set roughness of the surface, characterized by Sq , is slightly reduced by this smoothing. This must be taken into account if a very exact initial roughness is required. A brief relaxation of the substrate to a temperature of 300 K with repeated re-assignment of atomic velocities according to the Maxwell-Boltzmann velocity distribution leads to the initial state of the substrate used throughout this work which is shown in Fig. 1. When simulating two-body contact sliding, the abrasive particles are dragged across the surface in +x (max) (max) (max) (max) /8, resulting in a sliding velocity of v (max) ' vx = vx direction at vx and in +y direction at vy ◦ at an angle of 7.125 with the x-axis. This way the abrasives re-enter the simulation box at different y positions every time they pass through the periodic box boundary, thus scanning the surface more evenly and avoiding artificially amplified deep scratches. The relative positions and orientations of the abrasives are kept constant as would be the case for hard, rough counterbody such as a grindstone. The normal load is kept constant over time and distributed evenly between the 16 particles, which can change their z position collectively as required by the topography, but not individually. This may lead to loss of surface contact for individual abrasive particles for some time intervals, while the remaining particles share the additional load. In the simulation of three-body contact sliding, the abrasive particles are moved across the surface (max) parallel to the +x direction at vx . In contrast to the two-body contact simulations, here they may rotate about their respective centers of mass according to the atomic interactions with the substrate surface and the neighboring abrasive particles. The normal load is distributed evenly between the 16 particles, and each particle may individually change its y and z position as required according to the normal load, the sliding velocity, and the local topography, as would be the case for dust or debris in the lubrication gap of a tribosystem or in a polishing slurry. These additional degrees of freedom lead to a contact situation quite different from two-body contact sliding, since the number of load bearing regions is always equal to the number of abrasive particles, and the rolling motion of the abrasives leads to considerably less shear. We performed simulations of three-body wear processes with spherical abrasives and two-body wear processes with spherical or cubic abrasives at pressures of σz = 0.1, 0.5 and 1.0 GPa, and at sliding velocities of 8 and 16 m/s using the open source MD code LAMMPS [23]. The total simulation time was 3
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65
(a)
(b)
(c)
Figure 1: Initial simulation configurations. Abrasives are shown in green, the substrate up to the mean surface height in blue, and initial asperities in yellow. (a) 16 randomly rotated rigid cubic abrasive particles for two-body contact sliding simulations. (b) 16 rigid spherical abrasive particles for two-body and threebody contact sliding simulations. (c) Top view of initial surface topography, colored according to atomic z coordinate, where the peaks are red and the pits are blue. The abrasive particles are not shown.
4
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65
(max)
(max)
5 ns for v1 = 16 m/s and 10 ns for v2 of ∆x = 80 nm.
= 8 m/s to yield a total translation distance along the x axis
3. Identification of wear particles Fig. 2 shows top views of the system configurations without the abrasives after 5 ns of sliding at v (max) = 16 m/s, colored according to the atomic z coordinate. Although a comparison between the system in Fig. 2 (a), which was obtained with three-body contact sliding at low load, and the one in Fig. 2 (i), which is the result of two-body contact sliding with cubic abrasives at high load, clearly shows the formation of wear particles in the latter case, the automatic and time-resolved identification of the atoms belonging to this group is not so straightforward. The distinct horizontal grooves in Figs. 2 (b) and (c) occur due to interlocking of neighboring abrasives, leading to effective two-body contact sliding of one of the particles as its rolling motion is inhibited. The grooves are perfectly horizontal because in the three-body contact sliding simulations, the only kinematic constraint on the abrasives is the motion in x direction while the motion in y direction results mainly from the current local topography. If atoms were removed from the surface in such a fashion that there exists no Fe connection between them and the substrate, this detachment could serve as an indicator for wear debris. However, since this applies only to very few atoms, which cannot form a real wear particle, our proposed definition of wear particles is based on the advection velocity of the atoms. This method is akin to the contact atom counting procedure described in further detail in previous work by Eder et al. in Refs. [24] and [25], which is used to identify the contacting Fe atoms between two Fe sliders. Here, the main component of the advection velocity, v = vx , is calculated via the position differences ∆x evaluated at a sampling interval of ∆t = 20 ps. This signal is cleaned from the artifacts of the periodic boundary conditions and filtered using a Gauss filter with a total width of 0.8 ns and a standard deviation of 0.4 ns (for v (max) = 8 m/s) to yield continuous velocity profiles for each atom. Based on the instantaneous advection velocity thus calculated, each atom may be assigned to one of three groups: (a) atoms traveling at v > 0.9 v (max) are considered stuck to an abrasive and are therefore classified as wear particles, (b) atoms traveling at |v| < 0.1 v (max) are considered stationary and are classified as substrate atoms, and (c) all remaining atoms, which move faster than 0.1 v (max) (not necessarily with the same sign, e.g., in three-body contact sliding simulations), are classified as part of what we will henceforth refer to as the “shear zone”. Fig. 3 gives a simplified overview of the discussed zones. To justify this seemingly arbitrary atom identification scheme not only from an angle of computational ease, but also from a crystallographic one, we calculated and compared the time-averaged radial distribution functions (RDFs) of the wear particles, the shear zones, and the average bulk substrate, see Fig. 4. For the two-body contact sliding simulations, a comparison between the RDF for the average bulk, which exhibits the typical pattern of a bcc lattice with broadened peaks due to finite temperature, and that of the shear zone as identified using the above scheme, proves that the latter region has a considerably disturbed lattice. The coalescence of the first with the second and the third with the fourth peak suggest that this region is almost amorphous. The wear particles, on the other side, seem to crystallize onto the abrasives in such a way that they show almost the same pattern as the bulk substrate, the main difference being a height reduction of the second-neighbor peak. In the three-body contact sliding simulations, the shear zone is crystallographically much more similar to the wear particles. For both regions, crystallinity is generally retained, with only slight pairwise coalescence of neighboring peaks. However, these regions play a much smaller quantitative role in the three-body contact sliding simulations, since the rolling motion of the abrasives results in much less shear, and the non-cutting wear regime leads to few wear particles. We therefore consider our atom identification scheme sufficiently universal for application to all our simulations. Fig. 5 compares the sliding distance and the pressure dependence of the total number of atoms in wear particles and the shear zone for all considered wear processes. All pressure dependences are nearly linear. We can observe that the sliding velocity has little impact on the general shape of the curves. While simply monitoring the total number of wear particle atoms Nwp over time cannot prove that the formed wear particles are stable, a close scrutiny of simulation visualisations with coloring according to the advection velocity (not shown here) indicates that for configurations where appreciable numbers of wear particle 5
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65
(a)
(b)
(c)
(d)
(e)
(f)
(g)
(h)
(i)
Figure 2: Top views of surface topographies after 5 ns of sliding at 16 m/s. Colored according to atomic z coordinate, where the peaks are red and the pits are blue. Top: three-body contact sliding, center: twobody contact sliding with spherical abrasives, bottom: two-body contact sliding with cubic abrasives, left: 0.1 GPa, center: 0.5 GPa, right: 1.0 GPa.
6
Figure 3: Sketch of the zones to which the atoms can be attributed. atoms are produced, the respective wear particles are indeed stable and the shear zones (in which plastic flow occurs) are quite thin by comparison. What looks like a saturation of Nwp in Fig. 5 (a) actually marks the point where the running-in of the wear process has ended and the steady-state regime begins, a topic of ongoing research. In Fig. 5 (b), a sharp initial increase of the number of atoms in the shear zones takes place within the half-width of the employed Gauss filter and therefore cannot be resolved adequately. This is followed by a decrease over time as the asperity flanks gradually become less steep, with Nshear eventually becoming a constant proportional to the average shear force. 16 shear zone wear particle bulk average
14 12 10 g(r) [−]
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65
8 6 4 2 0 0
0.1
0.2
0.3
0.4 r [nm]
0.5
0.6
0.7
Figure 4: Comparison of time-averaged radial distribution functions calculated for the wear particles, the shear zone, and the bulk of the system with cubic abrasives, 1.0 GPa and 16 m/s.
4. Contact zone The smooth particle based estimation of the real contact area, as proposed by Eder et al. and discussed or applied in Refs. [24, 25, 26, 27], is quite time consuming for larger systems and currently reliable for singleasperity contacts only. We therefore quantify the contact zone via the number of atoms of one counterbody currently in contact with the other one, which according to Refs. [24, 25] is reasonably proportional to the contact area, irrespective of load and asperity shape. Here the problem of the definition of the counterbodies arises. In this work, we define the atoms of the abrasives together with the wear particles as the upper counterbodies. The lower counterbody consists of all other atoms which lie no more than 0.3 nm, the 7
4
4
x 10
2.5
2
2
1.5
1.5
Nwp [−]
Nwp [−]
2.5
1
1
0.5
0.5
0 0
20
40 ∆x [nm]
60
x 10
0 0
80
0.2
(a)
16000
16000
14000
14000
12000
12000
10000 8000
4000
4000
2000
2000 40 ∆x [nm]
0.8
1
0.8
1
8000 6000
20
0.6 σz [GPa]
10000
6000
0 0
0.4
(b)
Nshear [−]
Nshear [−]
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65
60
80
(c)
0 0
0.2
0.4
0.6 σz [GPa]
(d)
Figure 5: Sliding distance dependence (left) and averaged pressure dependence (right) of the number of atoms in wear particles (a–b) and in the shear zone (c–d). Green hues denote two-body contact sliding with cubic abrasives, red ones denote two-body contact sliding with spherical abrasives, and blue ones three-body contact sliding. The intensity of the color reflects the normal load. Dashed lines and small symbols represent sliding at 8 m/s, solid lines and large symbols at 16 m/s.
8
maximum contact distance, lower than the lowest atom of the upper counterbodies. These groups of atoms are well-defined at all times where there exists a filtered velocity profile. For computational efficiency, we have reduced the large number of atoms in the abrasives by considering only their bounding atoms, identified via their centro-symmetry parameter [28] which is larger than 5 in the initial configuration (before contact occurs). We then search for any atoms in the lower counterbody which lie within a distance dcont of an abrasive or a wear particle. This maximum contact distance dcont is chosen in such a way that the contact criterion includes all first (dcont = 0.265 nm), or first and second neighbors (dcont = 0.3 nm), cf. the radial distribution function in Fig. 4. The differences in the number of contact atoms between these two variants are small and nearly constant for constant load, although they do vary with load. In the two-body contact sliding simulations, almost all contact atoms lie within the shear zone, while in the three-body contact sliding simulations, quite a large proportion of the contact atoms are not sheared. This illustrates how free rolling abrasives hardly impose any shear on the counterbody, especially at low loads. Fig. 6 shows a comparison of the number of contact atoms with the number of shear zone atoms for three-body contact sliding at low load and two-body contact sliding at high load. Note that during three-body contact sliding, the number of contact atoms exceeds the number of shear zone atoms by a large margin. Some exemplary contact zone maps can be seen in Fig. 7, where the coloring according to height value gives an idea of the individual contact zone geometries. In Fig. 8 we compare the sliding distance and the pressure dependence of the total number of contact atoms (found for dcont = 0.265 nm) between the three considered configurations. The number of contact atoms initially decreases with the sliding distance, especially for medium and high pressures, but reaches a steady state towards the end of the simulations, and the pressure dependence is almost linear.
incl. 1st neighbors incl. 2nd neighbors all shear zone atoms
16000 14000
16000 14000 12000 Ncontact [−]
12000 Ncontact [−]
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65
10000 8000
10000 8000
6000
6000
4000
4000
2000
2000
0 0
20
40 ∆x [nm]
60
80
(a)
0 0
20
40 ∆x [nm]
60
80
(b)
Figure 6: The number of contact atoms compared with the shear zone atoms at 16 m/s. (a) Three-body contact sliding at 0.1 GPa, (b) two-body contact sliding with cubic abrasives at 1.0 GPa. By tracking which contact atom is in contact with which abrasive particle, it is possible to clusterize the contributions of the individual abrasives to the total contact zone. However, this only works while the abrasives are in direct contact with the lower counterbodies, which is the predominant case in three-body contact sliding simulations. During two-body contact sliding, the formation of wear particles soon shifts the main part of the contact area towards contact between wear particle and lower counterbody, see Fig. 9 (a). In this case, the direct affiliation of a contact atom touching a wear particle with one of the abrasives is not as straightforward as for three-body contact sliding. This can be overcome by iteratively searching for contact atoms which lie closer than dclust = 0.4 nm to other contact atoms whose abrasive affiliation is already known. The distance dclust can be chosen somewhat larger than dcont to reduce the number of required 9
(a)
(b)
Figure 7: Contact zone maps after '4 ns of sliding at 16 m/s. (a) Three-body contact sliding at 0.1 GPa, (b) two-body contact sliding with cubic abrasives at 1.0 GPa. Coloring is rainbow style according to height value of the contact atom, red is high, blue is low.
5000
5000
4000
4000 Ncontact [−]
Ncontact [−]
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65
3000
3000
2000
2000
1000
1000
0 0
20
40 ∆x [nm]
60
80
(a)
0 0
0.2
0.4
0.6 σz [GPa]
0.8
1
(b)
Figure 8: Sliding distance (a) and averaged pressure dependence (b) of the number of contact atoms. Color scheme as in Fig. 5.
10
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65
iterations, but at some point there is a risk of affiliations incorrectly spreading to neighboring clusters. The number of affiliated contact atoms may stop growing due to problems arising from the periodic boundary conditions, namely unaffiliated atoms with their affiliated neighbors at the other side of the simulation cell. This can be solved by once mapping all atoms in question which lie closer than dclust to the lower cell boundary across the upper cell boundary. Thus for the next iteration, a sufficient number of affiliated atoms is available on the correct side of the cell boundary so that more affiliations can be found. Finally, single outliers and small groups of outlying atoms can join the closest clusters with known affiliation. This last step is the only appreciable source of error, as isolated contact zones without a single initial particle affiliation will thus be falsely affiliated with neighboring abrasives.
(a)
(b)
(c)
Figure 9: Illustration of the clustering algorithm for the contact zone. Coloring is according to abrasive particle affiliation, and dark red is reserved for unaffiliated atoms in contact with wear particles. (a) Only atoms in direct contact with the numbered abrasive particles have known affiliations. (b) After clustering without taking the kinematics of the simulation into account, note the incorrect “backwards growth” of the light blue contact zone on the left. (c) Clustering corrected by taking simulation kinematics into account. Also observe that the light green area found at the top left due to periodic boundary conditions is correctly affiliated with the cluster at the top right. The above method is limited to non-coalescent contact clusters. As soon as the build-up of wear debris between two abrasives becomes sufficient to allow a single uninterrupted contact zone, it becomes a matter of definition which part of the contact should be identified with a particular abrasive. In this work, a wear particle related contact zone is attributed to the abrasive which formed the corresponding wear particle. In this case, the algorithm initially searches for affiliated atoms only in sliding direction, which effectively hinders the iterative “backwards growth” which would otherwise occur in coalescing contact zones, see Figs. 9 (b) and (c). The time-resolved knowledge of the size and the geometry of the individual contact zones allows an assessment of the influence of the process parameters on the contact situation. Fig. 10 shows that although the total number of contact atoms is very similar at constant normal load and sliding velocity, the distribution of the load and its time development varies considerably between the three considered surface finishing processes. A possible application of this clustering method is a time-resolved statistical analysis of abrasive indentation depths and orientations. Due to variations in surface topography and wear particle geometry, these parameters keep changing for every contact zone, and it is well known that they have considerable influence on real wear processes. These aspects as well as a sub-surface stress analysis are beyond the scope of this work. 11
3000
3500
2500
2500
3000
2000 1500 1000 500 0 0
Ncontact (per−particle) [−]
3000
Ncontact (per−particle) [−]
Ncontact (per−particle) [−]
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65
2000 1500 1000 500
20
40 ∆x [nm]
60
80
0 0
(a)
2500 2000 1500 1000 500
20
40 ∆x [nm]
60
80
(b)
0 0
20
40 ∆x [nm]
60
80
(c)
Figure 10: Break-down of the number of contact atoms by abrasive particle over sliding distance. Coloring is rainbow style with dark blue representing abrasive particle number 1 and red representing abrasive particle number 16. The three shown examples each represent one of the finishing methods at 0.5 GPa and 16 m/s. Left: three-body contact sliding, center: two-body contact sliding with spherical abrasives, right: two-body contact sliding with cubic abrasives. 5. Topography and wear evaluation For an efficient quantitative evaluation of the surface topography over time, the point cloud of atoms located on the surface must be determined and mapped to a regular grid. To achieve this, the simulation box was divided into square cells of side length dcell = 0.3966 nm, which is a trade-off between good lateral resolution and certainty of finding enough atoms in each cell. For every time step the substrate atom with the highest z value in each cell was then identified, i.e., all wear particle and shear zone atoms were excluded from this consideration. This procedure yields a time-resolved 2D representation zsubst (x, y; t) of the surface topography on a regular mesh of 72×72 points. Averaging along the x and y directions gives the mean substrate height for every time step, while the standard deviation of zsubst (x, y; t = const.) gives the Sq (t) roughness parameter at that time, both of which are shown as functions of sliding distance and normal load in Fig. 11. By monitoring the development of Sq (t), one can easily see the wear progress as the surface topography is flattened. Forming the pointwise difference ∆z(x, y; t) = zsubst (x, y; t0 ) − zsubst (x, y; t)
(1)
yields the wear height, which is positive where asperities are reduced in height, while it is negative where pits are filled up with material. This can be presented in the form of a wear map as shown in Fig. 12, where the former regions are shown in red and the latter ones in blue. The positive and negative contributions can be summed up individually and multiplied with the xy area of the mesh elements to yield the volumes of (−) (+) (−) (+) asperity height reduction Vasp and pit fill-up Vpit , see Fig. 13. Here we can see that Vasp and Vpit both saturate over the course of the simulation depending on the pressure. Note that the pressure dependence of (+) Vpit has a weak maximum at σz = 0.5 GPa for the two-body contact sliding simulations. The smoothing of the atomic velocity profiles, which determines the categorization of the atoms as substrate, shear or wear atoms, limits the number of time steps available for this analysis, so half the smoothing window width is lost at the beginning and the end of each dataset. The wear volume is estimated on the basis of the substrate atoms only, and there is usually an imbalance (−) (+) between Vasp and Vpit . The difference between the two, ∆V , can be decomposed into contributions from the atoms which form wear particles stuck to the abrasives, and those currently under shear between the wear particles and the substrate. As these two atomic categories are known by number, not by volume, an appropriate per-atom volume has to be found as a conversion factor. This per-atom volume should be in the range of the equilibrium volume of bcc Fe of 0.0116 nm3 . 12
3.9
3.8
3.8
zsubst [nm] (final)
zsubst [nm]
3.9
3.7 3.6
3.4 0
3.7 3.6 3.5
3.5
20
40 ∆x [nm]
60
3.4 0
80
0.2
0.4
(a)
0.7
0.6
0.6 Sq [nm] (final)
0.7
0.5 0.4
0.2
0.2 40 ∆x [nm]
1
0.8
1
0.4 0.3
20
0.8
0.5
0.3
0
0.6 σz [GPa] (b)
q
S [nm]
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65
60
80
(c)
0
0.2
0.4
0.6 σz [GPa] (d)
Figure 11: Sliding distance (left) and pressure dependence (right) of the mean surface height (a–b) and the Sq roughness parameter (c–d). Color scheme as in Fig. 5.
13
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65
(a)
(b)
Figure 12: Wear map after 4.8 ns of three-body contact sliding at 0.1 GPa and 16 m/s (a) and two-body contact sliding with cubic abrasives at 1.0 GPa and 16 m/s (b). Red areas denote regions of maximum height reduction (∆z (max) = +2.5 nm), blue areas denote regions of maximum fill-up (∆z (min) = −2.5 nm), while in the green areas there was no net change of height. However, in some cases, especially in three-body contact sliding simulations, ∆V cannot be neatly decomposed into these two contributions. There exists a small third, monotonically increasing and saturating contribution to the volumetric gap, which seems to arise from dynamic compression of the sub-surface substrate layer. To quantify this, we count the number of atoms in a 0.5 nm thick substrate layer immediately below the atom which defines zsubst (x, y; t), giving Nsubsurf (x, y; t) on the regular mesh introduced for the topography evaluation. As this layer has a constant volume over time, we can thus eliminate the error introduced by having to calculate the P total volume of the substrate. The total volume of the sub-surface layer should initially contain some x,y Nsubsurf (x, y; t0 ) ' 35000 atoms. However, the number actually obtained is closer to 34400 due to the influence of the topography and the somewhat crude resolution of the mesh. It may therefore be assumed that the exactness of this estimation method for the sub-surface compression is of the order ±3%. By subtracting the initial value from the current one, we obtain the number of surplus atoms in the sub-surface layer, X X Nsubsurf (x, y; t) − Nsubsurf (x, y; t0 ). (2) ∆Nsubsurf (t) = x,y
x,y
The fact that this sub-surface compression plays a non-negligible role predominantly during three-body contact sliding, where it can constitute the major contribution to ∆V , corresponds well to the wear regimes discussed in Ref. [15], where free rolling abrasives primarily impose pressure on the substrate. The sliding distance and pressure dependence of ∆Nsubsurf (t) is shown in Fig. 14. Here, we see evidence that the sharp-edged cubic abrasive particles lead to cutting wear, predominantly acting in the lateral direction and producing little normal compression, while the substrates underneath the spherical abrasive particles are compressed more. Now that the number of atoms are known which constitute the wear particles, Nwp (t), those that are in the shear zone, Nshear (t), and those additionally compressed into the sub-surface layer, ∆Nsubsurf (t) for every time step, we can fit per-atom volumes Va1 and Va2 as coefficients to these numbers to yield the difference between the total asperity volume reduction and the total pit fill-up volume such that ∆V (t) = Va1 [Nwp (t) + Nshear (t)] + Va2 ∆Nsubsurf (t).
(3)
All values obtained for Va1 lie close to the expected bulk value regardless of load, while the values of Va2 can vary substantially for small loads. However, a comparison with the time development of the per-atom 14
500
400
400 V(−) [nm3] asp
V(−) [nm3] asp
500
300
300
200
200
100
100
0 0
20
40 ∆x [nm]
60
0 0
80
0.2
250
250
200
200
150
100
0.8
1
0.8
1
0.8
1
150
100
20
40 ∆x [nm]
60
0 0
80
0.2
(c)
400
400
350
350
300
300 ∆V [nm3]
450
250 200
200 150
100
100
50
50 40 ∆x [nm]
0.6 σz [GPa]
250
150
20
0.4
(d)
450
0 0
0.6 σz [GPa]
50
50
0 0
0.4
(b)
V(+) [nm3] pit
V(+) [nm3] pit
(a)
∆V [nm3]
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65
60
80
(e)
0 0
0.2
0.4
0.6 σz [GPa] (f)
Figure 13: Sliding distance (left) and averaged pressure dependence (right) of the asperity height reduction volume (a–b), the pit fill-up volume (c–d) and the difference between the two (e–f). Color scheme as in Fig. 5. 15
3000
3000
2500
2500
2000
2000 ∆Nsubsurf [−]
∆Nsubsurf [−]
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65
1500 1000 500
1500 1000 500
0
0
−500
−500
−1000 0
20
40 ∆x [nm]
60
80
−1000 0
0.2
0.4
0.6 σ [GPa]
0.8
1
z
(a)
(b)
Figure 14: Sliding distance (a) and averaged pressure dependence (b) of the surplus atoms in the 0.5 nm thick sub-surface substrate layer compared with the initial state. Color scheme as in Fig. 5. volume for the entire substrate shows that overall number density changes are within a range of ±0.5% for all systems. With the knowledge of the coefficients Va1 and Va2 , the asperity volume reduction can now be decomposed into volumetric contributions of matter filled in the pits, wear particles stuck to the abrasives, matter currently under shear, and regions with sub-surface compression, see Fig. 15. This figure compares three-body contact sliding at low pressure with two-body contact sliding at high pressure with cubic abrasives. It is evident that the difference between the asperity reduction and the pit fill-up volumes varies considerably between these two processes. Virtually no wear particles are generated during three-body contact sliding, and ∆V consists of contributions from sub-surface compression and the shear zone. During two-body contact sliding, about half of the asperity reduction volume consists of wear particles, and the rest is evenly divided between material filled into the pits and the shear zone, while sub-surface compression is virtually negligible. By adapting the clustering algorithm described above to the identification of any wear particle with the abrasive that caused it, we can further break down the wear particles into contributions from the individual abrasives. The three examples shown in Fig. 16 give an impression of how this analysis method can give an overview over the distribution of wear volume formation. In the three-body contact sliding simulation, shown in Fig. 16 (a), wear particles are essentially produced by two abrasives only, while two-body contact sliding with cubic abrasives, shown in Fig. 16 (c), leads to a very even distribution of wear particles. Here, 12 out of 16 abrasives produce wear particles of similar volume which remain quite stable until the end of the simulation. In the data obtained for two-body contact sliding with spherical abrasives, shown in Fig. 16 (b), we can see a rather irregular distribution of the wear particles among the abrasives as well as evidence of the formation of transient wear particles, i.e., wear particles that are deposited in pits at some point and thus become part of the substrate again (e.g., the light blue area in the center portion). 6. Summary and Conclusion We have presented a method which allows the quantitative evaluation of the wear volume and the real contact zone resulting from two-body and three-body wear of nanoscopically rough surfaces using MD simulations. The focus of this work lay on implementing a multi-abrasive approach to this problem which enables us to study the effects of multiple abrasives in close proximity and changing surface topography. 16
500 400 3
400
asperity reduction pit fill−up wear particles shear zone sub−surface compression Vwear [nm ]
500
Vwear [nm3]
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65
300 200 100 0 0
300 200 100
20
40 ∆x [nm]
60
80
(a)
0 0
20
40 ∆x [nm]
60
80
(b)
Figure 15: Sliding distance dependence of the volumetric contributions to the total asperity reduction volume (shown in blue) at 16 m/s. (a) Three-body contact sliding with spherical abrasive particles, 0.1 GPa. (b) Two-body contact sliding with cubic abrasive particles, 1.0 GPa. In particular, we outlined a scheme to generate an initial surface topography with a pseudo-random Gaussian roughness distribution which can be easily configured and scaled while being periodically replicable. We discussed which constraints on the individual abrasive particles mimic certain surface finishing process conditions. An analysis of the atomic advection velocities allowed the dynamic identification of wear particles, the substrate and the shear zone in-between. We then presented a method for the time-resolved quantification of the real contact zone between the abrasives and the wear particles on the top, and the shear zone and substrate on the bottom. Using a clustering technique which takes into account the kinematics of the process, we could break down the multi-asperity contact zone into its single-asperity contributions. Using a mesh-based topography evaluation, we could generate time-dependent wear maps and break down the asperity reduction volume into its contributions, including the pit fill-up volume, the wear particles, the shear zone, and volume lost through substrate compression. By applying the clustering technique to the wear particles themselves, it was possible to analyze the contributions of the individual abrasives to wear particle formation. Acknowledgements This work was funded by the Austrian COMET-Program (Project K2 XTribology, Grant No. 824187), the ERDF as well as the province of Lower Austria (Onlab project), and it was carried out at the “Excellence Centre of Tribology”. S. J. Eder would like to thank Georg Vorlaufer for the helpful discussion about the sub-surface compression of the substrate. References [1] I. Szlufarska, M. Chandross, R. W. Carpick, Journal of Physics D: Applied Physics, 41 (2008) 123001. [2] E. Brinksmeier, J.C. Aurich, E. Govekar, C. Heinzel, H.-W. Hoffmeister, F. Klocke, J. Peters, R. Rentsch, D.J. Stephenson, E. Uhlmann, K. Weinert, M. Wittmann, CIRP Annals-Manufacturing Technology, 55 (2006) 667. [3] R. Rentsch, I. Inasaki, CIRP Annals-Manufacturing Technology, 43 (1994) 327. [4] K. Maekawa, A. Itoh, Wear, 188 (1995) 115. [5] L. Zhang, H. Tanaka, Wear, 211 (1997) 44. [6] R. Komanduri, N. Chandrasekaran, L. M. Raff, Wear, 240 (2000) 113. [7] Q. X. Pei, C. Lu, H. P. Lee, Y. W. Zhang, Nanoscale research letters, 4 (2009) 444.
17
3500
14000
600
3000
12000
500 400 300 200 100 0 0
Nwp (per−particle) [−]
700
Nwp (per−particle) [−]
Nwp (per−particle) [−]
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65
2500 2000 1500 1000 500
20
40 ∆x [nm]
60
80
0 0
10000 8000 6000 4000 2000
20
40 ∆x [nm]
60
80
0 0
20
40 ∆x [nm]
(a)
(b)
(c)
(d)
(e)
(f)
60
80
Figure 16: Top row: break-down of the number of wear particle atoms Nwp (t) by wear particle over sliding distance. Bottom row: top view of the wear particle distribution at ∆x = 40 nm (t = 2.5 ns). Coloring is rainbow style with dark blue representing wear particle number 1 and red representing wear particle number 16. The three shown examples each represent one of the finishing methods at 0.5 GPa and 16 m/s. left: three-body contact sliding, center: two-body contact sliding with spherical abrasives, right: two-body contact sliding with cubic abrasives.
18
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65
[8] [9] [10] [11] [12] [13] [14] [15] [16] [17] [18] [19] [20] [21] [22] [23] [24] [25] [26] [27] [28]
L. Si, D. Guo, J. Luo, X. Lu, Journal of Applied Physics, 107 (2010) 064310. P. M. Agrawal, L. M. Raff, S. Bukkapatnam, R. Komanduri, Applied Physics A, 100 (2010) 89. J. Zhong, R. Shakiba, J. B. Adams, Journal of Physics D: Applied Physics, 46 (2013) 055307. P. Spijker, G. Anciaux, J.-F. Molinari, Tribology International, 59 (2013) 222. X. Han, Y. Hu, S. Yu, Applied Physics A, 95 (2009) 899. R. Chen, R. Jiang, H. Lei, M. Liang, Applied Surface Science, 264 (2012) 148. L. Zhang, H. Zhao, Z. Ma, H. Huang, C. Shi, W. Zhang, AIP Advances, 2 (2012) 042116. L. Zhang, H. Tanaka, Tribology International, 31 (1998) 425. J. Sun, L. Fang, J. Han, Y. Han, H. Chen, K. Sun, Computational Materials Science, 82 (2014) 140. P. M. Agrawal, R. Narulkar, S. Bukkapatnam, L. M. Raff, R. Komanduri, Tribology International, 43 (2010) 100. M. I. Mendelev, S. Han, D. J. Srolovitz, G. J. Ackland, D. Y. Sun, M. Asta, Philosophical Magazine, 83 (2003) 3977. T.-H. Fang, C.-I. Weng, Nanotechnology, 11 (2000) 148. Q. X. Pei, C. Lu, F. Z. Fang, H. Wu, Computational materials science, 37 (2006) 434. ISO 25178-2:2012 Geometrical product specifications (GPS) – Surface texture: Areal – Part 2: Terms, definitions and surface texture parameters. August 2012. J. R. D’Errico, Understanding gridfit, 2006. URL: http://www.mathworks.at/matlabcentral/fileexchange/8998-surfacefitting-using-gridfit (Last visited: Feb. 7, 2014). S. J. Plimpton, J. Comput. Phys., 117 (1995) 1. S. Eder, A. Vernes, G. Vorlaufer, G. Betz, J. Phys.: Condens. Matter, 23 (2011) 175004. S. Eder, A. Vernes, G. Betz, Computer Physics Communications, 185 (2014) 217. S. J. Eder, A. Vernes, G. Betz, Langmuir, 29 (2013) 13760. A. Vernes, S. Eder, G. Vorlaufer, G. Betz, Faraday Discuss., 156 (2012) 173. C. L. Kelchner, S. J. Plimpton, J. C. Hamilton, Phys. Rev. B, 58 (1998) 11085.
19
S0010-4655(14)00183-0 http://dx.doi.org/10.1016/j.cpc.2014.05.018 COMPHY 5345
To appear in:
Computer Physics Communications
Received date: 7 February 2014 Revised date: 16 April 2014 Accepted date: 23 May 2014 Please cite this article as: S.J. Eder, D. Bianchi, U. Cihak-Bayr, A. Vernes, G. Betz, An analysis method for atomistic abrasion simulations featuring rough surfaces and multiple abrasive particles, Computer Physics Communications (2014), http://dx.doi.org/10.1016/j.cpc.2014.05.018 This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65
An analysis method for atomistic abrasion simulations featuring rough surfaces and multiple abrasive particles S. J. Edera,∗, D. Bianchia , U. Cihak-Bayra , A. Vernesa,b , G. Betzb a Austrian b Institute
Center of Competence for Tribology, Viktor-Kaplan-Straße 2, 2700 Wiener Neustadt, Austria of Applied Physics, Vienna University of Technology, Wiedner Hauptstraße 8–10 /134, 1040 Vienna, Austria
Abstract We present a molecular dynamics (MD) model system to quantitatively study nanoscopic wear of rough surfaces under two-body and three-body contact conditions with multiple abrasive particles. We describe how to generate a surface with a pseudo-random Gaussian topography which is periodically replicable, and we discuss the constraints on the abrasive particles that lead to certain wear conditions. We propose a postprocessing scheme which, based on advection velocity, dynamically identifies the atoms in the simulation as either part of a wear particle, the substrate, or the sheared zone in-between. This scheme is then justified from a crystallographic order point of view. We apply a distance-based contact zone identification scheme and outline a clustering algorithm which can associate each contact atom with the abrasive particle causing the respective contact zone. Finally, we show how the knowledge of each atom’s zone affiliation and a timeresolved evaluation of the substrate topography leads to a break-down of the asperity volume reduction into its components: the pit fill-up volume, the individual wear particles, the shear zone, and the sub-surface substrate compression. As an example, we analyze the time and pressure dependence of the wear volume contributions for two-body and three-body wear processes of a rough iron surface with rigid spherical and cubic abrasive particles. Keywords: atomistic abrasive wear, molecular dynamics, multiple particle abrasion PACS numbers: 46.55.+d, 31.15.xv, 68.35.Ct
1. Introduction The increasing miniaturization of electronic and mechanical devices requires certain surfaces to be finished to a high degree of smoothness. Even some macroscopic surfaces, such as those of molding parts for the optical industry, must be virtually scratch-free in order to produce high-quality optical products. As scratch dimensions enter the (sub-)nanometric length scale, experimental treatment and analysis of the workpieces becomes complicated, expensive, or even impossible. Similarly, with manufacturing tolerances becoming ever smaller, the wear mechanisms leading to geometrical changes in components have to be understood on the atomistic scale. Here molecular dynamics (MD) simulation can reveal important insights into the processes of interest [1]. In literature, MD simulations of two-body and three-body contact sliding are discussed either in terms of tribological phenomena such as abrasive wear, or as surface finishing processes such as grinding and polishing. An extensive overview of the simulation of the grinding process at several length scales can be found in Ref. [2], where it is noted that compared to nano-cutting, little MD literature is available on grinding. The first earnest attempts to explicitly simulate two-body and three-body wear at an atomistic level were undertaken in the mid-1990s [3, 4, 5]. Most MD simulations treat grinding as a nanometric ∗ Corresponding
author Email address: eder at ac2t.at (S. J. Eder)
Preprint submitted to Elsevier
April 16, 2014
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65
cutting process and focus on the interaction between one single abrasive particle and a usually flat substrate [6, 7, 8]. The atomistic wear of even a single roughness feature due to asperity-asperity interaction has only very recently received some attention, e.g., in Refs. [9, 10]. While these approaches can lead to an understanding of a process under controlled conditions, e.g., at constant penetration depth and orientation of the abrasive particle, both wear and finishing processes may involve circumstances such as changing rough topographies, random rake angle distributions or wear debris build-up between neighboring abrasives, which cannot be studied by considering only a single abrasive interacting with a flat surface. For these reasons, the implementation of a statistically based surface roughness, similar to what was recently used in sliding friction simulations [11], combined with multiple abrasive particles acting on the surface, is highly desirable. Literature on the atomistic modeling and simulation of three-body wear is scarce if one is not primarily interested in the chemo-mechanical planarization of silicon [8, 12, 13, 14]. Especially the extended degrees of freedom (rolling motion) of the particles have been considered in very few publications [15, 16]. An MD approach to the polishing of a silicon surface with several randomly distributed asperities and abrasive particles was proposed in Ref. [17], leading to a phenomenologically based material removal rate. In this work, we present an easily scalable and extendable MD model system which can be used to quantify the wear volume and the contact zone occurring during wear processes on a rough surface with multiple hard abrasive particles. The model systems were intentionally kept simple and do not include a medium or phenomena such as the wear of the abrasive particles. However, our approach is general enough to allow its application to the analysis of material removal rates in processes such as polishing and grinding as well. Note that a comprehensive discussion, from a materials science point of view, of the results obtained for our chosen model systems is beyond the scope of this publication. In Sec. 2, we discuss the model setup, the procedure for pseudo-random Gaussian topography generation as well as the constraints which lead to two-body or three-body contact sliding. Sec. 3 deals with the identification of wear particles based on their advection velocities. In Sec. 4 we define and quantify the contact zone between the abrasives plus the wear particles and the remaining substrate. Here we also discuss a clustering algorithm which allows the affiliation of every contact zone with an abrasive particle. Finally, in Sec. 5, we perform a mesh-based volumetric substrate topography analysis on which we base the wear volume estimation. The knowledge of the wear particles allows a break-down of the total asperity reduction volume into several contributions. By applying a modified version of the clustering algorithm, the wear particle volume is then broken down further into contributions due to single abrasive particles. 2. Model setup and simulation procedure The substrate consists of a block of bcc Fe with lateral dimensions of 28.5×28.5 nm2 and sufficient initial height to accommodate the pseudo-random Gaussian roughness distribution, the generation procedure of which will be discussed in more detail later. At its base, three monolayers of atoms form a rigid support, covered by two monolayers of Fe thermostatted to 300 K with a Langevin thermostat with a time constant of 0.5 ps, acting only in y direction in order not to overly interfere with the load acting in z direction and the movement of the abrasive particles, which is predominantly in x direction. Periodic boundary conditions are applied in the x and y directions, which means that if an atom leaves the simulation box through one boundary, it re-enters at the opposing one. No matter is lost during a simulation run. The abrasive particles are modeled as rigid bodies arranged in a staggered 4×4 grid, rotated randomly between 0 and 90 degrees about all three Cartesian axes, and are initially placed several nanometers above the surface. Two abrasive particle geometries are considered: spherical with a radius of 3.25 nm and cubic with a side length of 4.2 nm. While these extreme geometries are unlikely to occur in real processes, they were chosen for illustrative purposes due to their likelihood of producing clear differences in the results. For similar reasons, the areal density of the abrasive particles was substantially increased compared to realistic systems in order to speed up the generation rate of wear debris and thus reduce computational expense. Snapshots of the initial configurations are shown in Fig. 1. A Finnis-Sinclair potential is used for the Fe–Fe interactions [18]. The interactions between abrasive particles and the Fe surface is modeled using a Lennard-Jones potential with ε = 0.125 eV and σ = 2
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65
0.2203 nm, a parameter set comparable to Refs. [4, 5, 7, 19, 20]. As the abrasive spheres in the threebody contact simulations have the possibility to touch each other, we assumed a very weakly attractive Lennard-Jones potential for this interaction with the parameters ε = 0.03 eV and σ = 0.2203 nm. We built a periodically replicated surface topography with Gaussian height distribution. The parameters which influence the final topography are the roughness parameter Sq , which is the standard deviation from (Asp) the mean surface height zsubst [21] as well as the nominal asperity width dx/y . First, these nominal values are replaced with ones of similar size so that the simulation box may be filled in x and y direction with an integer number of the resulting cells. The cell centers constitute the control points of the topography and are assigned a random height (z) value such that the total distribution is as close to Gaussian as can be achieved with the relatively small number of control points. The resulting 2D-matrix of z values is then padded at each boundary with two rows (or columns) from the opposite side to account for periodicity. The respective x and y values of the support points are offset by random values between 0 and a maximum (offset) (Asp) offset dx/y < dx/y , and then also padded to produce a periodic boundary. For the initial topography used as an example throughout this work, we assigned the parameter values zsubst = 4 nm, Sq = 0.8 nm, (Asp) (offset) dx/y = 2 nm and dx/y = 1 nm. On a realistic surface, the typical height of a roughness feature would be approximately 2–3 orders of magnitude smaller than its lateral dimensions. As we are primarily interested in (Asp) the illustration of our analysis method, in this example Sq and dx/y were chosen to be of similar proportion, leading to considerably steeper asperity flanks. The scattered point cloud generated in this fashion constitutes a very coarse representation of the topography. In order to filter an initial set of substrate atoms according to this surface, it is necessary to know the z value at any given lateral position. To achieve this, the scattered data was smoothed and mapped to a fine mesh (of the order of the lattice parameter of Fe) using a modified ridge estimator algorithm [22]. The application of the resulting highly resolved topography as a criterion for surface generation yields a perfectly periodic atomistic topography with smooth transitions between roughness features. It should be noted that the originally set roughness of the surface, characterized by Sq , is slightly reduced by this smoothing. This must be taken into account if a very exact initial roughness is required. A brief relaxation of the substrate to a temperature of 300 K with repeated re-assignment of atomic velocities according to the Maxwell-Boltzmann velocity distribution leads to the initial state of the substrate used throughout this work which is shown in Fig. 1. When simulating two-body contact sliding, the abrasive particles are dragged across the surface in +x (max) (max) (max) (max) /8, resulting in a sliding velocity of v (max) ' vx = vx direction at vx and in +y direction at vy ◦ at an angle of 7.125 with the x-axis. This way the abrasives re-enter the simulation box at different y positions every time they pass through the periodic box boundary, thus scanning the surface more evenly and avoiding artificially amplified deep scratches. The relative positions and orientations of the abrasives are kept constant as would be the case for hard, rough counterbody such as a grindstone. The normal load is kept constant over time and distributed evenly between the 16 particles, which can change their z position collectively as required by the topography, but not individually. This may lead to loss of surface contact for individual abrasive particles for some time intervals, while the remaining particles share the additional load. In the simulation of three-body contact sliding, the abrasive particles are moved across the surface (max) parallel to the +x direction at vx . In contrast to the two-body contact simulations, here they may rotate about their respective centers of mass according to the atomic interactions with the substrate surface and the neighboring abrasive particles. The normal load is distributed evenly between the 16 particles, and each particle may individually change its y and z position as required according to the normal load, the sliding velocity, and the local topography, as would be the case for dust or debris in the lubrication gap of a tribosystem or in a polishing slurry. These additional degrees of freedom lead to a contact situation quite different from two-body contact sliding, since the number of load bearing regions is always equal to the number of abrasive particles, and the rolling motion of the abrasives leads to considerably less shear. We performed simulations of three-body wear processes with spherical abrasives and two-body wear processes with spherical or cubic abrasives at pressures of σz = 0.1, 0.5 and 1.0 GPa, and at sliding velocities of 8 and 16 m/s using the open source MD code LAMMPS [23]. The total simulation time was 3
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65
(a)
(b)
(c)
Figure 1: Initial simulation configurations. Abrasives are shown in green, the substrate up to the mean surface height in blue, and initial asperities in yellow. (a) 16 randomly rotated rigid cubic abrasive particles for two-body contact sliding simulations. (b) 16 rigid spherical abrasive particles for two-body and threebody contact sliding simulations. (c) Top view of initial surface topography, colored according to atomic z coordinate, where the peaks are red and the pits are blue. The abrasive particles are not shown.
4
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65
(max)
(max)
5 ns for v1 = 16 m/s and 10 ns for v2 of ∆x = 80 nm.
= 8 m/s to yield a total translation distance along the x axis
3. Identification of wear particles Fig. 2 shows top views of the system configurations without the abrasives after 5 ns of sliding at v (max) = 16 m/s, colored according to the atomic z coordinate. Although a comparison between the system in Fig. 2 (a), which was obtained with three-body contact sliding at low load, and the one in Fig. 2 (i), which is the result of two-body contact sliding with cubic abrasives at high load, clearly shows the formation of wear particles in the latter case, the automatic and time-resolved identification of the atoms belonging to this group is not so straightforward. The distinct horizontal grooves in Figs. 2 (b) and (c) occur due to interlocking of neighboring abrasives, leading to effective two-body contact sliding of one of the particles as its rolling motion is inhibited. The grooves are perfectly horizontal because in the three-body contact sliding simulations, the only kinematic constraint on the abrasives is the motion in x direction while the motion in y direction results mainly from the current local topography. If atoms were removed from the surface in such a fashion that there exists no Fe connection between them and the substrate, this detachment could serve as an indicator for wear debris. However, since this applies only to very few atoms, which cannot form a real wear particle, our proposed definition of wear particles is based on the advection velocity of the atoms. This method is akin to the contact atom counting procedure described in further detail in previous work by Eder et al. in Refs. [24] and [25], which is used to identify the contacting Fe atoms between two Fe sliders. Here, the main component of the advection velocity, v = vx , is calculated via the position differences ∆x evaluated at a sampling interval of ∆t = 20 ps. This signal is cleaned from the artifacts of the periodic boundary conditions and filtered using a Gauss filter with a total width of 0.8 ns and a standard deviation of 0.4 ns (for v (max) = 8 m/s) to yield continuous velocity profiles for each atom. Based on the instantaneous advection velocity thus calculated, each atom may be assigned to one of three groups: (a) atoms traveling at v > 0.9 v (max) are considered stuck to an abrasive and are therefore classified as wear particles, (b) atoms traveling at |v| < 0.1 v (max) are considered stationary and are classified as substrate atoms, and (c) all remaining atoms, which move faster than 0.1 v (max) (not necessarily with the same sign, e.g., in three-body contact sliding simulations), are classified as part of what we will henceforth refer to as the “shear zone”. Fig. 3 gives a simplified overview of the discussed zones. To justify this seemingly arbitrary atom identification scheme not only from an angle of computational ease, but also from a crystallographic one, we calculated and compared the time-averaged radial distribution functions (RDFs) of the wear particles, the shear zones, and the average bulk substrate, see Fig. 4. For the two-body contact sliding simulations, a comparison between the RDF for the average bulk, which exhibits the typical pattern of a bcc lattice with broadened peaks due to finite temperature, and that of the shear zone as identified using the above scheme, proves that the latter region has a considerably disturbed lattice. The coalescence of the first with the second and the third with the fourth peak suggest that this region is almost amorphous. The wear particles, on the other side, seem to crystallize onto the abrasives in such a way that they show almost the same pattern as the bulk substrate, the main difference being a height reduction of the second-neighbor peak. In the three-body contact sliding simulations, the shear zone is crystallographically much more similar to the wear particles. For both regions, crystallinity is generally retained, with only slight pairwise coalescence of neighboring peaks. However, these regions play a much smaller quantitative role in the three-body contact sliding simulations, since the rolling motion of the abrasives results in much less shear, and the non-cutting wear regime leads to few wear particles. We therefore consider our atom identification scheme sufficiently universal for application to all our simulations. Fig. 5 compares the sliding distance and the pressure dependence of the total number of atoms in wear particles and the shear zone for all considered wear processes. All pressure dependences are nearly linear. We can observe that the sliding velocity has little impact on the general shape of the curves. While simply monitoring the total number of wear particle atoms Nwp over time cannot prove that the formed wear particles are stable, a close scrutiny of simulation visualisations with coloring according to the advection velocity (not shown here) indicates that for configurations where appreciable numbers of wear particle 5
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65
(a)
(b)
(c)
(d)
(e)
(f)
(g)
(h)
(i)
Figure 2: Top views of surface topographies after 5 ns of sliding at 16 m/s. Colored according to atomic z coordinate, where the peaks are red and the pits are blue. Top: three-body contact sliding, center: twobody contact sliding with spherical abrasives, bottom: two-body contact sliding with cubic abrasives, left: 0.1 GPa, center: 0.5 GPa, right: 1.0 GPa.
6
Figure 3: Sketch of the zones to which the atoms can be attributed. atoms are produced, the respective wear particles are indeed stable and the shear zones (in which plastic flow occurs) are quite thin by comparison. What looks like a saturation of Nwp in Fig. 5 (a) actually marks the point where the running-in of the wear process has ended and the steady-state regime begins, a topic of ongoing research. In Fig. 5 (b), a sharp initial increase of the number of atoms in the shear zones takes place within the half-width of the employed Gauss filter and therefore cannot be resolved adequately. This is followed by a decrease over time as the asperity flanks gradually become less steep, with Nshear eventually becoming a constant proportional to the average shear force. 16 shear zone wear particle bulk average
14 12 10 g(r) [−]
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65
8 6 4 2 0 0
0.1
0.2
0.3
0.4 r [nm]
0.5
0.6
0.7
Figure 4: Comparison of time-averaged radial distribution functions calculated for the wear particles, the shear zone, and the bulk of the system with cubic abrasives, 1.0 GPa and 16 m/s.
4. Contact zone The smooth particle based estimation of the real contact area, as proposed by Eder et al. and discussed or applied in Refs. [24, 25, 26, 27], is quite time consuming for larger systems and currently reliable for singleasperity contacts only. We therefore quantify the contact zone via the number of atoms of one counterbody currently in contact with the other one, which according to Refs. [24, 25] is reasonably proportional to the contact area, irrespective of load and asperity shape. Here the problem of the definition of the counterbodies arises. In this work, we define the atoms of the abrasives together with the wear particles as the upper counterbodies. The lower counterbody consists of all other atoms which lie no more than 0.3 nm, the 7
4
4
x 10
2.5
2
2
1.5
1.5
Nwp [−]
Nwp [−]
2.5
1
1
0.5
0.5
0 0
20
40 ∆x [nm]
60
x 10
0 0
80
0.2
(a)
16000
16000
14000
14000
12000
12000
10000 8000
4000
4000
2000
2000 40 ∆x [nm]
0.8
1
0.8
1
8000 6000
20
0.6 σz [GPa]
10000
6000
0 0
0.4
(b)
Nshear [−]
Nshear [−]
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65
60
80
(c)
0 0
0.2
0.4
0.6 σz [GPa]
(d)
Figure 5: Sliding distance dependence (left) and averaged pressure dependence (right) of the number of atoms in wear particles (a–b) and in the shear zone (c–d). Green hues denote two-body contact sliding with cubic abrasives, red ones denote two-body contact sliding with spherical abrasives, and blue ones three-body contact sliding. The intensity of the color reflects the normal load. Dashed lines and small symbols represent sliding at 8 m/s, solid lines and large symbols at 16 m/s.
8
maximum contact distance, lower than the lowest atom of the upper counterbodies. These groups of atoms are well-defined at all times where there exists a filtered velocity profile. For computational efficiency, we have reduced the large number of atoms in the abrasives by considering only their bounding atoms, identified via their centro-symmetry parameter [28] which is larger than 5 in the initial configuration (before contact occurs). We then search for any atoms in the lower counterbody which lie within a distance dcont of an abrasive or a wear particle. This maximum contact distance dcont is chosen in such a way that the contact criterion includes all first (dcont = 0.265 nm), or first and second neighbors (dcont = 0.3 nm), cf. the radial distribution function in Fig. 4. The differences in the number of contact atoms between these two variants are small and nearly constant for constant load, although they do vary with load. In the two-body contact sliding simulations, almost all contact atoms lie within the shear zone, while in the three-body contact sliding simulations, quite a large proportion of the contact atoms are not sheared. This illustrates how free rolling abrasives hardly impose any shear on the counterbody, especially at low loads. Fig. 6 shows a comparison of the number of contact atoms with the number of shear zone atoms for three-body contact sliding at low load and two-body contact sliding at high load. Note that during three-body contact sliding, the number of contact atoms exceeds the number of shear zone atoms by a large margin. Some exemplary contact zone maps can be seen in Fig. 7, where the coloring according to height value gives an idea of the individual contact zone geometries. In Fig. 8 we compare the sliding distance and the pressure dependence of the total number of contact atoms (found for dcont = 0.265 nm) between the three considered configurations. The number of contact atoms initially decreases with the sliding distance, especially for medium and high pressures, but reaches a steady state towards the end of the simulations, and the pressure dependence is almost linear.
incl. 1st neighbors incl. 2nd neighbors all shear zone atoms
16000 14000
16000 14000 12000 Ncontact [−]
12000 Ncontact [−]
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65
10000 8000
10000 8000
6000
6000
4000
4000
2000
2000
0 0
20
40 ∆x [nm]
60
80
(a)
0 0
20
40 ∆x [nm]
60
80
(b)
Figure 6: The number of contact atoms compared with the shear zone atoms at 16 m/s. (a) Three-body contact sliding at 0.1 GPa, (b) two-body contact sliding with cubic abrasives at 1.0 GPa. By tracking which contact atom is in contact with which abrasive particle, it is possible to clusterize the contributions of the individual abrasives to the total contact zone. However, this only works while the abrasives are in direct contact with the lower counterbodies, which is the predominant case in three-body contact sliding simulations. During two-body contact sliding, the formation of wear particles soon shifts the main part of the contact area towards contact between wear particle and lower counterbody, see Fig. 9 (a). In this case, the direct affiliation of a contact atom touching a wear particle with one of the abrasives is not as straightforward as for three-body contact sliding. This can be overcome by iteratively searching for contact atoms which lie closer than dclust = 0.4 nm to other contact atoms whose abrasive affiliation is already known. The distance dclust can be chosen somewhat larger than dcont to reduce the number of required 9
(a)
(b)
Figure 7: Contact zone maps after '4 ns of sliding at 16 m/s. (a) Three-body contact sliding at 0.1 GPa, (b) two-body contact sliding with cubic abrasives at 1.0 GPa. Coloring is rainbow style according to height value of the contact atom, red is high, blue is low.
5000
5000
4000
4000 Ncontact [−]
Ncontact [−]
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65
3000
3000
2000
2000
1000
1000
0 0
20
40 ∆x [nm]
60
80
(a)
0 0
0.2
0.4
0.6 σz [GPa]
0.8
1
(b)
Figure 8: Sliding distance (a) and averaged pressure dependence (b) of the number of contact atoms. Color scheme as in Fig. 5.
10
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65
iterations, but at some point there is a risk of affiliations incorrectly spreading to neighboring clusters. The number of affiliated contact atoms may stop growing due to problems arising from the periodic boundary conditions, namely unaffiliated atoms with their affiliated neighbors at the other side of the simulation cell. This can be solved by once mapping all atoms in question which lie closer than dclust to the lower cell boundary across the upper cell boundary. Thus for the next iteration, a sufficient number of affiliated atoms is available on the correct side of the cell boundary so that more affiliations can be found. Finally, single outliers and small groups of outlying atoms can join the closest clusters with known affiliation. This last step is the only appreciable source of error, as isolated contact zones without a single initial particle affiliation will thus be falsely affiliated with neighboring abrasives.
(a)
(b)
(c)
Figure 9: Illustration of the clustering algorithm for the contact zone. Coloring is according to abrasive particle affiliation, and dark red is reserved for unaffiliated atoms in contact with wear particles. (a) Only atoms in direct contact with the numbered abrasive particles have known affiliations. (b) After clustering without taking the kinematics of the simulation into account, note the incorrect “backwards growth” of the light blue contact zone on the left. (c) Clustering corrected by taking simulation kinematics into account. Also observe that the light green area found at the top left due to periodic boundary conditions is correctly affiliated with the cluster at the top right. The above method is limited to non-coalescent contact clusters. As soon as the build-up of wear debris between two abrasives becomes sufficient to allow a single uninterrupted contact zone, it becomes a matter of definition which part of the contact should be identified with a particular abrasive. In this work, a wear particle related contact zone is attributed to the abrasive which formed the corresponding wear particle. In this case, the algorithm initially searches for affiliated atoms only in sliding direction, which effectively hinders the iterative “backwards growth” which would otherwise occur in coalescing contact zones, see Figs. 9 (b) and (c). The time-resolved knowledge of the size and the geometry of the individual contact zones allows an assessment of the influence of the process parameters on the contact situation. Fig. 10 shows that although the total number of contact atoms is very similar at constant normal load and sliding velocity, the distribution of the load and its time development varies considerably between the three considered surface finishing processes. A possible application of this clustering method is a time-resolved statistical analysis of abrasive indentation depths and orientations. Due to variations in surface topography and wear particle geometry, these parameters keep changing for every contact zone, and it is well known that they have considerable influence on real wear processes. These aspects as well as a sub-surface stress analysis are beyond the scope of this work. 11
3000
3500
2500
2500
3000
2000 1500 1000 500 0 0
Ncontact (per−particle) [−]
3000
Ncontact (per−particle) [−]
Ncontact (per−particle) [−]
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65
2000 1500 1000 500
20
40 ∆x [nm]
60
80
0 0
(a)
2500 2000 1500 1000 500
20
40 ∆x [nm]
60
80
(b)
0 0
20
40 ∆x [nm]
60
80
(c)
Figure 10: Break-down of the number of contact atoms by abrasive particle over sliding distance. Coloring is rainbow style with dark blue representing abrasive particle number 1 and red representing abrasive particle number 16. The three shown examples each represent one of the finishing methods at 0.5 GPa and 16 m/s. Left: three-body contact sliding, center: two-body contact sliding with spherical abrasives, right: two-body contact sliding with cubic abrasives. 5. Topography and wear evaluation For an efficient quantitative evaluation of the surface topography over time, the point cloud of atoms located on the surface must be determined and mapped to a regular grid. To achieve this, the simulation box was divided into square cells of side length dcell = 0.3966 nm, which is a trade-off between good lateral resolution and certainty of finding enough atoms in each cell. For every time step the substrate atom with the highest z value in each cell was then identified, i.e., all wear particle and shear zone atoms were excluded from this consideration. This procedure yields a time-resolved 2D representation zsubst (x, y; t) of the surface topography on a regular mesh of 72×72 points. Averaging along the x and y directions gives the mean substrate height for every time step, while the standard deviation of zsubst (x, y; t = const.) gives the Sq (t) roughness parameter at that time, both of which are shown as functions of sliding distance and normal load in Fig. 11. By monitoring the development of Sq (t), one can easily see the wear progress as the surface topography is flattened. Forming the pointwise difference ∆z(x, y; t) = zsubst (x, y; t0 ) − zsubst (x, y; t)
(1)
yields the wear height, which is positive where asperities are reduced in height, while it is negative where pits are filled up with material. This can be presented in the form of a wear map as shown in Fig. 12, where the former regions are shown in red and the latter ones in blue. The positive and negative contributions can be summed up individually and multiplied with the xy area of the mesh elements to yield the volumes of (−) (+) (−) (+) asperity height reduction Vasp and pit fill-up Vpit , see Fig. 13. Here we can see that Vasp and Vpit both saturate over the course of the simulation depending on the pressure. Note that the pressure dependence of (+) Vpit has a weak maximum at σz = 0.5 GPa for the two-body contact sliding simulations. The smoothing of the atomic velocity profiles, which determines the categorization of the atoms as substrate, shear or wear atoms, limits the number of time steps available for this analysis, so half the smoothing window width is lost at the beginning and the end of each dataset. The wear volume is estimated on the basis of the substrate atoms only, and there is usually an imbalance (−) (+) between Vasp and Vpit . The difference between the two, ∆V , can be decomposed into contributions from the atoms which form wear particles stuck to the abrasives, and those currently under shear between the wear particles and the substrate. As these two atomic categories are known by number, not by volume, an appropriate per-atom volume has to be found as a conversion factor. This per-atom volume should be in the range of the equilibrium volume of bcc Fe of 0.0116 nm3 . 12
3.9
3.8
3.8
zsubst [nm] (final)
zsubst [nm]
3.9
3.7 3.6
3.4 0
3.7 3.6 3.5
3.5
20
40 ∆x [nm]
60
3.4 0
80
0.2
0.4
(a)
0.7
0.6
0.6 Sq [nm] (final)
0.7
0.5 0.4
0.2
0.2 40 ∆x [nm]
1
0.8
1
0.4 0.3
20
0.8
0.5
0.3
0
0.6 σz [GPa] (b)
q
S [nm]
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65
60
80
(c)
0
0.2
0.4
0.6 σz [GPa] (d)
Figure 11: Sliding distance (left) and pressure dependence (right) of the mean surface height (a–b) and the Sq roughness parameter (c–d). Color scheme as in Fig. 5.
13
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65
(a)
(b)
Figure 12: Wear map after 4.8 ns of three-body contact sliding at 0.1 GPa and 16 m/s (a) and two-body contact sliding with cubic abrasives at 1.0 GPa and 16 m/s (b). Red areas denote regions of maximum height reduction (∆z (max) = +2.5 nm), blue areas denote regions of maximum fill-up (∆z (min) = −2.5 nm), while in the green areas there was no net change of height. However, in some cases, especially in three-body contact sliding simulations, ∆V cannot be neatly decomposed into these two contributions. There exists a small third, monotonically increasing and saturating contribution to the volumetric gap, which seems to arise from dynamic compression of the sub-surface substrate layer. To quantify this, we count the number of atoms in a 0.5 nm thick substrate layer immediately below the atom which defines zsubst (x, y; t), giving Nsubsurf (x, y; t) on the regular mesh introduced for the topography evaluation. As this layer has a constant volume over time, we can thus eliminate the error introduced by having to calculate the P total volume of the substrate. The total volume of the sub-surface layer should initially contain some x,y Nsubsurf (x, y; t0 ) ' 35000 atoms. However, the number actually obtained is closer to 34400 due to the influence of the topography and the somewhat crude resolution of the mesh. It may therefore be assumed that the exactness of this estimation method for the sub-surface compression is of the order ±3%. By subtracting the initial value from the current one, we obtain the number of surplus atoms in the sub-surface layer, X X Nsubsurf (x, y; t) − Nsubsurf (x, y; t0 ). (2) ∆Nsubsurf (t) = x,y
x,y
The fact that this sub-surface compression plays a non-negligible role predominantly during three-body contact sliding, where it can constitute the major contribution to ∆V , corresponds well to the wear regimes discussed in Ref. [15], where free rolling abrasives primarily impose pressure on the substrate. The sliding distance and pressure dependence of ∆Nsubsurf (t) is shown in Fig. 14. Here, we see evidence that the sharp-edged cubic abrasive particles lead to cutting wear, predominantly acting in the lateral direction and producing little normal compression, while the substrates underneath the spherical abrasive particles are compressed more. Now that the number of atoms are known which constitute the wear particles, Nwp (t), those that are in the shear zone, Nshear (t), and those additionally compressed into the sub-surface layer, ∆Nsubsurf (t) for every time step, we can fit per-atom volumes Va1 and Va2 as coefficients to these numbers to yield the difference between the total asperity volume reduction and the total pit fill-up volume such that ∆V (t) = Va1 [Nwp (t) + Nshear (t)] + Va2 ∆Nsubsurf (t).
(3)
All values obtained for Va1 lie close to the expected bulk value regardless of load, while the values of Va2 can vary substantially for small loads. However, a comparison with the time development of the per-atom 14
500
400
400 V(−) [nm3] asp
V(−) [nm3] asp
500
300
300
200
200
100
100
0 0
20
40 ∆x [nm]
60
0 0
80
0.2
250
250
200
200
150
100
0.8
1
0.8
1
0.8
1
150
100
20
40 ∆x [nm]
60
0 0
80
0.2
(c)
400
400
350
350
300
300 ∆V [nm3]
450
250 200
200 150
100
100
50
50 40 ∆x [nm]
0.6 σz [GPa]
250
150
20
0.4
(d)
450
0 0
0.6 σz [GPa]
50
50
0 0
0.4
(b)
V(+) [nm3] pit
V(+) [nm3] pit
(a)
∆V [nm3]
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65
60
80
(e)
0 0
0.2
0.4
0.6 σz [GPa] (f)
Figure 13: Sliding distance (left) and averaged pressure dependence (right) of the asperity height reduction volume (a–b), the pit fill-up volume (c–d) and the difference between the two (e–f). Color scheme as in Fig. 5. 15
3000
3000
2500
2500
2000
2000 ∆Nsubsurf [−]
∆Nsubsurf [−]
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65
1500 1000 500
1500 1000 500
0
0
−500
−500
−1000 0
20
40 ∆x [nm]
60
80
−1000 0
0.2
0.4
0.6 σ [GPa]
0.8
1
z
(a)
(b)
Figure 14: Sliding distance (a) and averaged pressure dependence (b) of the surplus atoms in the 0.5 nm thick sub-surface substrate layer compared with the initial state. Color scheme as in Fig. 5. volume for the entire substrate shows that overall number density changes are within a range of ±0.5% for all systems. With the knowledge of the coefficients Va1 and Va2 , the asperity volume reduction can now be decomposed into volumetric contributions of matter filled in the pits, wear particles stuck to the abrasives, matter currently under shear, and regions with sub-surface compression, see Fig. 15. This figure compares three-body contact sliding at low pressure with two-body contact sliding at high pressure with cubic abrasives. It is evident that the difference between the asperity reduction and the pit fill-up volumes varies considerably between these two processes. Virtually no wear particles are generated during three-body contact sliding, and ∆V consists of contributions from sub-surface compression and the shear zone. During two-body contact sliding, about half of the asperity reduction volume consists of wear particles, and the rest is evenly divided between material filled into the pits and the shear zone, while sub-surface compression is virtually negligible. By adapting the clustering algorithm described above to the identification of any wear particle with the abrasive that caused it, we can further break down the wear particles into contributions from the individual abrasives. The three examples shown in Fig. 16 give an impression of how this analysis method can give an overview over the distribution of wear volume formation. In the three-body contact sliding simulation, shown in Fig. 16 (a), wear particles are essentially produced by two abrasives only, while two-body contact sliding with cubic abrasives, shown in Fig. 16 (c), leads to a very even distribution of wear particles. Here, 12 out of 16 abrasives produce wear particles of similar volume which remain quite stable until the end of the simulation. In the data obtained for two-body contact sliding with spherical abrasives, shown in Fig. 16 (b), we can see a rather irregular distribution of the wear particles among the abrasives as well as evidence of the formation of transient wear particles, i.e., wear particles that are deposited in pits at some point and thus become part of the substrate again (e.g., the light blue area in the center portion). 6. Summary and Conclusion We have presented a method which allows the quantitative evaluation of the wear volume and the real contact zone resulting from two-body and three-body wear of nanoscopically rough surfaces using MD simulations. The focus of this work lay on implementing a multi-abrasive approach to this problem which enables us to study the effects of multiple abrasives in close proximity and changing surface topography. 16
500 400 3
400
asperity reduction pit fill−up wear particles shear zone sub−surface compression Vwear [nm ]
500
Vwear [nm3]
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65
300 200 100 0 0
300 200 100
20
40 ∆x [nm]
60
80
(a)
0 0
20
40 ∆x [nm]
60
80
(b)
Figure 15: Sliding distance dependence of the volumetric contributions to the total asperity reduction volume (shown in blue) at 16 m/s. (a) Three-body contact sliding with spherical abrasive particles, 0.1 GPa. (b) Two-body contact sliding with cubic abrasive particles, 1.0 GPa. In particular, we outlined a scheme to generate an initial surface topography with a pseudo-random Gaussian roughness distribution which can be easily configured and scaled while being periodically replicable. We discussed which constraints on the individual abrasive particles mimic certain surface finishing process conditions. An analysis of the atomic advection velocities allowed the dynamic identification of wear particles, the substrate and the shear zone in-between. We then presented a method for the time-resolved quantification of the real contact zone between the abrasives and the wear particles on the top, and the shear zone and substrate on the bottom. Using a clustering technique which takes into account the kinematics of the process, we could break down the multi-asperity contact zone into its single-asperity contributions. Using a mesh-based topography evaluation, we could generate time-dependent wear maps and break down the asperity reduction volume into its contributions, including the pit fill-up volume, the wear particles, the shear zone, and volume lost through substrate compression. By applying the clustering technique to the wear particles themselves, it was possible to analyze the contributions of the individual abrasives to wear particle formation. Acknowledgements This work was funded by the Austrian COMET-Program (Project K2 XTribology, Grant No. 824187), the ERDF as well as the province of Lower Austria (Onlab project), and it was carried out at the “Excellence Centre of Tribology”. S. J. Eder would like to thank Georg Vorlaufer for the helpful discussion about the sub-surface compression of the substrate. References [1] I. Szlufarska, M. Chandross, R. W. Carpick, Journal of Physics D: Applied Physics, 41 (2008) 123001. [2] E. Brinksmeier, J.C. Aurich, E. Govekar, C. Heinzel, H.-W. Hoffmeister, F. Klocke, J. Peters, R. Rentsch, D.J. Stephenson, E. Uhlmann, K. Weinert, M. Wittmann, CIRP Annals-Manufacturing Technology, 55 (2006) 667. [3] R. Rentsch, I. Inasaki, CIRP Annals-Manufacturing Technology, 43 (1994) 327. [4] K. Maekawa, A. Itoh, Wear, 188 (1995) 115. [5] L. Zhang, H. Tanaka, Wear, 211 (1997) 44. [6] R. Komanduri, N. Chandrasekaran, L. M. Raff, Wear, 240 (2000) 113. [7] Q. X. Pei, C. Lu, H. P. Lee, Y. W. Zhang, Nanoscale research letters, 4 (2009) 444.
17
3500
14000
600
3000
12000
500 400 300 200 100 0 0
Nwp (per−particle) [−]
700
Nwp (per−particle) [−]
Nwp (per−particle) [−]
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65
2500 2000 1500 1000 500
20
40 ∆x [nm]
60
80
0 0
10000 8000 6000 4000 2000
20
40 ∆x [nm]
60
80
0 0
20
40 ∆x [nm]
(a)
(b)
(c)
(d)
(e)
(f)
60
80
Figure 16: Top row: break-down of the number of wear particle atoms Nwp (t) by wear particle over sliding distance. Bottom row: top view of the wear particle distribution at ∆x = 40 nm (t = 2.5 ns). Coloring is rainbow style with dark blue representing wear particle number 1 and red representing wear particle number 16. The three shown examples each represent one of the finishing methods at 0.5 GPa and 16 m/s. left: three-body contact sliding, center: two-body contact sliding with spherical abrasives, right: two-body contact sliding with cubic abrasives.
18
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65
[8] [9] [10] [11] [12] [13] [14] [15] [16] [17] [18] [19] [20] [21] [22] [23] [24] [25] [26] [27] [28]
L. Si, D. Guo, J. Luo, X. Lu, Journal of Applied Physics, 107 (2010) 064310. P. M. Agrawal, L. M. Raff, S. Bukkapatnam, R. Komanduri, Applied Physics A, 100 (2010) 89. J. Zhong, R. Shakiba, J. B. Adams, Journal of Physics D: Applied Physics, 46 (2013) 055307. P. Spijker, G. Anciaux, J.-F. Molinari, Tribology International, 59 (2013) 222. X. Han, Y. Hu, S. Yu, Applied Physics A, 95 (2009) 899. R. Chen, R. Jiang, H. Lei, M. Liang, Applied Surface Science, 264 (2012) 148. L. Zhang, H. Zhao, Z. Ma, H. Huang, C. Shi, W. Zhang, AIP Advances, 2 (2012) 042116. L. Zhang, H. Tanaka, Tribology International, 31 (1998) 425. J. Sun, L. Fang, J. Han, Y. Han, H. Chen, K. Sun, Computational Materials Science, 82 (2014) 140. P. M. Agrawal, R. Narulkar, S. Bukkapatnam, L. M. Raff, R. Komanduri, Tribology International, 43 (2010) 100. M. I. Mendelev, S. Han, D. J. Srolovitz, G. J. Ackland, D. Y. Sun, M. Asta, Philosophical Magazine, 83 (2003) 3977. T.-H. Fang, C.-I. Weng, Nanotechnology, 11 (2000) 148. Q. X. Pei, C. Lu, F. Z. Fang, H. Wu, Computational materials science, 37 (2006) 434. ISO 25178-2:2012 Geometrical product specifications (GPS) – Surface texture: Areal – Part 2: Terms, definitions and surface texture parameters. August 2012. J. R. D’Errico, Understanding gridfit, 2006. URL: http://www.mathworks.at/matlabcentral/fileexchange/8998-surfacefitting-using-gridfit (Last visited: Feb. 7, 2014). S. J. Plimpton, J. Comput. Phys., 117 (1995) 1. S. Eder, A. Vernes, G. Vorlaufer, G. Betz, J. Phys.: Condens. Matter, 23 (2011) 175004. S. Eder, A. Vernes, G. Betz, Computer Physics Communications, 185 (2014) 217. S. J. Eder, A. Vernes, G. Betz, Langmuir, 29 (2013) 13760. A. Vernes, S. Eder, G. Vorlaufer, G. Betz, Faraday Discuss., 156 (2012) 173. C. L. Kelchner, S. J. Plimpton, J. C. Hamilton, Phys. Rev. B, 58 (1998) 11085.
19