We discuss large-scale non-equilibrium molecular dynamics (NEMD) simulations of ductile metal sliding comprising up to 1.8 × 10^{9} atoms over time scales of 100 ns. The results of these simulations have identified a variety of physical mechanisms that are important in determining the steady-state frictional force for a wide range of velocities at compressed metal–metal interfaces. These include grain growth and refinement, the evolution of large plastic strains and strain rates, material mixing, and melting. These phenomena can be included in a strain, strain rate, and grain size model that gives good agreement with the NEMD simulations and can be applied to macroscopic continua.

## I. INTRODUCTION

The physics of frictional interactions at ductile metal interfaces is complex, involving material mixing, grain structure transformation, and the formation of graded microstructures (Rigney, 1988; Rigney and Hammerberg, 1998; Rigney, 2000; Rigney and Karthikeyan, 2010; Hughes and Hansen, 2001). At high sliding speeds, the frictional force typically exhibits velocity weakening, i.e., a decrease in the frictional force with increasing velocity with melting occurring at very high velocities (Yuan *et al.,* 2009; Bowden and Freitag, 1958). Large-scale non-equilibrium molecular dynamics (NEMD) simulations provide a window into the microscopic world of frictional dissipation (Hammerberg and Holian, 2004). Early NEMD simulations for ductile metals such as Cu (Hammerberg *et al.*, 1998) for single crystal two-dimensional systems showed all of the above phenomena. Atomistic simulations studies of fully three-dimensional compressed single crystal aluminum interfaces showed that, regardless of sliding orientation, the frictional force increases as a function of sliding velocity and has a maximum, which coincides with the temperature at the interface reaching the melt temperature of the material at the appropriate pressure. Above this critical velocity, the frictional force decreases with sliding velocity following a universal power law with an exponent close to ¾ (Hammerberg *et al.*, 2017). Polycrystalline systems, while following the same behavior above a critical velocity at which the frictional force is a maximum, exhibit a more complicated evolution to steady state due primarily to grain growth, grain rotation, and grain boundary accommodation (Hammerberg *et al.*, 2018). Computational studies of such systems require much larger length scales than single crystal studies. This is primarily due to the minimum grain size required to realistically simulate plastic deformation within a grain as well as to be able to reproduce dislocation pile up processes at grain boundaries, which are absent in grain sizes below the Hall–Petch maximum. For nanocrystalline aluminum, large-scale computer simulations put the Hall–Petch maximum around 25 nm (Kadau *et al.*, 2004; and Xu *et al.*, 2018).

Due to advances in hardware architectures and massively parallel codes, it is now possible to investigate polycrystalline systems at much larger length scales with particle numbers of several billion atoms, time scales approaching 100 ns and grain sizes of the order of tenths of microns. In this work, we will discuss a series of simulations of sliding friction of polycrystalline Aluminum interfaces compressed to a nominal pressure of 15 GPa, typical of plate impact or explosively driven experiments, and for sliding velocities in the range of 10–2000 m/s. The effect of initial grain size was systematically studied by utilizing (initially) mono-dispersive grain distributions with initial mean grain sizes of 4, 13, 20, and 50 nm, which translates into atomistic systems containing between a few hundred million atoms and up to 2 × 10^{9} atoms. We describe the dynamics of the grain structure leading to a steady sliding state very different from the initial structure, characterized by grain growth and refinement at high rates of plastic deformation and large values of plastic strain. The organization of the paper is as follows: in Sec. II, we detail the simulation methodology, boundary conditions, and sample geometry. In Sec. III, we discuss some main results of our simulations, contrasting the effect of initial grain size, strain buildup leading to grain refinement, time evolution of initial grain distribution, and its effect on the measured frictional force. In Sec. IV, we present a constitutive model that quantitatively reproduces the velocity dependence of the frictional force seen in the simulations. We provide our summary and conclusions in Sec. V.

## II. SIMULATION METHODOLOGY

