- Email: [email protected]
Sliding and deformation of metallic glass: experiments and MD simulations
Journal of Non-Crystalline Solids 317 (2003) 206–214 www.elsevier.com/locate/jnoncrysol
Sliding and deformation of metallic glass: experiments and MD simulations Xi-Yong Fu a, D.A. Rigney
b,*
, M.L. Falk
c
a Merck Research Laboratories, WP78-304, Sumneytown Pike, West Point, PA 19486, USA Materials Science and Engineering, The Ohio State University, 2041 College Rd., Columbus, OH 43210, USA Materials Science and Engineering, and Applied Physics, University of Michigan, Ann Arbor, MI 48109-2136, USA b
c
Abstract The sliding behavior of metallic glass (MG) was studied using pin-on-disk tests and molecular dynamics (MD) simulations. Friction coefficients and wear rates of a Zr–Ti–Cu–Ni–Be alloy were similar to those reported for ductile materials, e.g., normal crystalline metals; i.e., exceptional friction and wear characteristics of bulk MG were not observed. Sliding caused plastic deformation, transfer, mechanical mixing and reamorphizing of devitrified material. In vacuum, a softer layer developed adjacent to the interface. In air, an additional harder layer appeared when oxidation products were incorporated into the near-surface material. MD simulations involved a two-component 2D amorphous system. Simulations of tensile tests showed elastic/perfectly plastic response, strain rate dependence, void formation and shear bands. Simulations of sliding showed decreased density near the interface, suggesting an increase in free volume during shear, but neither voids nor shear bands. Subsurface displacement profiles were similar to those reported in experiments on crystalline materials and were consistent with flow patterns expected for flow near a boundary. The MD results on mechanical mixing suggest relevance to other processes, including mechanical alloying, friction welding, formation of nanocrystals, erosion and deformation at high strain rates. 2003 Elsevier Science B.V. All rights reserved. PACS: 81.05.Kf; 02.70.Ns; 46.55.+d; 83.50.)v; 47.32.)y; 61.43.Dq; 62.20.Qp; 62.20.Fe; 81.40.Pq; 81.40.Lm; 81.05.Ys; 61.46.+w; 61.16.Bg
1. Introduction Early studies of the tribological behavior of metallic glass used thin specimens produced by rapid quenching from the melt [1,2]. The recent availability of bulk metallic glass (BMG) [3,4]
*
Corresponding author. Tel.: +1-614 292 1775; fax: +1-614 292 1537. E-mail address: [email protected] (D.A. Rigney).
simplifies tribological testing and allows more thorough chemical and structural characterization of worn specimens [5]. The experimental work reported here was undertaken to confirm reports, e.g., in [4], of favorable friction and wear characteristics of metallic glass. To complement the experimental work, molecular dynamics (MD) simulations of the deformation and sliding of a model amorphous solid were also performed. Others have used MD techniques to study unlubricated friction involving crystalline materials
0022-3093/03/$ - see front matter 2003 Elsevier Science B.V. All rights reserved. doi:10.1016/S0022-3093(02)01999-3
X.-Y. Fu et al. / Journal of Non-Crystalline Solids 317 (2003) 206–214
[6–11]. Those earlier simulations showed a reduction of coefficient of friction (COF) with increasing sliding speed, mechanical mixing and the development of nanocrystals adjacent to the sliding interface [10].
2. Experiments The alloy chosen was Zr41:2 Ti13:8 Cu12:5 Ni10:0 Be22:5 (composition in at.%) cast by Howmet. Its hardness was HV530 with a load of 0.3 kg and the shear modulus has been reported to be 33 GPa [12]. X-ray diffraction (XRD) and high resolution transmission electron microscope (TEM) imaging confirmed that the structure was amorphous. Sliding tests were done with a simple pin/disk arrangement, shown schematically in Fig. 1. The chamber could be evacuated, so tests could be done in air or vacuum (typically 4 103 Pa). Not shown in Fig. 1 is a Kelvin probe that was used in situ to monitor the combined effects of structural and chemical changes on the surface of the disk. Test pins were 3 3 12 mm rods of BMG or 6 mm diameter balls of 52 100 bearing steel. Disks were 31 mm in diameter and were 3 mm thick. Surfaces were prepared by metallographic polishing, using 1 lm diamond paste for the final stage and ultrasonic cleaning in acetone, methanol and distilled water. XRD confirmed that the prepared surfaces remained amorphous. Normal loads ran-
207
ged from 0.1 to 1.2 kg and sliding speeds were typically 0.05 m/s but some tests were run up to 1 m/s. Average wear rate was determined from measurements of weight changes caused by sliding. Structural and chemical characterization of worn surfaces and debris was accomplished by optical microscopy, scanning electron microscopy (SEM), energy dispersive X-ray Spectroscopy (EDS), XRD, TEM and microhardness. Other details of the equipment and procedure are presented in references [5] and [13].
3. MD simulations No attempt was made to model in three dimensions the five-component BMG material used in the experimental part of this study. Such a simulation would require interatomic potentials that are not available and would involve excessive computation time. Instead, a Lennard–Jones potential and a simple 2D amorphous system was used. It was anticipated that these simulations would exhibit some of the major features that a more realistic model would generate. The energy and length scale parameters in the potential were those reported in [14]. The composition was chosen such that the ratio of the number of L or ÔlargeÕ atoms to S or ÔsmallÕ atoms was ð1 þ 51=2 Þ=4, one half the golden mean, because this ratio is characteristic of a Penrose tiling of space by rhombic
Fig. 1. Experimental set-up for pin/disk tests.
208
X.-Y. Fu et al. / Journal of Non-Crystalline Solids 317 (2003) 206–214
units consisting of these particles (atoms) [14]. This choice optimizes for the existence of a quasicrystalline ground state that corresponds to this tiling. The chosen system is easily quenched to a metastable glass having a glass transition temperature Tg of 0.3 (rescaled dimensionless units) [15]. All simulation values reported in this paper are presented in these rescaled units. Details of the potential and the rescaled units are given in Appendix A and in [13–16]. The model system contained 4472 atoms designated L and 5528 atoms designated S, for a total of 10 000 atoms. These were contained in a square simulation cell with periodic boundary conditions. Temperature and pressure control were accomplished by using a Nose–Hoover thermostat and a Parinello–Rahman barostat [17]. After the system was equilibrated at high temperature (T ¼ 2:3) and pressure, the system was quenched below Tg . The quench sequence assured that voids were suppressed and that the product was a metastable 2D glass at a very low temperature (T ¼ 0:01) and zero pressure. The mechanical properties of the 10 000 atom model system were checked by running two kinds of tensile test simulations, one with periodic cell boundaries on all sides in order to eliminate surface effects and the other with this constraint relaxed along the edges parallel to the tensile axis. In both cases, zero stress boundary conditions were maintained perpendicular to the tensile axis, resulting in comparable uniaxial tests. For simulations of sliding, an identical set of 10 000 atoms was placed in contact with the first set, but shifted along the x axis, i.e., parallel to the interface by half of the cell dimension. This process generated a sliding pair consisting of 20 000 atoms. Periodic boundary conditions were retained at the cell boundaries normal to x but those conditions were relaxed at the other two boundaries, those normal to y. Also, an important slab of ÔreservoirÕ material was defined adjacent to each of those boundaries. Each of these reservoir regions was maintained at a low temperature of T ¼ 0:01, allowing frictional heat to conduct to these thermostat regions. Also, external forces were only applied to the atoms in those regions, thus simulating experiments in which the forces responsible for sliding
are imposed at some distance from the interface. To initiate sliding (t ¼ 0), all atoms in the upper block were given a velocity to the right ðþxÞ and all atoms in the lower block were given the same velocity to the left ðxÞ, providing a net interface sliding speed of vs . At t > 0, temperature and friction force control operated in the reservoirs. Force control involved a feedback loop so that sliding speed was maintained at vs . To provide a normal load, a pressure of 1.0 was applied in the y direction. Three aspects of this procedure are similar to those used earlier by Hammerberg et al. [7–9,11]: the use of two reservoirs, the periodic boundary conditions in the x-direction and the way in which sliding was initiated. The following aspects differed from those reported in [7–9,11]: the system composition and structure, the interatomic potentials, the method used to apply the shear force and the use of a constant pressure constraint rather than constant volume during sliding.
4. Results 4.1. Experimental results Sliding behavior of the metallic glass (MG) depended on load, speed and atmosphere. In slow speed tests, typically at 0.05 m/s, the steady state coefficient of friction (COF) in air for self-mated MG sliding decreased smoothly from 0:65 0:05 to 0:38 0:05 as the load increased from 0.2 to 1.1 kg. For the same load range, but in vacuum, the COF was somewhat higher, 0:82 0:05 to 0:60 0:05, but followed a similar trend. When the pin was changed to 52 100 steel, the values of COF in air and vacuum and the trends with load were similar to those found with MG pins, except that COF for 52 100 on MG in air was higher than for MG on MG at light loads. For comparison, a limited number of tests of 52 100 on 52 100 were done in air and vacuum. The COF values were close to the values found for tests with MG. In a series of tests at 0.5 kg normal load and sliding speeds ranging from 0.05 to 1.0 m/s, the COF decreased from 0:50 0:05 to 0:29 0:05. The highest wear rates were for MG on MG and 52 100 on MG in air. They ranged from 15 4
X.-Y. Fu et al. / Journal of Non-Crystalline Solids 317 (2003) 206–214
to 120 10 (expressed in 1010 kg/m) as the load varied from 0.2 kg to 1.1 kg. Wear rates in air with 52 100 as one or both components varied similarly with load, but were lower. The lowest wear rates were for MG on MG in vacuum (8 5 to 5 3, both times 1010 kg/m, for loads of 0.2–0.5 kg). SEM images and EDS maps of wear tracks after tests in air showed patches of material with enhanced oxygen content. These were consistent with Kelvin probe data obtained during sliding [18]. The groove and ridge topography was typical of worn surfaces of ductile materials. On most of the wear track, hardness had increased to HV580, but the regions with higher oxygen content were harder (HV640). Oxygen rich patches were absent after tests in vacuum and hardness values on the wear track were similar to those of the unworn material. Worn surfaces of pins were similar to those of disks. When the pin was 52 100 steel, transfer of MG from the disk was readily detected by SEM/EDS. The amount transferred was larger at higher loads or in vacuum. Debris particles were irregular flakes and plates ranging in size from 5 to 100 lm. Smaller particles with enhanced oxygen content were present after tests in air. Some cutting chips were also found, especially after tests in air. Sliding tests were also done on devitrified material obtained by annealing for 2 h at 500 C. The hardness after annealing was HV720 and XRD confirmed that the material was crystalline. Steadystate COFs and wear rates were similar to those obtained with amorphous material and XRD of the debris showed only broad peaks typical of amorphous material. SEM/EDS observations of longitudinal and transverse cross-sections (parallel and perpendicular to the sliding direction) revealed that sliding caused changes well below the surface, typically tens of microns. For specimens tested in air, surface patches with enhanced oxygen content were detected by back-scattered electron images and by EDS. The hardness in these regions was also enhanced (600–700 Knoop hardness, KHN). For specimens tested in vacuum, the material adjacent to the surface had decreased hardness compared with bulk material. That is, the material showed work-softening rather than work-hardening. The
209
softer region extended for 30–70 lm from the surface. After tests in air, a corresponding soft layer was also detected between the locally hard patches and the underlying base material. No differences in appearance or in chemical composition were detected between the softer material and the base material. 4.2. Results of MD simulations Fig. 2(a) shows the results of simulated tensile tests on the model MG system for a range of strain rates. For these simulations, periodic boundaries were maintained along the surfaces parallel to the tensile axis. The response is typical of an elasticperfectly-plastic material. The transition to plastic flow occurs at higher values of stress as strain rate is increased. The YoungÕs modulus varies from 28.4 to 32.9 (rescaled units) and the 2D PoissonÕs ratio from 0.55 to 0.50 as the strain rate is increased from 1 105 to 2 104 . Fig. 2(b) (lower strain rate) and Fig. 2(c) (higher strain rate) show that voids developed during the tensile test simulations. These were smaller and more randomly distributed for the high strain rate case. For the larger voids, it is clear that they nucleated where atoms of the same kind were clustered in the original material. Fig. 3(a) shows the results of a tensile simulation in which the periodic constraints at the cell walls parallel to the tensile axis were removed. This allowed these surfaces to roughen. In this case, the shading indicates regions in which strain is concentrated. The corresponding stress-strain curve has the same general form as those shown in Fig. 2(a), but with larger excursions from the general trend. Fig. 3(b) uses the same shading convention, but for the simulation already shown in Fig. 2(b). In that case, in addition to seeing the voids, one can see regions with localized strain, but they are not as strongly developed as in Fig. 3(a). Simulations of sliding of MG against MG showed COF rising quickly to a maximum value of 0:50 0:05 and then falling to a steady-state value that decreased from 0:31 0:05 to 0:06 0:03 as sliding speed increased from 0.1 to 1.0 (rescaled units). Fig. 4 shows a sequence in the interface region for a sliding speed of 1.0. The material above the
210
X.-Y. Fu et al. / Journal of Non-Crystalline Solids 317 (2003) 206–214
Fig. 3. (a) Structure produced by tensile strain at slower strain rate (1 105 , rescaled units), showing shear bands when surface roughening is allowed. Shaded regions have higher strain. (b) Same as 3a, but periodic boundary conditions imposed at sides of cell, showing voids and incipient shear bands.
Fig. 2. (a) Tensile stress–strain curves at different strain rates (1 105 to 2 104 , rescaled units) from MD simulations. Periodic boundary conditions imposed at sides of cell. (b) Structure produced by tensile strain at slower strain rate (1 105 , rescaled units), showing large voids. (c) Similar to 2b, but at high strain rate (2 104 , rescaled units), showing smaller voids.
interface is color coded differently from the material below the interface but is otherwise the same. The initially smooth interface first roughens and then a mixing process is observed. A movie version of this sequence shows that mechanical mixing
near the interface proceeds by a random process involving rotations of irregular clusters of atoms. That is, it is not a simple atomic diffusion process. Closer to the reservoir regions, the material flow approximates laminar flow parallel to the sliding direction. The color coding allows one to monitor the progress of mixing, and, in particular, the growth kinetics of the mixed layer. This was done by calculating the net displacements along the y direction (normal to the interface) relative to the center of mass for atoms in thin slabs parallel to the sliding direction. Each resulting profile was fit to a Gaussian function centered at the interface, allowing an objective determination of the thickness of the mechanically mixed material at different times. As shown in Fig. 6, the thickness of this material was proportional to t1=2 . A plot of the
X.-Y. Fu et al. / Journal of Non-Crystalline Solids 317 (2003) 206–214
211
Fig. 4. Sequence showing interface roughening and mixing during sliding for t ¼ 0 and two later times (rescaled).
strain rate vs. distance from the interface can also be fit to a Gaussian function centered at the interface, and the thickness of the high strain rate region is also proportional to t1=2 . The transients observed in simulations of COF suggest that changes in the near-surface material occur during sliding. The simplest change to investigate would be changes in density. Therefore, the density in thin layers parallel to the interface was compared with the original density. The density difference profile is shown in Fig. 5. This shows a clear decrease in density in the most highly deformed material. Part of this density decrease is Fig. 6. Density difference profile relative to average density.
Fig. 5. Growth kinetics of the thickness of the mixed layer, from MD simulations.
from thermal expansion, but when this smaller contribution is subtracted, the decrease in density near the interface remains. The slight increase in density elsewhere in the cell arises because of the use of constant volume conditions. The highest local temperature at t ¼ 25, determined from the random part of the kinetic energy in each thin layer, was 0.1, well below the Tg of 0.3. In experimental work on crystalline materials, various features associated with the microstructure can be used as markers to measure displacement profiles adjacent to sliding interfaces. These are not available in MG, but color coding provides convenient markers in MD simulations. Fig. 7 shows how an initially straight marker bends over
212
X.-Y. Fu et al. / Journal of Non-Crystalline Solids 317 (2003) 206–214
Fig. 8. Displacement profiles at four different times, from MD simulations.
Fig. 7. Simulation of bending of an originally straight marker during sliding.
chosen time. The results show that profiles from MD simulations and from Navier–Stokes analysis are in agreement [16].
5. Discussion in response to the applied shear stresses for sliding of MG on MG. This profile is similar to displacement profiles observed for a range of size scales in experiments on various single and multiphase crystalline materials. In fact, they all resemble the profiles expected for the case of a viscous fluid near a moving boundary. This would be described by the Navier–Stokes equation with suitable boundary conditions. The familiar solution for the fluid velocity as a function of distance from the interface is a complementary error function and this yields a similar displacement profile. If one accounts for the elastic to plastic transition in the MG, one can generate a time sequence of displacement profiles, as shown in Fig. 8. After an initial transient, the elastic strain saturates and the plastically deformed zone reaches a steady state depth of about 50 length units. If one subtracts the elastic contributions to these curves, one can fit the simulation data for a chosen time to the profile predicted by Navier–Stokes analysis. This gives a value for the effective viscosity (assuming it remains constant), and that parameter in turn allows prediction of the displacement profile for any
The unlubricated friction and wear results reported here suggest that bulk metallic glass has tribological behavior similar to that of more familiar crystalline materials. In both cases the loading conditions during sliding are compression and shear; these are conditions that inhibit fracture and encourage ductility. The results may be quite general for systems in which plastic deformation continues during sliding. Blau has recently reported similar results and conclusions for a Zr–Cu–Ni–Ti– Al BMG [19]. Earlier expectations that BMGÕs might have exceptional sliding characteristics were based largely on experiments with thin ribbons of Fe-based MG [2,20]. It may be that MG materials will exhibit a range of sliding behavior, just as crystalline materials do, depending on composition and the presence or absence of other phases. It is clear that some oxidation occurs during sliding in air, but this does not produce an oxide layer that yields oxide wear debris directly. Instead, the oxide that forms is blended into the deforming MG, and it is this mixed material that generates the oxide content in wear debris. Different amounts of
X.-Y. Fu et al. / Journal of Non-Crystalline Solids 317 (2003) 206–214
mixed in oxide cause differences in local hardness on the wear surface and in transfer material on the pin. If the local hardness difference at contacting asperities becomes greater than 20%, cutting debris can be generated [21]. No evidence of sliding-induced devitrification was found. In fact, the data from tests on annealed BMG showed the opposite effect: the deformation produced by sliding reamorphized the material. This could be a useful characteristic of BMG – if any local region did crystallize, perhaps in response to temporary thermal excursions, continuing sliding would restore the amorphous structure. The soft layer created by sliding is very interesting for at least two reasons. First, it assures that, in the absence of outside influences that can raise the surface hardness, sliding on BMG produces and maintains a soft surface layer on a harder base. This controls friction and correlates with relatively smooth sliding. Second, the data are consistent with suggestions that increases in free volume are associated with plastic deformation of MG [21,22]. The regions of decreased density predicted by the MD simulations would also be consistent with this view. The shape of the tensile stress–strain curves shown in Fig. 2(a) is expected for a material that lacks the microstructural features that would allow work-hardening. The smoother curves at higher strain rates would be expected for a structure that is more homogeneous, i.e., having smaller and more evenly distributed voids. The voids shown in Fig. 2(b) are larger because they have more time to grow at slower strain rates. The nucleation of the voids is affected by local fluctuations in composition. They occur more easily where there is a local concentration of weaker bonds. In this simulation, bonds between two L atoms or two S atoms are weaker than those between unlike atoms. The results of tensile test simulations depend also on conditions at the surface. When the surface is allowed to roughen, shear bands develop and void formation is inhibited. The corresponding stressstrain curve is consistent with the appearance of shear bands. A transient peak followed by a decrease to a steady state value are common features of experimental COF vs. time plots. The MD results predict
213
similar behavior. Experimental measurements show a decrease in COF as sliding speed increases. A similar pattern results from the MD simulations. MD simulations with crystalline copper show similar Ôvelocity softeningÕ [10]. Even more striking are the observations that mechanical mixing occurs in sliding experiments and in MD simulations for both MG and crystalline material [10]. In the simulations on MG, this mixing is associated with nanoscale rotational flow in material where the strain rates are highest. In this region, the collective effect of continuing shear-induced microstructural rearrangements in the amorphous material can lead to long-range mass transport normal to the principal shear direction. Rotational flow and mixing may also be responsible for the reamorphizing observed for sliding of devitrified MG.
6. Conclusions This project began with an emphasis on the sliding friction and wear behavior of bulk metallic glass. The experimental observations, combined with the results of molecular dynamics simulations, suggest that the flow processes involved are more general than originally thought. It is suggested here that amorphous materials are well suited to serve as model materials, both for experiments and simulations, allowing identification of generic properties and behavior in a broader range of materials and processes. The insight obtained from these studies may be helpful for understanding the behavior of both amorphous and crystalline materials when large plastic strains are imposed at high strain rates. Potential areas of application include friction and wear, high strain rate deformation, mechanical alloying and friction stir welding. Earlier reports [2,20] that BMGÕs might provide exceptional friction and wear characteristics are not confirmed by the results of the present work.
Acknowledgements We are pleased to acknowledge research support from the National Science Foundation
214
X.-Y. Fu et al. / Journal of Non-Crystalline Solids 317 (2003) 206–214
(NSF), Surface Engineering and Tribology, most recently under Grant no. CMS-9812854, NSF, Division of Materials Research under Grant no. DMR-0135009, and from Howmet Corp. We are grateful to Sharon Glotzer, then at NIST), The Ohio Supercomputer Center and the University of MichiganÕs Center for Parallel Computing for the use of computer facilities, and we thank the following for helpful discussions: M.J. Mills, P.K. Gupta and A. Markworth at The Ohio State University, A. Zharin at the Belarussian State Powder Metallurgy Concern, Minsk, Belarus, and J. Hammerberg at Los Alamos National Laboratory. Appendix A
If the system simulated were pure Cu, the relative ; e ¼ 0:35279 coefficients would be r ¼ 2:5172 A 25 kg. eV; m ¼ 1:0551 10 Based on these values, the following relationship can be calculated ; r^ ¼ 1 () 2:5172 A
^t ¼ 1 () 3:4423 1013 s; P^ ¼ 1 () 3:5392 GPað3DÞ or P^ ¼ 1 () 0:89087N=m ð2DÞ; v^ ¼ 1 () 731:42m=s:
References
The Lennard–Jones pair potential has the form r 12 r 6 U ¼ 4e ; r r where r is the distance between two atoms; r is the distance at which the interatomic potential U is zero and e is the bond strength. The MD simulation used Lennard–Jones potentials with units rescaled such that all quantities are expressed relative to the scaling factor. After rescaling, the simulation results are not associated with any particular system; different systems would only affect the scaling factors, leaving the rescaled simulation results unchanged. The rescaling factors for some of the most commonly used quantities are given below:
[1] [2] [3] [4] [5] [6] [7] [8] [9] [10] [11] [12] [13] [14]
• BoltzmannÕs constant and each atomic mass are set equal to 1.0. • Length is rescaled such that ^ r ¼ r=r where ^ r is the rescaled length and r is the actual length. • Energy is rescaled such that e^ ¼ e=e. pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 2 =ÞÞ. • Time is rescaled such that ^t ¼ t=ð ðmr pffiffiffiffiffiffiffiffi • Velocity is rescaled such that v^ ¼ v=ð e=mÞ. • Pressure, stress and modulus are rescaled such that P^ ¼
P for 3D or ðe=r3 Þ
e^ ¼ 1 () 0:35279 eV;
P^ ¼
P for 2D: ðe=r2 Þ
[15] [16] [17] [18]
[19] [20] [21] [22]
P.G. Boswell, J. Mater. Sci. 14 (1979) 1505. K. Miyoshi, D.H. Buckley, ASLE Trans. 27 (1983) 295. W.L. Johnson, MRS Bull. 24 (1999) 42. A. Inoue, Acta Mater. 48 (2000) 279. X.-Y. Fu, T. Kasai, M.L. Falk, D.A. Rigney, Wear 250 (2001) 409. G. He, M.O. Robbins, Science 284 (1999) 1650. J.E. Hammerberg, B.L. Holian, J. Roeder, A.R. Bishop, S.J. Zhou, Phys. D 123 (1998) 330. B.L. Holian, J.E. Hammerberg, P.S. Lomdahl, J. Comp. Aided Mater. Design 5 (1998) 207. J.E. Hammerberg, B.L. Holian, S.J. Zhou, Chem. Phys. Rep. 17 (1998) 343. D.A. Rigney, J.E. Hammerberg, MRS Bull. 23 (1998) 32. J. Roeder, A.R. Bishop, B.L. Holian, J.E. Hammerberg, R.P. Mikulla, Phys. D 142 (2000) 306. W.L. Johnson, Mater. Sci. Forum 225–227 (1996) 35. X.-Y. Fu, PhD dissertation, The Ohio State University, Columbus, OH, 2001. F. Lancon, L. Billard, P. Chaudhari, Europhys. Lett. 2 (1986) 625. M. Falk, J. Langer, Phys. Rev. E 57 (1998) 7192. X.-Y. Fu, M.L. Falk, D.A. Rigney, Wear 250 (2001) 420. S. Melchionna, G. Ciccotti, B.L. Holian, Mol. Phys. 78 (1993) 533. T. Kasai, X.-Y. Fu, T. Honda, A. Zharin, D.A. Rigney, Proceedings of International Tribology Conference, Nagasaki, 2000, Japan Society of Tribologists, Tokyo, 2001, p. 333. P.J. Blau, Wear 250 (2001) 431. T. Imura et al., Mater. Sci. Eng. A 133 (1991) 332. T. Akagaki, D.A. Rigney, Wear 149 (1991) 353. D. Deng, A.S. Argon, S. Yip, Philos. Trans. R. Soc. Lond. A 329 (1989) 613.
Sliding and deformation of metallic glass: experiments and MD simulations Xi-Yong Fu a, D.A. Rigney
b,*
, M.L. Falk
c
a Merck Research Laboratories, WP78-304, Sumneytown Pike, West Point, PA 19486, USA Materials Science and Engineering, The Ohio State University, 2041 College Rd., Columbus, OH 43210, USA Materials Science and Engineering, and Applied Physics, University of Michigan, Ann Arbor, MI 48109-2136, USA b
c
Abstract The sliding behavior of metallic glass (MG) was studied using pin-on-disk tests and molecular dynamics (MD) simulations. Friction coefficients and wear rates of a Zr–Ti–Cu–Ni–Be alloy were similar to those reported for ductile materials, e.g., normal crystalline metals; i.e., exceptional friction and wear characteristics of bulk MG were not observed. Sliding caused plastic deformation, transfer, mechanical mixing and reamorphizing of devitrified material. In vacuum, a softer layer developed adjacent to the interface. In air, an additional harder layer appeared when oxidation products were incorporated into the near-surface material. MD simulations involved a two-component 2D amorphous system. Simulations of tensile tests showed elastic/perfectly plastic response, strain rate dependence, void formation and shear bands. Simulations of sliding showed decreased density near the interface, suggesting an increase in free volume during shear, but neither voids nor shear bands. Subsurface displacement profiles were similar to those reported in experiments on crystalline materials and were consistent with flow patterns expected for flow near a boundary. The MD results on mechanical mixing suggest relevance to other processes, including mechanical alloying, friction welding, formation of nanocrystals, erosion and deformation at high strain rates. 2003 Elsevier Science B.V. All rights reserved. PACS: 81.05.Kf; 02.70.Ns; 46.55.+d; 83.50.)v; 47.32.)y; 61.43.Dq; 62.20.Qp; 62.20.Fe; 81.40.Pq; 81.40.Lm; 81.05.Ys; 61.46.+w; 61.16.Bg
1. Introduction Early studies of the tribological behavior of metallic glass used thin specimens produced by rapid quenching from the melt [1,2]. The recent availability of bulk metallic glass (BMG) [3,4]
*
Corresponding author. Tel.: +1-614 292 1775; fax: +1-614 292 1537. E-mail address: [email protected] (D.A. Rigney).
simplifies tribological testing and allows more thorough chemical and structural characterization of worn specimens [5]. The experimental work reported here was undertaken to confirm reports, e.g., in [4], of favorable friction and wear characteristics of metallic glass. To complement the experimental work, molecular dynamics (MD) simulations of the deformation and sliding of a model amorphous solid were also performed. Others have used MD techniques to study unlubricated friction involving crystalline materials
0022-3093/03/$ - see front matter 2003 Elsevier Science B.V. All rights reserved. doi:10.1016/S0022-3093(02)01999-3
X.-Y. Fu et al. / Journal of Non-Crystalline Solids 317 (2003) 206–214
[6–11]. Those earlier simulations showed a reduction of coefficient of friction (COF) with increasing sliding speed, mechanical mixing and the development of nanocrystals adjacent to the sliding interface [10].
2. Experiments The alloy chosen was Zr41:2 Ti13:8 Cu12:5 Ni10:0 Be22:5 (composition in at.%) cast by Howmet. Its hardness was HV530 with a load of 0.3 kg and the shear modulus has been reported to be 33 GPa [12]. X-ray diffraction (XRD) and high resolution transmission electron microscope (TEM) imaging confirmed that the structure was amorphous. Sliding tests were done with a simple pin/disk arrangement, shown schematically in Fig. 1. The chamber could be evacuated, so tests could be done in air or vacuum (typically 4 103 Pa). Not shown in Fig. 1 is a Kelvin probe that was used in situ to monitor the combined effects of structural and chemical changes on the surface of the disk. Test pins were 3 3 12 mm rods of BMG or 6 mm diameter balls of 52 100 bearing steel. Disks were 31 mm in diameter and were 3 mm thick. Surfaces were prepared by metallographic polishing, using 1 lm diamond paste for the final stage and ultrasonic cleaning in acetone, methanol and distilled water. XRD confirmed that the prepared surfaces remained amorphous. Normal loads ran-
207
ged from 0.1 to 1.2 kg and sliding speeds were typically 0.05 m/s but some tests were run up to 1 m/s. Average wear rate was determined from measurements of weight changes caused by sliding. Structural and chemical characterization of worn surfaces and debris was accomplished by optical microscopy, scanning electron microscopy (SEM), energy dispersive X-ray Spectroscopy (EDS), XRD, TEM and microhardness. Other details of the equipment and procedure are presented in references [5] and [13].
3. MD simulations No attempt was made to model in three dimensions the five-component BMG material used in the experimental part of this study. Such a simulation would require interatomic potentials that are not available and would involve excessive computation time. Instead, a Lennard–Jones potential and a simple 2D amorphous system was used. It was anticipated that these simulations would exhibit some of the major features that a more realistic model would generate. The energy and length scale parameters in the potential were those reported in [14]. The composition was chosen such that the ratio of the number of L or ÔlargeÕ atoms to S or ÔsmallÕ atoms was ð1 þ 51=2 Þ=4, one half the golden mean, because this ratio is characteristic of a Penrose tiling of space by rhombic
Fig. 1. Experimental set-up for pin/disk tests.
208
X.-Y. Fu et al. / Journal of Non-Crystalline Solids 317 (2003) 206–214
units consisting of these particles (atoms) [14]. This choice optimizes for the existence of a quasicrystalline ground state that corresponds to this tiling. The chosen system is easily quenched to a metastable glass having a glass transition temperature Tg of 0.3 (rescaled dimensionless units) [15]. All simulation values reported in this paper are presented in these rescaled units. Details of the potential and the rescaled units are given in Appendix A and in [13–16]. The model system contained 4472 atoms designated L and 5528 atoms designated S, for a total of 10 000 atoms. These were contained in a square simulation cell with periodic boundary conditions. Temperature and pressure control were accomplished by using a Nose–Hoover thermostat and a Parinello–Rahman barostat [17]. After the system was equilibrated at high temperature (T ¼ 2:3) and pressure, the system was quenched below Tg . The quench sequence assured that voids were suppressed and that the product was a metastable 2D glass at a very low temperature (T ¼ 0:01) and zero pressure. The mechanical properties of the 10 000 atom model system were checked by running two kinds of tensile test simulations, one with periodic cell boundaries on all sides in order to eliminate surface effects and the other with this constraint relaxed along the edges parallel to the tensile axis. In both cases, zero stress boundary conditions were maintained perpendicular to the tensile axis, resulting in comparable uniaxial tests. For simulations of sliding, an identical set of 10 000 atoms was placed in contact with the first set, but shifted along the x axis, i.e., parallel to the interface by half of the cell dimension. This process generated a sliding pair consisting of 20 000 atoms. Periodic boundary conditions were retained at the cell boundaries normal to x but those conditions were relaxed at the other two boundaries, those normal to y. Also, an important slab of ÔreservoirÕ material was defined adjacent to each of those boundaries. Each of these reservoir regions was maintained at a low temperature of T ¼ 0:01, allowing frictional heat to conduct to these thermostat regions. Also, external forces were only applied to the atoms in those regions, thus simulating experiments in which the forces responsible for sliding
are imposed at some distance from the interface. To initiate sliding (t ¼ 0), all atoms in the upper block were given a velocity to the right ðþxÞ and all atoms in the lower block were given the same velocity to the left ðxÞ, providing a net interface sliding speed of vs . At t > 0, temperature and friction force control operated in the reservoirs. Force control involved a feedback loop so that sliding speed was maintained at vs . To provide a normal load, a pressure of 1.0 was applied in the y direction. Three aspects of this procedure are similar to those used earlier by Hammerberg et al. [7–9,11]: the use of two reservoirs, the periodic boundary conditions in the x-direction and the way in which sliding was initiated. The following aspects differed from those reported in [7–9,11]: the system composition and structure, the interatomic potentials, the method used to apply the shear force and the use of a constant pressure constraint rather than constant volume during sliding.
4. Results 4.1. Experimental results Sliding behavior of the metallic glass (MG) depended on load, speed and atmosphere. In slow speed tests, typically at 0.05 m/s, the steady state coefficient of friction (COF) in air for self-mated MG sliding decreased smoothly from 0:65 0:05 to 0:38 0:05 as the load increased from 0.2 to 1.1 kg. For the same load range, but in vacuum, the COF was somewhat higher, 0:82 0:05 to 0:60 0:05, but followed a similar trend. When the pin was changed to 52 100 steel, the values of COF in air and vacuum and the trends with load were similar to those found with MG pins, except that COF for 52 100 on MG in air was higher than for MG on MG at light loads. For comparison, a limited number of tests of 52 100 on 52 100 were done in air and vacuum. The COF values were close to the values found for tests with MG. In a series of tests at 0.5 kg normal load and sliding speeds ranging from 0.05 to 1.0 m/s, the COF decreased from 0:50 0:05 to 0:29 0:05. The highest wear rates were for MG on MG and 52 100 on MG in air. They ranged from 15 4
X.-Y. Fu et al. / Journal of Non-Crystalline Solids 317 (2003) 206–214
to 120 10 (expressed in 1010 kg/m) as the load varied from 0.2 kg to 1.1 kg. Wear rates in air with 52 100 as one or both components varied similarly with load, but were lower. The lowest wear rates were for MG on MG in vacuum (8 5 to 5 3, both times 1010 kg/m, for loads of 0.2–0.5 kg). SEM images and EDS maps of wear tracks after tests in air showed patches of material with enhanced oxygen content. These were consistent with Kelvin probe data obtained during sliding [18]. The groove and ridge topography was typical of worn surfaces of ductile materials. On most of the wear track, hardness had increased to HV580, but the regions with higher oxygen content were harder (HV640). Oxygen rich patches were absent after tests in vacuum and hardness values on the wear track were similar to those of the unworn material. Worn surfaces of pins were similar to those of disks. When the pin was 52 100 steel, transfer of MG from the disk was readily detected by SEM/EDS. The amount transferred was larger at higher loads or in vacuum. Debris particles were irregular flakes and plates ranging in size from 5 to 100 lm. Smaller particles with enhanced oxygen content were present after tests in air. Some cutting chips were also found, especially after tests in air. Sliding tests were also done on devitrified material obtained by annealing for 2 h at 500 C. The hardness after annealing was HV720 and XRD confirmed that the material was crystalline. Steadystate COFs and wear rates were similar to those obtained with amorphous material and XRD of the debris showed only broad peaks typical of amorphous material. SEM/EDS observations of longitudinal and transverse cross-sections (parallel and perpendicular to the sliding direction) revealed that sliding caused changes well below the surface, typically tens of microns. For specimens tested in air, surface patches with enhanced oxygen content were detected by back-scattered electron images and by EDS. The hardness in these regions was also enhanced (600–700 Knoop hardness, KHN). For specimens tested in vacuum, the material adjacent to the surface had decreased hardness compared with bulk material. That is, the material showed work-softening rather than work-hardening. The
209
softer region extended for 30–70 lm from the surface. After tests in air, a corresponding soft layer was also detected between the locally hard patches and the underlying base material. No differences in appearance or in chemical composition were detected between the softer material and the base material. 4.2. Results of MD simulations Fig. 2(a) shows the results of simulated tensile tests on the model MG system for a range of strain rates. For these simulations, periodic boundaries were maintained along the surfaces parallel to the tensile axis. The response is typical of an elasticperfectly-plastic material. The transition to plastic flow occurs at higher values of stress as strain rate is increased. The YoungÕs modulus varies from 28.4 to 32.9 (rescaled units) and the 2D PoissonÕs ratio from 0.55 to 0.50 as the strain rate is increased from 1 105 to 2 104 . Fig. 2(b) (lower strain rate) and Fig. 2(c) (higher strain rate) show that voids developed during the tensile test simulations. These were smaller and more randomly distributed for the high strain rate case. For the larger voids, it is clear that they nucleated where atoms of the same kind were clustered in the original material. Fig. 3(a) shows the results of a tensile simulation in which the periodic constraints at the cell walls parallel to the tensile axis were removed. This allowed these surfaces to roughen. In this case, the shading indicates regions in which strain is concentrated. The corresponding stress-strain curve has the same general form as those shown in Fig. 2(a), but with larger excursions from the general trend. Fig. 3(b) uses the same shading convention, but for the simulation already shown in Fig. 2(b). In that case, in addition to seeing the voids, one can see regions with localized strain, but they are not as strongly developed as in Fig. 3(a). Simulations of sliding of MG against MG showed COF rising quickly to a maximum value of 0:50 0:05 and then falling to a steady-state value that decreased from 0:31 0:05 to 0:06 0:03 as sliding speed increased from 0.1 to 1.0 (rescaled units). Fig. 4 shows a sequence in the interface region for a sliding speed of 1.0. The material above the
210
X.-Y. Fu et al. / Journal of Non-Crystalline Solids 317 (2003) 206–214
Fig. 3. (a) Structure produced by tensile strain at slower strain rate (1 105 , rescaled units), showing shear bands when surface roughening is allowed. Shaded regions have higher strain. (b) Same as 3a, but periodic boundary conditions imposed at sides of cell, showing voids and incipient shear bands.
Fig. 2. (a) Tensile stress–strain curves at different strain rates (1 105 to 2 104 , rescaled units) from MD simulations. Periodic boundary conditions imposed at sides of cell. (b) Structure produced by tensile strain at slower strain rate (1 105 , rescaled units), showing large voids. (c) Similar to 2b, but at high strain rate (2 104 , rescaled units), showing smaller voids.
interface is color coded differently from the material below the interface but is otherwise the same. The initially smooth interface first roughens and then a mixing process is observed. A movie version of this sequence shows that mechanical mixing
near the interface proceeds by a random process involving rotations of irregular clusters of atoms. That is, it is not a simple atomic diffusion process. Closer to the reservoir regions, the material flow approximates laminar flow parallel to the sliding direction. The color coding allows one to monitor the progress of mixing, and, in particular, the growth kinetics of the mixed layer. This was done by calculating the net displacements along the y direction (normal to the interface) relative to the center of mass for atoms in thin slabs parallel to the sliding direction. Each resulting profile was fit to a Gaussian function centered at the interface, allowing an objective determination of the thickness of the mechanically mixed material at different times. As shown in Fig. 6, the thickness of this material was proportional to t1=2 . A plot of the
X.-Y. Fu et al. / Journal of Non-Crystalline Solids 317 (2003) 206–214
211
Fig. 4. Sequence showing interface roughening and mixing during sliding for t ¼ 0 and two later times (rescaled).
strain rate vs. distance from the interface can also be fit to a Gaussian function centered at the interface, and the thickness of the high strain rate region is also proportional to t1=2 . The transients observed in simulations of COF suggest that changes in the near-surface material occur during sliding. The simplest change to investigate would be changes in density. Therefore, the density in thin layers parallel to the interface was compared with the original density. The density difference profile is shown in Fig. 5. This shows a clear decrease in density in the most highly deformed material. Part of this density decrease is Fig. 6. Density difference profile relative to average density.
Fig. 5. Growth kinetics of the thickness of the mixed layer, from MD simulations.
from thermal expansion, but when this smaller contribution is subtracted, the decrease in density near the interface remains. The slight increase in density elsewhere in the cell arises because of the use of constant volume conditions. The highest local temperature at t ¼ 25, determined from the random part of the kinetic energy in each thin layer, was 0.1, well below the Tg of 0.3. In experimental work on crystalline materials, various features associated with the microstructure can be used as markers to measure displacement profiles adjacent to sliding interfaces. These are not available in MG, but color coding provides convenient markers in MD simulations. Fig. 7 shows how an initially straight marker bends over
212
X.-Y. Fu et al. / Journal of Non-Crystalline Solids 317 (2003) 206–214
Fig. 8. Displacement profiles at four different times, from MD simulations.
Fig. 7. Simulation of bending of an originally straight marker during sliding.
chosen time. The results show that profiles from MD simulations and from Navier–Stokes analysis are in agreement [16].
5. Discussion in response to the applied shear stresses for sliding of MG on MG. This profile is similar to displacement profiles observed for a range of size scales in experiments on various single and multiphase crystalline materials. In fact, they all resemble the profiles expected for the case of a viscous fluid near a moving boundary. This would be described by the Navier–Stokes equation with suitable boundary conditions. The familiar solution for the fluid velocity as a function of distance from the interface is a complementary error function and this yields a similar displacement profile. If one accounts for the elastic to plastic transition in the MG, one can generate a time sequence of displacement profiles, as shown in Fig. 8. After an initial transient, the elastic strain saturates and the plastically deformed zone reaches a steady state depth of about 50 length units. If one subtracts the elastic contributions to these curves, one can fit the simulation data for a chosen time to the profile predicted by Navier–Stokes analysis. This gives a value for the effective viscosity (assuming it remains constant), and that parameter in turn allows prediction of the displacement profile for any
The unlubricated friction and wear results reported here suggest that bulk metallic glass has tribological behavior similar to that of more familiar crystalline materials. In both cases the loading conditions during sliding are compression and shear; these are conditions that inhibit fracture and encourage ductility. The results may be quite general for systems in which plastic deformation continues during sliding. Blau has recently reported similar results and conclusions for a Zr–Cu–Ni–Ti– Al BMG [19]. Earlier expectations that BMGÕs might have exceptional sliding characteristics were based largely on experiments with thin ribbons of Fe-based MG [2,20]. It may be that MG materials will exhibit a range of sliding behavior, just as crystalline materials do, depending on composition and the presence or absence of other phases. It is clear that some oxidation occurs during sliding in air, but this does not produce an oxide layer that yields oxide wear debris directly. Instead, the oxide that forms is blended into the deforming MG, and it is this mixed material that generates the oxide content in wear debris. Different amounts of
X.-Y. Fu et al. / Journal of Non-Crystalline Solids 317 (2003) 206–214
mixed in oxide cause differences in local hardness on the wear surface and in transfer material on the pin. If the local hardness difference at contacting asperities becomes greater than 20%, cutting debris can be generated [21]. No evidence of sliding-induced devitrification was found. In fact, the data from tests on annealed BMG showed the opposite effect: the deformation produced by sliding reamorphized the material. This could be a useful characteristic of BMG – if any local region did crystallize, perhaps in response to temporary thermal excursions, continuing sliding would restore the amorphous structure. The soft layer created by sliding is very interesting for at least two reasons. First, it assures that, in the absence of outside influences that can raise the surface hardness, sliding on BMG produces and maintains a soft surface layer on a harder base. This controls friction and correlates with relatively smooth sliding. Second, the data are consistent with suggestions that increases in free volume are associated with plastic deformation of MG [21,22]. The regions of decreased density predicted by the MD simulations would also be consistent with this view. The shape of the tensile stress–strain curves shown in Fig. 2(a) is expected for a material that lacks the microstructural features that would allow work-hardening. The smoother curves at higher strain rates would be expected for a structure that is more homogeneous, i.e., having smaller and more evenly distributed voids. The voids shown in Fig. 2(b) are larger because they have more time to grow at slower strain rates. The nucleation of the voids is affected by local fluctuations in composition. They occur more easily where there is a local concentration of weaker bonds. In this simulation, bonds between two L atoms or two S atoms are weaker than those between unlike atoms. The results of tensile test simulations depend also on conditions at the surface. When the surface is allowed to roughen, shear bands develop and void formation is inhibited. The corresponding stressstrain curve is consistent with the appearance of shear bands. A transient peak followed by a decrease to a steady state value are common features of experimental COF vs. time plots. The MD results predict
213
similar behavior. Experimental measurements show a decrease in COF as sliding speed increases. A similar pattern results from the MD simulations. MD simulations with crystalline copper show similar Ôvelocity softeningÕ [10]. Even more striking are the observations that mechanical mixing occurs in sliding experiments and in MD simulations for both MG and crystalline material [10]. In the simulations on MG, this mixing is associated with nanoscale rotational flow in material where the strain rates are highest. In this region, the collective effect of continuing shear-induced microstructural rearrangements in the amorphous material can lead to long-range mass transport normal to the principal shear direction. Rotational flow and mixing may also be responsible for the reamorphizing observed for sliding of devitrified MG.
6. Conclusions This project began with an emphasis on the sliding friction and wear behavior of bulk metallic glass. The experimental observations, combined with the results of molecular dynamics simulations, suggest that the flow processes involved are more general than originally thought. It is suggested here that amorphous materials are well suited to serve as model materials, both for experiments and simulations, allowing identification of generic properties and behavior in a broader range of materials and processes. The insight obtained from these studies may be helpful for understanding the behavior of both amorphous and crystalline materials when large plastic strains are imposed at high strain rates. Potential areas of application include friction and wear, high strain rate deformation, mechanical alloying and friction stir welding. Earlier reports [2,20] that BMGÕs might provide exceptional friction and wear characteristics are not confirmed by the results of the present work.
Acknowledgements We are pleased to acknowledge research support from the National Science Foundation
214
X.-Y. Fu et al. / Journal of Non-Crystalline Solids 317 (2003) 206–214
(NSF), Surface Engineering and Tribology, most recently under Grant no. CMS-9812854, NSF, Division of Materials Research under Grant no. DMR-0135009, and from Howmet Corp. We are grateful to Sharon Glotzer, then at NIST), The Ohio Supercomputer Center and the University of MichiganÕs Center for Parallel Computing for the use of computer facilities, and we thank the following for helpful discussions: M.J. Mills, P.K. Gupta and A. Markworth at The Ohio State University, A. Zharin at the Belarussian State Powder Metallurgy Concern, Minsk, Belarus, and J. Hammerberg at Los Alamos National Laboratory. Appendix A
If the system simulated were pure Cu, the relative ; e ¼ 0:35279 coefficients would be r ¼ 2:5172 A 25 kg. eV; m ¼ 1:0551 10 Based on these values, the following relationship can be calculated ; r^ ¼ 1 () 2:5172 A
^t ¼ 1 () 3:4423 1013 s; P^ ¼ 1 () 3:5392 GPað3DÞ or P^ ¼ 1 () 0:89087N=m ð2DÞ; v^ ¼ 1 () 731:42m=s:
References
The Lennard–Jones pair potential has the form r 12 r 6 U ¼ 4e ; r r where r is the distance between two atoms; r is the distance at which the interatomic potential U is zero and e is the bond strength. The MD simulation used Lennard–Jones potentials with units rescaled such that all quantities are expressed relative to the scaling factor. After rescaling, the simulation results are not associated with any particular system; different systems would only affect the scaling factors, leaving the rescaled simulation results unchanged. The rescaling factors for some of the most commonly used quantities are given below:
[1] [2] [3] [4] [5] [6] [7] [8] [9] [10] [11] [12] [13] [14]
• BoltzmannÕs constant and each atomic mass are set equal to 1.0. • Length is rescaled such that ^ r ¼ r=r where ^ r is the rescaled length and r is the actual length. • Energy is rescaled such that e^ ¼ e=e. pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 2 =ÞÞ. • Time is rescaled such that ^t ¼ t=ð ðmr pffiffiffiffiffiffiffiffi • Velocity is rescaled such that v^ ¼ v=ð e=mÞ. • Pressure, stress and modulus are rescaled such that P^ ¼
P for 3D or ðe=r3 Þ
e^ ¼ 1 () 0:35279 eV;
P^ ¼
P for 2D: ðe=r2 Þ
[15] [16] [17] [18]
[19] [20] [21] [22]
P.G. Boswell, J. Mater. Sci. 14 (1979) 1505. K. Miyoshi, D.H. Buckley, ASLE Trans. 27 (1983) 295. W.L. Johnson, MRS Bull. 24 (1999) 42. A. Inoue, Acta Mater. 48 (2000) 279. X.-Y. Fu, T. Kasai, M.L. Falk, D.A. Rigney, Wear 250 (2001) 409. G. He, M.O. Robbins, Science 284 (1999) 1650. J.E. Hammerberg, B.L. Holian, J. Roeder, A.R. Bishop, S.J. Zhou, Phys. D 123 (1998) 330. B.L. Holian, J.E. Hammerberg, P.S. Lomdahl, J. Comp. Aided Mater. Design 5 (1998) 207. J.E. Hammerberg, B.L. Holian, S.J. Zhou, Chem. Phys. Rep. 17 (1998) 343. D.A. Rigney, J.E. Hammerberg, MRS Bull. 23 (1998) 32. J. Roeder, A.R. Bishop, B.L. Holian, J.E. Hammerberg, R.P. Mikulla, Phys. D 142 (2000) 306. W.L. Johnson, Mater. Sci. Forum 225–227 (1996) 35. X.-Y. Fu, PhD dissertation, The Ohio State University, Columbus, OH, 2001. F. Lancon, L. Billard, P. Chaudhari, Europhys. Lett. 2 (1986) 625. M. Falk, J. Langer, Phys. Rev. E 57 (1998) 7192. X.-Y. Fu, M.L. Falk, D.A. Rigney, Wear 250 (2001) 420. S. Melchionna, G. Ciccotti, B.L. Holian, Mol. Phys. 78 (1993) 533. T. Kasai, X.-Y. Fu, T. Honda, A. Zharin, D.A. Rigney, Proceedings of International Tribology Conference, Nagasaki, 2000, Japan Society of Tribologists, Tokyo, 2001, p. 333. P.J. Blau, Wear 250 (2001) 431. T. Imura et al., Mater. Sci. Eng. A 133 (1991) 332. T. Akagaki, D.A. Rigney, Wear 149 (1991) 353. D. Deng, A.S. Argon, S. Yip, Philos. Trans. R. Soc. Lond. A 329 (1989) 613.