The NEMD simulations for metal–metal sliding have been carried out using the high-performance parallel molecular dynamics code SPaSM (Scalable Parallel Short Range Molecular dynamics) developed at the Los Alamos National Laboratory (Lomdahl *et al.*, 1993). The computational configuration uses a time-dependent boundary condition to measure the frictional force required to establish an interface sliding with a constant relative velocity between two materials. A schematic of the computational geometry is shown in Fig. 1. For discussions related to the geometry, the sliding direction is the x-direction and the vertical orientation will be taken in the y-direction. Periodic boundary conditions are applied in the lateral directions, but the top and bottom surfaces are constrained from expanding by perfectly reflecting walls, which reflect the velocity of atoms crossing the walls back into the computational box in the normal direction (y). Reservoir regions at the top and bottom of the computational box are kept at a constant temperature. The thickness of the reservoir regions, $\Delta y$, is small relative to the overall dimension of the work pieces but thick enough to contain a large number of lattice planes, typically of order 20. For instance, for computational cells with about 27 × 10^{6} atoms, the reservoir regions contain about 200 000 atoms. For computational cells of 1.8 × 10^{9} atoms discussed below, the reservoir regions contain 25.8 × 10^{6} atoms. The sliding boundary condition is enforced in these reservoir regions located at positions $(yres\xb1,yres\xb1\u2213\Delta y)$, where $yres\xb1=\xb1L/2$ parallel to the initial interface at y = 0 at the initial time t = 0, where L is the vertical dimension of the computational box. The constant velocity of sliding (v = 2u_{p}) is maintained by applying a frictional force, uniformly to each atom in the reservoir regions such that the average velocity is $\xb1up$ in the upper and lower reservoir regions.

In addition to the boundary condition with reflecting walls applied in the y-direction as described above (constant volume), we have also utilized a constant normal load boundary condition, in which a constant normal force is applied to the reservoir atoms which maintain a constant longitudinal stress. We will discuss here the results using the former fixed volume boundary condition. The initial sample is equilibrated at a temperature of 300 K and pressure of 15 GPa, although other initial temperatures were also investigated. At the beginning of the simulation, a velocity ±u_{p} is applied for the upper (y > 0) and lower regions (y < 0), respectively.

Polycrystalline samples were constructed employing a constraint Voronoi tessellation method (Weaire *et al.*, 1986) to create an initial mono-dispersive grain size distribution with a nominal grain size. The number of grains in the samples varied from N_{g} = 250 for a 50 nm initial grain size to over 30 000 grains for a 4 nm initial grain size. This initial configuration is subsequently relaxed by an annealing procedure, in which the sample is strained and annealed at 800 K for 500–800 ps and relaxed to a hydrostatic state at a pressure of 15 GPa and a temperature of 300 K. Figure 2 shows a typical grain map for a large sample with a 50 nm average grain size and comprised of 1.8 × 10^{9} atoms. The grain orientation is indicated by the color. This “orientation imaging map (OIM)” grain coloring scheme is based on the crystal orientation of the grain relative to the axis normal to the sliding plane (y axis) with red = (100), blue = (111), and green = (110) (Ravelo *et al.*, 2013). This initial grain morphology is to be contrasted with the grain distribution at later times, in which the sample goes from mono-dispersive to poly-dispersive. The atomic interactions are described by an embedded atom method (EAM) potential of aluminum (Angelo *et al*., 1995; and Ravelo *et al*., 1998). This model potential describes well the elastic moduli, gamma surface (relevant to dislocation nucleation and motion), and stress–strain curves up to 29 GPa.

## III. RESULTS OF POLYCRYSTALLINE Al–Al SLIDING SIMULATIONS

Figure 3 shows a sequence of snapshots from the 50 nm grain size simulation for an initially flat Al–Al interface with a sliding velocity of v = 2u_{p} = 60 m/s. The system size in this case was L_{x} = 236.6 nm, L_{y} = 473.6 nm, and L_{z} = 235.6 nm comprising 1.8 × 10^{9} atoms. The reservoir regions, located at ±L_{y}/2 each contained 25.8 × 10^{6} atoms. The times shown are at 10, 20, 30, 40, 50, and 60 ns. By 10 ns, the near interface region has developed very significant plastic deformation and grain distortion. This is mirrored in the velocity and temperature fields normal to the interface. These are shown in Figs. 4 and 5. In the steady state, characterized by a statistically time independent frictional force and temperature profile which occurs after approximately 40 ns, a layer of thickness D, approximately 150 nm wide forms at the sliding interface corresponding to a strain rate of approximately 2 × 10^{8} s^{−1}. The thickness D is determined from the velocity profiles as shown in Fig. 4 and corresponds to the region with a nearly linear variation in the sliding velocity (from −u_{p} to+u_{p}). For these constant volume simulations, the pressure in the steady state is approximately 16.7 GPa, for which the aluminum EAM potential melting temperature is 1200 K. Insight into the plastic deformation in this region can be obtained by quenching a configuration and monitoring a central symmetry order parameter used to identify defective (non-centro-symmetric) atoms (Kelchner *et al*., 1998). This can be seen at early times in Fig. 6, where a sequence of centro-symmetry images at t = 1.0, 3.0, and 5.0 ns show the development of a complex array of stacking faults and dislocations near the sliding surface. The upper images are OIM grain maps and the lower images are the corresponding centro-symmetry images. Figure 7 shows a snapshot taken at t = 62 ns well into the steady-state regime. Similar results have been obtained with other initial grain size distributions; in particular, 4, 13, and 20 nm grain sizes have been studied (Hammerberg *et al.* 2018; 2017). An example of a system with an initial grain size 4 nm is shown in Fig. 8 (Hammerberg *et al*., 2018) for times 0, 2.5, and 5.4 ns. The 100 × 200 × 100 nm sample with sliding velocity v = 100 m/s, initial pressure 15 GPa, and boundary temperature 300 K shows very similar grain evolution and coarsening. The polycrystalline simulations have demonstrated that the initial grain distribution evolves to a steady-state distribution very different from the initial grain distribution, showing grain growth and grain refinement in the steady state. Thus, in these simulations, there is no Hall–Petch relation based on the initial grain size.

## IV. CONSTITUTIVE MODEL OF FRICTION IN SLIDING METALLIC INTERFACES

For the particular boundary condition employed in these simulations, it is also the case that, both for polycrystalline and single crystal simulations, the velocity dependence of the frictional force depends on two critical velocities, a lower critical velocity, v_{c}, and an upper critical velocity v_{c1}. For velocities below v_{c}, integrating the equation for thermal flow in the steady state gives the approximate relationship for a symmetric interface,

where T(0) is the interface temperature, T_{0} is the boundary (reservoir) temperature, L/2 is the distance from the sliding surface to the boundary, and $\kappa \xaf$ is the average thermal conductivity between temperatures T(0) and T_{0} at pressure P,

We have found Eq. (1) to be sufficiently accurate in our simulations since the steady-state temperature profile (cf. Fig. 5) is nearly linear. If more accuracy is required, the full thermal transport equation may be solved. The critical velocity v_{c} is the velocity at which the interface is at the melting temperature and a quasi-two-dimensional molten layer has formed. For the 50 nm grain size simulations shown in Figs. 2–7, v_{c} is approximately 62 m/s. From Eq. (1), the lower critical velocity is given by

where T_{m} is the melting temperature at pressure P and f_{c} is a critical force per unit area. For $v\u2264vc$, the highly deformed region near the interface seen in Figs. 3, 6, and 8 has an average thickness D (see Fig. 4). This is a dynamical quantity that reaches a steady-state value. Our assumption is that, within this region, our sample size is large enough that we may express the critical force per unit area, f_{c}, as the transverse component of a flow stress, $\sigma t(T,P;\Psi ,\Psi \u02d9)$ determined by the local values of temperature T, pressure P, plastic strain $\psi $, and plastic strain rate $\psi \u02d9$. We also assume that $\Psi $ reaches a critical value $\Psi c$ that limits grain growth. Figure 9 shows an example of the time dependent grain growth for the largest grain for 13 and 20 nm initial grain size samples showing the transition to grain refinement events at later times.

We denote the plastic strain rate at v_{c} by λ with

where $\gamma $ is a numerical constant equal to 1/√3 for a Prandtl–Reuss flow rule and von Mises flow surface. Then,

The right-hand side of Eq. (5) denotes a rate dependent model for the transverse component of the flow stress. We have taken the model of Preston *et al*., (2003), which is accurate to high rates of plastic deformation for our analysis. For velocities greater than the critical velocity, we have found that the frictional force per unit area, f(v), scales as a power law in (v/v_{c}) and that

with $\beta $ close to the value $1\u2212\alpha $, where $\alpha $ is the scaling exponent for the high strain rate branch of the rate dependent plasticity model of Preston *et al*. (2003). A discussion of this scaling for $v\u2265vc$ is given in the Appendix. The model is then specified by

In principle, the interface temperature, T(0), may be obtained from a solution of the thermal conduction equation but we have found that a good approximation is, from Eqs. (1) and (3),

with f_{c} given by Eq. (5). When this model is compared with the NEMD data for initial grain sizes of 20 and 50 nm (Hammerberg *et al*., 2017), the agreement is within 15% of the average NEMD steady-state values and generally within the temporal statistical fluctuations of f(v). Figure 10 shows a comparison between NEMD data and the constitutive model for two different initial grain sizes. The NEMD data correspond to polycrystalline Al–Al with 20 and 50 nm initial grains and the frictional force is the measured averaged steady-state value. The values for the parameters D and $\Psi c$ used, taken from the simulations, are D = 60 nm, $\Psi c=3$ for the 20 nm grain size simulations and D = 150 nm, $\Psi c=2$ for the 50 nm grain size simulations. The values for $\Psi c$ were calculated from the strain rate by multiplying $\Psi \u02d9=\gamma vD$ by $\Delta t$, where $\Delta t$ is the temporal periodicity of coarsening–refinement transitions seen both in the grain size analysis and the frictional force time dependence, $\Psi c=13vD\Delta t$ (Hammerberg *et al*., 2017). The corresponding critical velocities, calculated from Eqs. (3)–(5), i.e.,

are v_{c} = 110 m/s for the 20 nm grain size simulations and v_{c} = 62 m/s for the 50 nm grain size simulations.

It is of interest that although the frictional force per unit area has a self-similar form as a function of v/v_{c} for v ≥ v_{c} the interface temperature does not. Rather, for v > v_{c}, there is a second critical velocity, v_{c1}, for which when v_{c} ≤ v ≤ v_{c1}, $T(0)\u2245Tm$ and when $v\u2265vc1$, $T(0)\u2265Tm$ (Hammerberg *et al*., 2014) and the sliding interface develops a fully three-dimensional fluid layer. This behavior is analogous to melting along a path in the T-P thermodynamic plane that crosses the melt line, where the system resides in a two-phase region at the melt temperature as pressure increases until the entropy change is such that the phase is fully fluid and the temperature rises above the melt temperature.

The estimates of plastic deformation that develop in the region near the interface, calculated from the strain rate assuming that the strain is almost completely plastic, are very large, of order 2–3, similar to those seen in experiments for severe plastic deformation (Orlov *et al*., 2012; and Terhune *et al*., 2002). The centrosymmetry image in Fig. 7, for a sliding velocity of 60 m/s and initial pressure of 15 GPa is a demonstration of the very large plastic deformation in the steady sliding state.

## V. DISCUSSION

We have described a series of large scale NEMD simulations for compressed Al–Al interfaces for initial grain sizes from 4 to 50 nm, the largest of which, for 50 nm polycrystalline Al–Al interfaces at initial pressures of 15 GPa have been for 1.8 × 10^{9} atoms.

These have shown the following features:

The initial grain structure evolves to a larger scale steady-state grain structure characterized by large plastic strains (of order 2 and above) in the steady sliding state.

In the steady sliding state, both grain growth and refinement occur.

The velocity dependence of the average frictional force per unit area is characterized by two critical velocities v

_{c}and v_{c1}with v_{c}< v_{c1}.

For v ≤ v_{c}, the interface temperature is below the melt temperature. For v_{c} ≤ v ≤ v_{c1}, the interface temperature is near the melting temperature with a quasi-two-dimensional fluid layer. For $v\u2265vc1$, the interface temperature is above the melting temperature and a three-dimensional fluid Couette flow region is formed with a thickness, $d=v/\epsilon \u02d9c$, characterized by a critical strain rate $\epsilon \u02d9c$ and a temperature increase that is supra-linear in velocity. The results of these simulations can be described by a constitutive model that takes into account the formation of a transformation layer at the sliding interface characterized by very large plastic strains and high plastic strain rates. It assumes that in this region the simulation size is large enough that continuum plasticity models may be applied. For the very large sample sizes reported here, we believe this to be valid. The rate-dependent plasticity model that we incorporate is the model of Preston *et al*. (2003) which is able to treat the high rates of deformation seen in these simulations. The model for the frictional force requires thermal conductivities, melting temperatures and subsumes the rate-dependent plasticity model. It treats the transition from solid to fluid interfacial behavior by means of a second critical velocity. This second critical velocity enters into the calculation of the fluid interface temperature and thickness but not into the model for the frictional force where only v_{c} appears [Eq. (7)].

In large-scale NEMD simulations, the thermal conductivity is due only to lattice and defect scattering, not to electron–electron and electron–phonon scattering. In applications, the total thermal conductivity would enter into the definition of v_{c}. The model has been applied to simulations with a particular boundary condition at ±L/2, where the temperature is held fixed at T_{0}. In a continuum code implementation L becomes a time dependent quantity L(t) and the boundary temperature also time dependent T_{0}(L(t)/2) with L(t)/2 determined through a solution of the thermal transport equation. The model requires the thickness, D, of the strongly plastically deformed layer in the steady sliding state. A more complete theory requires a micro- and meso-scale theory of the grain dynamics, including material mixing, grain growth and refinement and dislocation evolution at large strains and strain rates that would allow this quantity to be calculated.

## ACKNOWLEDGMENT

This work was supported by the U.S. Department of Energy (DOE) through the Los Alamos National Laboratory. Los Alamos National Laboratory is operated by Triad National Security, LLC, for the National Nuclear Security Administration of U.S. Department of Energy (Contract No. 89233218CNA000001). The support of the LANL Advanced Technology Computing Campaign and the DOE ASC-PEM program is gratefully acknowledged. It is a pleasure to acknowledge the advice and support of Dr. B. L. Holian during the course of this work.

## AUTHOR DECLARATIONS

### Conflict of Interest

The authors have no conflicts to disclose.

## DATA AVAILABILITY

The data that support the findings of this study are available from the corresponding author upon reasonable request.

### APPENDIX: FRICTIONAL FORCE: HIGH-VELOCITY SCALING

A short derivation and justification of the high velocity ($v\u2265vc$) behavior given in Eq. (7) follows from considering a volume $\Omega $ at the sliding interface (Hammerberg *et al*., 2002). We assume that, in this volume, the frictional work per unit time is dissipated by plastic deformation,

where $\sigma $ is the flow stress and $\psi \u02d9$ is the time derivative of the plastic strain. If we assume that the interface undergoes a phase transformation when a critical plastic strain, $\psi c$, is achieved then,

where A is the surface area and *l* is the normal dimension to the interface. For very large dislocation densities, the thickness of the region, $\Omega (t)$, is assumed to be diffusive with

where $D$ is a dislocation diffusion coefficient. Then,

or

with $vc=\psi cDlc$. For very high strain rates, the model of Preston *et al*. (2003) scales as $\psi \u02d9\alpha $ giving a velocity scaling for f, from Eq. (A5), of $(vvc)\u2212(1\u2212\alpha )$.

## REFERENCES

_{3}Al