1 Introduction

Metal interfaces in relative motion are ubiquitous of many relevant machine components, such as bearings, piston rings, and gears. The mated metallic surfaces are usually only in contact at the level of individual roughness features, assuming a focus on dry or boundary-lubricated situations [1, 2]. The sliding surface is therefore repeatedly exposed to non-homogeneous pressure distributions, with fluctuating stress conditions applying near the contact regions as the asperities of the counterbody pass [3]. The applied stresses can lead to irreversible processes that alter the material’s microstructure in the near-surface zone significantly, ultimately resulting in a reduced or increased grain size depending on the loading conditions and the initial microstructure [4,5,6]. These microstructural changes have as a consequence that the mechanical properties of these near-surface zones such as toughness or hardness are altered, which can change the material’s reaction to external loads and thus directly affect its friction and wear performance [4, 7,8,9]. Tailoring materials and their interacting surfaces to meet specific requirements during their service life is becoming increasingly important to improve the performance or extend the lifetime of machine parts. A comprehensive description of the microstructural changes that occur in the near-surface of metals cannot be described using macroscopic bulk material properties, although their local properties will greatly influence the friction behavior between two surfaces [10,11,12]. In general, for conventional coarse-grained metals, sliding leads to grain refinement and the formation of a well-defined nanocrystalline layer [13,14,15]. Grain coarsening has also been experimentally observed but in connection with ultrafine-grained metals [10]. The reason for the clear and distinct separation between deformed and undeformed regions in sliding metals is still the subject of ongoing research. Single-pass experiments performed on pure copper have shown the formation of a sharp and narrow band of dislocations, dubbed dislocation trace line [16]. The presence of dislocation trace lines has been recently observed independently of the applied load and size of the counterbody [17]. For high applied contact pressures, grain refinement is initiated in the region between the dislocation trace line and the sliding interface. This observation supports the thesis that the formation of an array of dislocations might be the first step that determines the extent of the nanocrystalline layer formed after a higher number of sliding cycles. Surprisingly, the authors found that the use of identical contact pressures applied via the combination of different loads and diameter of the spherical indenters led to radically different microstructural features. The latter is crucial, as contact models relying on conventional continuum approaches may be unable to capture the complexity of the dislocation phenomena occurring beneath a sliding asperity.

The analysis through modeling of the highly dynamic processes occurring in the near-surface zones of a tribologically loaded system requires a simulation approach that can handle the occurring non-homogeneous loads to describe the deformation behavior in the first couple of nanometers from the sliding interface. Progress in high-performance computing has made molecular dynamics (MD) simulations and other meshless simulation methods viable tools for studying processes occurring in sliding or abrasive contacts [18, 19]. Mesoscopic approaches such as smooth particle hydrodynamics or the material point method are computationally cheap and have already been successfully used to study tribological aspects of scratching [20], milling [21], and grinding [22]. However, their strengths lie mainly in the stable treatment of large deformations, but they fall short of being able to resolve microstructural developments near the mated surfaces. Conversely, it has already been demonstrated in the early 2000s that plastic deformation processes like dislocation generation in nanocrystals can be modeled using polycrystalline MD [23,24,25]. One important advantage of the atomistic approach using accurate interatomic potentials is that it does not require parameters such as activation energies for structural changes like grain growth, or the assumption of constitutive material laws as in continuum mechanics [26]. Polycrystalline MD models featuring tens of millions of atoms have now matured to the level of being able to make predictions about the outcome of nanoscopic sliding processes in terms of the microstructural evolution during plastic deformation [27, 28]. These predictions can be qualitatively translated to the microscale because the simulated grains are sufficiently large to correctly reproduce a realistic material response, thus generating results of surface engineering relevance [29, 30].

Changes of the microstructure, especially subtle ones, already occur during running-in and sometimes even after a single sliding pass of the counterbody [17, 31, 32]. Since it is known that running-in greatly determines the subsequent tribological behavior of the material, it is highly important to understand the microstructural evolution happening during this stage. In particular, for subtle changes in the microstructure it is beneficial to look at the changes in crystallographic orientation/microstructure with respect to the initial state rather than the typical orientation imaging, which is experimentally obtained via electron backscatter diffraction (EBSD) and shows the current state of crystallographic orientation. Prasad et al. have impressively demonstrated such a technique to visualize microstructural changes happening in single crystal nickel subjected to sliding [33]. However, this type of visualization is usually almost impossible to carry out experimentally, as the exact initial microstructure cannot be known for a polycrystalline material after the tribological experiment. Using computer simulations, however, the same initial microstructure can be used for varying conditions and material compositions, making it easy to compare the current microstructure of a tribologically loaded surface to its initial state via exact superposition. Additionally, experimental studies at this length scale are costly and time-consuming, which is why computer simulations have become an important means to study and quantify complex phenomena as well as to identify the associated mechanisms [34, 35]. Finally, the time span observable with MD is perfectly reasonable to study the phenomena happening during the early stages of the sliding process, which greatly influences its later stages.

As this article and the entire special issue are dedicated to the memory of Mark Robbins, this may be the right place to relate a brief appropriate story. I (Stefan J. Eder) was discussing some of the groundwork for this article [32] with him standing in front of my poster at the Gordon Research Conference 2018 in Lewiston, ME, both with a cold New England pale ale in our hands. As his poster was right next to mine, we had each just spent 15 minutes “walking” the other through our work. He had many helpful suggestions, several of which made it into my subsequent work. He mentioned David Rigney, who he thought might get a kick out of this kind of simulation and visualization approach, and suggested that I should just look him up, write him a note and attach the above paper (which had just been accepted). Upon my return to Vienna I wrote to Dave, who answered right away although it was the 4th of July and he was about to turn 80 that year. Two days later, he sent me an annotated list of all his publications with little notes below the papers that I might find useful, which absolutely made my day. Although Mark was a busy man, he never kept me waiting for an answer. Two days before he passed away, I invited him to contribute to a special issue that I am guest editing for an open-access journal. Needless to say, he was the first person to reply, kindly declining because he said he preferred to publish in more classical journals, such as Tribology Letters.

In this contribution, we will first provide a concise overview of our computational model and simulation setup, followed by a description of how to produce, read, and quantify differential computational orientation tomographs (“dCOT”). We then showcase the capabilities of the analysis approach at the example of 8 representative CuNi systems, varying the alloy composition, the contact pressure, as well as the temperature to sample the relevant regions of deformation mechanism map introduced in [30] and extended in [36]. The last part of the paper focuses on the interpretation of the time- and depth-resolved maps of the dCOT color saturation and its standard deviation, allowing us to quantify plastic deformation and its lateral inhomogeneities. Wherever relevant, we attempt to make clear the limitations of the molecular dynamics approach, especially with respect to size and shear rate, as well as size-related thermodynamics issues.

2 Simulation, Visualization, and Analysis

All our MD simulations were performed using the open-source code LAMMPS [37], the de-facto standard for meshless simulations in academic tribological research. The polycrystalline model of the CuNi base body (see Fig. 1) measures \(85\times 85\times 40\) nm\(^3\) and contains approximately 25 million atoms. The grain structure was prepared using two 3D Voronoi tessellations, the main one with a mean grain size of 40 nm, and one with 10 nm grains acting as a sacrificial “handle” at the lower end to impose boundary constraints without overly interfering with the microstructural evolution in the near-surface region. Further details on sample preparation can be found in Refs. [32, 38]. While until recently the size restrictions on MD systems meant that only nanocrystalline deformation behavior could be studied, our chosen grain size is comfortably located on the macroscopic side of the Hall–Petch breakdown, so that our simulation results stand a reasonable chance of being transferable to realistic tribosystems. The FCC CuNi samples (Ni content: 0, 5, 25, 60, and 100at%) were created based on a single “master Cu microstructure” via random swapping of Cu with Ni atoms, and they interact according to an embedded atom potential with parameters taken from [39]. The system was subsequently subjected to heat treatment (650 K for systems containing Cu, 1050 K for pure Ni) and cooling to 300 K/600 K as described in [30, 36]. Rigid boundary constraints were applied to the lower 3 Å of the model, and periodic boundary conditions were applied in the lateral (x and y) directions.

Fig. 1
figure 1

Exploded 3D view of the MD system consisting of a plastically deformed FCC CuNi base body and a rigid rough counterbody. The colors represent different grain orientations, and grain and twin boundaries have been emphasized using an edge filter

Full size image

The counterbody is a rough BCC Fe monocrystal with a Gaussian root-mean square surface roughness of 0.5 nm, a fractal dimension of 2.186, and a characteristic lateral asperity extent of 33 nm, values that are similar to what has previously been applied [40], in an attempt to keep asperity flank angles low. We reduced the counterbody thickness to several monolayers and kept it rigid to have most of the computational resources available for the polycrystalline CuNi sample. The fact that the atoms in the counterbody are constrained to be fixed with respect to each other simulates a much harder material than that of the base body. Lennard-Jones potentials controlled the interactions between the counterbody and the sample, which implies a third body, with the global energy parameter \(\varepsilon =0.095\) eV obtained as described in [41], while \(\sigma _\mathrm {Fe-Cu}=0.224\) nm and \(\sigma _\mathrm {Fe-Ni}=0.221\) nm were calculated from interatomic first-neighbor distances using the Lorentz–Berthelot mixing rules.

The counterbody was moved across the surface of the sample at a sliding velocity of \(v_x=\)80 m/s and a small orthogonal component of \(v_y=9\) m/s to prevent roughness features from coming into repeated contact with their own sliding grooves. This sliding velocity may be some two orders of magnitude above what would be expected in a tribological (as opposed to a surface processing) application. However, to be able to reasonably scan the parameter space along the composition, pressure, and temperature dimensions, this is a necessary condition to ensure the availability of sufficient data. In a yet unpublished study featuring a single alloy composition at one load and temperature, but with the sliding speed ranging from 10 m/s to 1280 m/s, we found that the microstructure behaves comparably up to speeds of approximately 100 m/s, beyond which the quality of the plastic deformation changes drastically. During the simulation runs, the total applied force in \(-z\) direction was kept constant during sliding up to a total sliding distance of 400 nm. For ease of comparison, we present our results in terms of the mean normal pressure \(\sigma _{z}\) (defined as the constant force divided by the lateral cross-section of the periodic simulation box), being aware that the local stresses will vary due to the time evolution of the contact area and the subsequently increasing conformity of contact. While stress amplification due to a rigid counterbody might promote a more localized deformation pathway than in a real system, our fixed applied force condition (as opposed to a constant displacement between base body and counterbody) should provide some compliance to the rigid layer that constitutes a reasonable analog for atomic-level compliance. The normal pressure was adapted to the composition, the system temperature, and the desired deformation mechanism regime as discussed in depth in [30, 36], focusing on parameter sets where (partial) lattice rotation, twinning, and grain refinement dominate.

A Langevin thermostat with a time constant of 0.5 ps acted on all the non-rigid atoms to drain frictional heat from the system. This is equivalent to the phonons in the system being coupled to the electrons, which act as a heat bath to mimic the electronic contribution to the thermal conductivity in a metal [42] in an attempt to approximate the macroscopic heat conductivity. The thermostat acts only perpendicular to the directions of sliding and normal pressure in order not to interfere with these external constraints. In his foundational studies on shearing liquids [43], Mark Robbins observed that if the thermostat acted on degrees of freedom in the direction of shear, artifacts such as shear localization were enhanced. These artifacts are minimized when the thermostat operates on degrees of freedom not directly coupled to the shear force. Mark referred to this as a profile unbiased thermostat. It should be kept in mind, however, that without the possibility of having the lower (bulk) end of the system heat up as the unlubricated sliding simulation progresses, the temperature gradients can be kept realistic, but the total temperature will likely be too low. This is why we performed separate simulation runs at higher thermostat target temperatures to also obtain results for temperature ranges that cannot be achieved self-consistently.

Many of the computational tomographs of the CuNi alloy system discussed throughout this work are colored according to grain orientations, analogous to electron backscatter diffraction (EBSD), using the inverse pole figure (IPF) coloring standard. The orientations were calculated using polyhedral template matching [44] as implemented in OVITO [45], and the color rendition was carried out using the MTEX toolbox [46, 47] for Matlab, with the orientations projected onto the axis parallel to the normal loading direction (z). We will refer to such visualizations of our simulation data as “computational orientation tomography” (COT), being aware that, as in EBSD, the color code representing the orientations is not unique but still depends on the axis that the orientations are projected onto.

In this work, we discuss some of the more subtle microstructural changes in the system that occur during the onset of plasticity, especially migrating grain boundaries and slight lattice rotations. These may be visualized and analyzed using differential COT (referred to as “dCOT” henceforth), which compares the structure of a given time step with the initial one. In Ref. [48], this approach was introduced in a rather pedestrian fashion using an image processor, which allowed the analysis of only a handful of frames.

Here, we produced some code that can quickly produce a time-resolved set of dCOT images from a set of COT images, bringing out the full benefit of dCOT analysis by allowing the production of movies. We initially filtered out the gray representation of the rigid counterbody from the initial time step which otherwise would lead to confusing artifacts. Then, the differences between the initial and the current time step are calculated as

$$\begin{aligned} A_\mathrm {diff}(t) = 255 - |A_0 - A(t)| \end{aligned}$$
(1)

for all three color channels, where \(A_0\), A(t), and \(A_\mathrm {diff}(t)\) are the initial, the current, and the differential RGB values (0–255), respectively. As can be seen, this filter is symmetrical with respect to \(A_0\) and A(t). Care must be taken that the data types of \(A_0\) and A(t) allow negative values (int16), as most software will by default import RGB data as an unsigned integer (uint8), which leads to saturation artifacts in the difference term. To be able to export the dCOT image, \(A_\mathrm {diff}(t)\) has to be reconverted to uint8.

Fig. 2
figure 2

On how to produce and read a dCOT image. The left column shows the initial COT\(_0\) image and a subsequent one at time t together with the dCOT image at time t. The right side shows the color legend, which gives a sense of the non-uniqueness of the dCOT color scheme. The dCOT colors in the trapezoidal matrix result from applying Eq. (1) to the color bars aligned with the axes. The horizontal color bar corresponds to the path followed by the white arrow in the COT (EBSD-IPF style) orientation map, spiraling inwards from low-index to higher-index orientations (Color figure online)

Full size image

The resulting representation of the system is nearly white in all regions where no changes occurred and colorful in regions that changed by a large amount, with pastel hues where changes were slight, see Fig. 2. It is therefore perfectly suited for identifying small deviations in the courses of grain boundaries and especially slight lattice rotations. Once noteworthy changes have been identified using dCOT, the respective COT images should be revisited for a thorough discussion of the results, as the differential image may obscure which feature belongs to which configuration, and the dCOT color scheme is not unambiguous per definition. This can be seen in the bottom right of Fig. 2, where the resulting dCOT colors are shown as a trapezoidal matrix as a function of the initial and final COT colors aligned with the matrix axes.

It makes sense to compare our dCOT analysis with a small selection of existing approaches to highlight deformation in grains, most of which were developed from an experimental point of view. Quite common is the kernel average misorientation (KAM), a measure of local grain misorientation that is easily derived from EBSD data and can reveal local dislocation structures of a deformed material [49]. As the local misorientations are also related to the density of geometrically necessary dislocations (GND), KAM has also successfully been exploited for calculating the GND density using experimental EBSD data [50]. Without proper pre-conditioning, KAM maps are usually dominated by orientation gradients at subgrain boundaries, which may be removed by specifying a reasonably small threshold angle (e.g., \(2^{\circ }\)). KAM maps may be sensitive to measurement errors, making them noisy, which can be addressed by either including higher-order neighbors into the analysis or by denoising the EBSD maps from which they are derived. A more recent approach is the grain reference orientation deviation (GROD), which is the misorientation between each pixel orientation and the mean or a reference grain orientation, see Ref. [51] for a good overview. While GROD maps represent the orientation field with respect to a fixed point, KAM maps show the magnitude of the gradient, which may be interpreted as the first derivative of the orientation field. So these two approaches can be considered complementary. There exist some notable differences between the above approaches and dCOT: for GROD, a reference orientation must be available or wisely selected, which can only be done if the initial grain orientation is known from the start or still visible (e.g., because only part of the grain has deformed). As, for our dCOT approach, the initial grain orientation configuration is always known, cases such as homogeneous lattice rotation of entire grains or short-distance grain boundary migration can be identified, which is impossible using KAM or GROD. However, it is clear that in experiments, KAM and GROD may constitute the only available tools for visualizing and analyzing plastic grain deformation.

To extend the benefit from our dCOT analysis approach and furnish it with some quantification, we produced time- and depth-resolved maps representing the saturation value of the HSV (hue-saturation-value) color model, \(S_\mathrm {HSV}\), of the dCOT images. This quantity is defined as

$$\begin{aligned} S_\mathrm {HSV} = \frac{\mathrm {max}(R,G,B)-\mathrm {min}(R,G,B)}{\mathrm {max}(R,G,B)} \qquad \forall \quad \mathrm {max}(R,G,B)\ne 0 \end{aligned}$$
(2)

and \(S_\mathrm {HSV}=0\) otherwise. \(S_\mathrm {HSV}\) was laterally averaged over 20 evenly spaced tomographic sections, and the standard deviation of the value distribution was also recorded to retain some information about its lateral variance. The saturation value itself was then analyzed to study the time-dependent extent of subtle and stronger microstructural changes into the depth of the sample, while its standard deviation was used to quantify lateral equilibration processes within the evolving microstructure.

3 Results and Discussion

Out of the large amount of available data, we have selected 8 representative systems with parameter sets that showcase the influence of composition, normal pressure, and temperature on the microstructural development. These influences correspond to the variation of the stacking fault energy, the mechanical driving force, as well as a diffusive component acting on the system, respectively. The selected systems cover the regimes of incipient twinning and lattice rotation as well as grain refinement (as detailed in [30, 36]), with one system exhibiting almost no plastic deformation (pure Ni) and two systems bordering the regime that features strong plastic deformation with the formation of a shear layer (highest load/Cu content). In Fig. 3, we have arranged representative dCOT slices in such a way that time progresses from left to right, and the variation of an individual parameter is denoted by grouping within a colored frame. The system CuNi60 at 300 K and 0.5 GPa, framed in black, is chosen as a base system as it features incipient plastic deformation and can therefore be compared to all other systems. All systems of other alloy compositions, but otherwise identical loading parameters, are grouped together in a green frame. Those where only the normal pressure is varied are surrounded by a blue frame, and a variation of temperature is denoted by a pink frame.

At first sight, the dCOT images show the influence of changing composition (green frame), where both the depth and the severity of microstructural change decreases with increasing Ni content. The increase in Ni results in a higher stacking fault energy [52], which in combination with a higher shear modulus results in a greater resistance to emit partial dislocations from the GBs, possibly being the dominant deformation mechanism for the simulated grain size of 40 nm as the typical formation of dislocations from intragranular sources are suppressed [53, 54]. Therefore, the effect of a higher stacking fault energy is particularly relevant for the investigated nanocrystalline FCC compositions. The higher stress required to emit partial dislocations increases the resistance of the polycrystal to plastic deformation, ultimately leading to a negligible plastic deformation as observed for pure Ni at 300 K and 0.5 GPa [55]. The rise in stacking fault energy changes the deformation mechanisms from shear layer formation, via grain refinement and incipient twinning to purely elastic deformation with increasing Ni content [30, 36]. For the reference system CuNi60, the dCOT images clearly show incipient twinning, partial lattice rotation in the grain to the left, where the upper half of the grain is sheared and rotated, as well as slight grain boundary migration (stronger hues near the grain boundaries). Partial lattice rotation potentially leading to subgrain formation and eventually grain refinement is also frequently being observed in experimental studies [56]. In particular, grain boundary migration would not be apparent in the normal visualization with COT, demonstrating the strength of the dCOT visualization technique. When going to higher Cu content or higher loads and, therefore, stronger plastic deformation, the dCOT visualization comes to its limits since the images become very busy, and in particular the onset of shearing leads to a highly incongruous superposition of micrographs. This becomes obvious when looking at the composition with the lowest Ni content (CuNi5 at 0.5 GPa) or CuNi60 at the highest load of 1.1 GPa, where it seems that a new grain is developing in the middle until 2–3 ns. However, when looking at the conventional COT color scheme in Fig. 4, whose arrangement is identical to that in Fig. 3, it becomes obvious that the upper layer starts shearing to the right, therefore the grains tilt slightly in sliding direction. This shearing of an upper layer in direction of the sliding direction is typical for a frictional contact stemming from reorientation of dislocations and subsequent increase of misorientation, leading to elongated grains [57]. In the dCOT, this grain boundary migration and shearing make it seem as if a new grain is developing. Similar observations can be made for CuNi60 at the higher temperature, where a small grain seems to develop in the middle of the slice. Apart from this, twinning, partial lattice rotation, and dislocation activity can be very nicely seen in the upper part of the slice with a clear separation between deformed and undeformed zone [16]. In the dCOT it thus becomes clear that merely the part above the separating line rotates while the bottom part maintains its original orientation. Therefore, it is important to stress that dCOT is a tool to make subtle microstructural changes visible, but their analysis should always be performed with the aid of the original COT visualization.

Fig. 3
figure 3

Time-resolved comparison of dCOT images. The system CuNi60 at 300 K and 0.5 GPa (black frame) can be compared to all other systems in terms of composition (green frame), pressure (blue frame), or temperature (pink frame), with all other parameters kept constant (Color figure online)

Full size image

Additionally, as was mentioned earlier, the coloring scheme for dCOT cannot be unique by definition, so that it is sensible to review the times and points of interest in the COT color scheme known to every materials scientist because it is equivalent to EBSD-IPF, see Fig. 4. An explicit display of the initial configuration was not necessary, as the pure Ni system in row 4 exhibits such little microstructural change that it can double as a base for comparison.

Fig. 4
figure 4

Time-resolved comparison of COT images. The system CuNi60 at 300 K and 0.5 GPa (black frame) can be compared to all other systems in terms of composition (green frame), pressure (blue frame), or temperature (pink frame), with all other parameters kept constant (Color figure online)

Full size image

Figure 5 visualizes the trajectories of the orientations of the four marked grains shown in the COT panel over time under selected loading conditions. As before, the system CuNi60 at 300 K and 0.5 GPa is central in this comparison, as it is subjected to variations in composition (green axis “nickel content”, top), normal pressure (blue axis, bottom), and temperature (red axis, right). The orientation maps were arranged along the composition and load axes according to the microstructural effect on the system. Note that only those orientations that feature sufficient population densities were considered in this visualization to keep it as simple as possible.

The systems with varying composition show that lattice rotation is mostly constrained to the three grains marked with white, green, and blue pentagons in the figure being closest to the surface, which has similarly been shown in earlier studies [58]. Thereby, the two smaller grains marked with white and green pentagons are more prone to lattice rotation than the larger grain (blue pentagon), which can be traced back to the dependence of the rotation rate on the grain diameter [59]. The sub-surface grain marked with a red pentagon only shows minor lattice changes for CuNi25 and CuNi60, while it remains entirely unaltered in the pure Ni system. In general, pure Ni shows hardly any microstructural evolution, as already mentioned in the section discussing the micrographs. A decrease in Ni content (CuNi60) leads to a higher activity in the grains marked with white and green symbols, which rotate noticeably, particularly during the first 2–3 ns of the simulation. The surface grain (blue) and the sub-surface one (red) start deforming plastically for this composition, albeit weakly. A further reduction of the Ni content (CuNi25) results in a more pronounced lattice rotation of the two grains marked in white and green when compared to higher Ni contents, with the respective trajectories following more elaborate paths. The major feature of this composition in terms of texture is that the grain marked in blue that remains unchanged for higher Ni contents, shows in this case considerable twinning. However, as observed for higher Ni contents, the grain marked in red hardly deforms as it has no direct contact to the sliding interface. At the highest Cu content (CuNi5), we observe lattice rotations that populate large portions of the orientation triangle, including bifurcations that might be associated with grain refinement (potentially resulting from partial rotation of grains) as well as intermittent orientations that may arise from constant lattice rotation within a shear layer. In this case, also the sub-surface grain marked with red symbols rotates towards a higher-index orientation.

Fig. 5
figure 5

Grain orientation trajectories for various compositions, loads, and temperatures. The colors of the symbols/arrows in the orientation maps correspond to those of the pentagon symbols marking the respective grains in the COT panel on the left. Symbols occurring more than once at a given time indicate either a bifurcation due to splitting of a grain or a gradient due to partial lattice rotation (Color figure online)

Full size image

If we allow ourselves to fix the composition (CuNi60) and move along the pressure axis, we see similar trends as those observed for the composition variation. An increase of the contact pressure to 0.9 GPa leads to lattice rotations with trajectories similar to those observed for CuNi25 at 0.5 GPa, with apparent twinning and clear lattice rotations in the grains marked with green and white. Increasing the pressure further to 1.1 GPa leads to a complex set of trajectories reminiscent of those for CuNi5 at 0.5 GPa, with the surface grain marked with white symbols rotating first in the direction of [1 0 0] and then on to [1 1 1] as a shear layer starts evolving.

Finally, when increasing the temperature to 600 K, the map again looks reminiscent of that representing CuNi25 at 300 K, featuring twinning and the following of elaborate orientation trajectories mentioned earlier, with some additional twinning in the grain marked with white symbols. The latter seems to be more particular to this diffusive loading of the system, leading to more dislocation activity inside the grain. Despite similar trajectories shown here independent of load, composition, or temperature is varied, where higher-index orientations evolve towards lower-index ones, a preferential grain orientation due to mechanical loading cannot be observed. This may probably be attributed to the limitations of time inherent to MD simulations as well as the small number of grains (around 15), making it difficult to obtain any reasonable statistics that would imply an evolving texture.

In recent experimental work, Cai et al. performed pin-on-disc experiments with a carefully selected copper nickel bronze [60]. They designed their experiment in such a way that the plastically deformed zone is much smaller than the initial grain size, and they optimized their alloy and the loading conditions such that lattice rotation dominates over grain subdivision and subgrain formation. They found that the majority of the grains rotated around an axis parallel to the transverse direction of the pin. In some additional analyses on our simulation data, we found a definite trend corresponding to these experimental results, i.e., the majority of our grains also rotate about axes that are nearly perpendicular to both the loading and the sliding directions. In Fig. 5, this applies to the entire grain marked by the green pentagon, as well as the upper portions (upon formation of a grain boundary) of the other two surface grains (marked white and blue).

In Fig. 6, we compare the time- and depth-resolved, laterally averaged HSV saturation value of the dCOT images for our representative selection of systems. This quantity gives a good overview over when, to which depth into the sample, and how strongly microstructural changes (including lattice rotation and twinning) occur. Note that the pure Ni system was omitted in this comparison, as it did not feature much contrast regardless of time and depth. The lower and upper bounds of the logarithmic color axis were chosen ex post to yield the highest contrast without the need to re-scale between systems. Although the HSV color scheme is periodic (it begins and ends at red), it is clear which occurrences of red denote an HSV saturation of 0.08 or below (lower left) and which ones denote 0.6 or above (top right). As before, the base system for this comparison is CuNi60 at 300 K and 0.5 GPa, which exhibits the weakest microstructural response in this representation. As either the copper content, the normal pressure, or the temperature are increased, we observe that the onset of plasticity occurs earlier. At the highest pressures or Cu contents, plasticity can be observed only about 200 ps after sliding starts. Furthermore, the rate at which the plastic deformation propagates into the depth of the sample increases, as does the equilibrium depth of the plastically deformed zone, which may or may not have been reached after 5 ns depending on the composition and boundary conditions. As with the trajectories in orientation space discussed in Fig. 5, pairwise similarities exist between some systems. In the order of increasing deformation intensity, CuNi60 at 300 K and 0.7 GPa behaves like CuNi60 at 600 K and 0.5 GPa; CuNi60 at 300 K and 0.9 GPa corresponds to CuNi25 at 300 K and 0.5 GPa; and CuNi60 at 300 K and 1.1 GPa shares similarities with CuNi5 at 300 K and 0.5 GPa. That said, there seems to exist a consistent difference between the onset of plasticity induced by higher pressure or temperature compared to an increased Cu content. In the former two cases, we observe very early slight increases in the HSV saturation value that permeate large portions of the sample, whereas the latter case keeps deeper lying regions undeformed for longer periods of time, but plastic deformation then sets in much faster and to a deeper extent.

Fig. 6
figure 6

Time- and depth-resolved HSV saturation value of laterally averaged dCOT images for a set of representative systems. Note the logarithmic color scale (Color figure online)

Full size image

After studying the time- and depth development of the mean saturation value, it is also highly instructive to compare the corresponding maps of the lateral standard deviation of the HSV saturation, see Fig. 7. This quantity reflects how homogeneously the plastic deformation is distributed laterally at a given time and depth. As before, pure Ni was excluded from this figure, as it does not exhibit enough plastic deformation to produce any contrast. Again, the color range was adjusted for maximum contrast and is the same for all the maps. While the saturation maps in Fig. 6 generally showed a monotonic increase of the depth to which the system is affected by plastic deformation, possibly saturating to some constant value, these maps feature interesting fluctuations that correspond to dynamic equilibration processes. All maps except those for CuNi60 at 0.5 GPa reflect such behavior, which is basically a lateral homogenization of the plastically deformed zone. This homogenization process is particularly illustrative for CuNi60 at 0.7 GPa. In this case, the HSV saturation value (Fig. 6) exhibits a monotonic growth of the deformation depth with sliding distance/time. However, this microstructural refinement does not proceed homogeneously, as highlighted by the lateral standard deviation of the HSV saturation. The latter features a decrease in value after 3.5 ns, suggesting a homogenization of the deformed region.

Fig. 7
figure 7

Time- and depth-resolved lateral standard deviation of the HSV saturation value of dCOT images for a set of representative systems

Full size image

To discuss this effect with an example, the HSV saturation for pure Cu at 0.4 GPa is shown in Fig. 8, along with its standard deviation and corresponding representative dCOT images. The tomographs are particularly suitable for visualizing the homogenization of the microstructure as after 5 ns sliding, all cross-sections show a well-defined deformed region comprising the 10 nm right below the surface. However, previous time steps reveal that the development of this discontinuity is inhomogeneous.

Fig. 8
figure 8

Lateral dynamic equilibration of plastic deformation in pure Cu sliding at 300 K and 0.4 GPa. At the top, the time- and depth-resolved values of the dCOT HSV saturation value and its standard deviation are mapped. For illustration of the lateral inhomogeneity of the deformation process, five dCOT slices, equally spaced along the y-direction perpendicular to normal pressure and sliding, are shown at 1 ns intervals

Full size image

According to the tomographic images, the inhomogeneity is caused by two main factors. One contribution is given by the diverse deformation stages that can be observed in the slices of the polycrystalline aggregate as a consequence of the varying deformation imparted by the rough counterbody. Due to its surface topography, the imparted deformation has a stochastic component that tends to homogenize for sufficiently large sliding distances. However, an inhomogeneous evolution of the deformed region can be also observed within a single tomograph, even though the amount of strain is identical in this case (or similar, if the asperity of the counterbody is not sliding fully parallel to the slice). Let us consider the third slice (at \(y=40\) nm) as an example. Initially, after 1 ns, only the upper part of one grain is deformed, while the neighboring regions start following after 2 ns. A similar observation can be made in the first slice (\(y=6\) nm), where until 4 ns, the entire plastic deformation is confined to a single grain. The reason for the inhomogeneous deformation within one slice seems to be related to the different lattice orientations of the grains, as certain orientations will be oriented more favorably for undergoing slip by having slip systems with a higher Schmid factor. We calculated the Schmid factor for the initial grain configuration (which is essentially the same for all systems) using the MTEX toolbox assuming uniaxial stress in sliding (x) direction. In case of the slice in the top row (\(y=6\) nm), the right grain has a Schmid factor of over 0.45 while the Schmid factor of the left grain is below 0.35. As can be seen, deformation starts in the grain with the higher Schmid factor and is much more pronounced after 5 ns. By contrast, for the third slice (\(y=40\) nm), the Schmid factors of the surface grains are more similar, with all values exceeding 0.45, leading to a more homogeneously distributed deformation. However, also in this case, the right grain that carries most of the initial deformation has a slightly higher Schmid factor.

As deformation progresses, the previously deformed grains undergo strain hardening, raising their critical resolved shear stress and making them progressively more resistant against plastic deformation. This allows the neighboring undeformed grains to catch up with the deformation, which ultimately leads to a well-defined discontinuity with increasing sliding time. Even though this homogenization of the deformed sub-surface region was observed in a number of testing conditions, it seems to be more evident at intermediate loads with respect to the varying resistance to plastic deformation for the respective compositions (region of grain boundary dominated processes). For instance, the most exemplary cases were Cu and CuNi5 at 0.4 GPa, CuNi25 at 0.5 GPa, CuNi60 at 0.7 GPa, and Ni at 1.5 GPa. At lower normal pressures, where plasticity was confined to the immediate vicinity of the interface, or at high loads, where plastic deformation immediately spread throughout the polycrystalline aggregate, such a homogenization process could not be observed. This irregular formation of the deformed layer seems to contradict the experimental observations made in Ref. [17]. However, the authors found that the dislocation trace line was wavy and rather difficult to identify in some cases. These cases corresponded to those obtained with the highest contact pressure (1953 MPa), but achieved using a combination of low load and small indenter diameter. For higher loads, maintaining the contact pressure with a correspondingly larger indenter diameter or smaller contact pressures, the dislocation trace line was again present and defined below the surface. This could be an indication that particular loading conditions lead to an initial inhomogeneous formation of the discontinuity between the deformed and the undeformed sub-surface region, even though further research would be required to verify this statement.

4 Conclusion and Outlook

The early time evolution of near-surface microstructure in metal sliding contacts play a significant role in the long-time tribological behavior of bearings, piston rings, and gears. For example, near-surface microstructures that harden to the point of brittleness can spall off and result in enhanced wear. Here, we addressed the question of early time microstructure evolution in sliding contacts through a series of large-scale MD simulations in nanoscale polycrystalline FCC metals. Atoms in the sliding counterbody are held fixed, mimicking a much harder slider. In order to understand the influence of process conditions (temperature and normal pressure) and alloy composition, we fixed the initial microstructure for all simulations and correlated the evolving microstructure to this initial microstructure using a novel dCOT visual analysis. While EBSD has become a standard tool to quantify a microstructure, dCOT compares microstructures at two different times to expose changes in microstructure, such as grain boundary motion, lattice rotation, twin formation, and dislocation plasticity.

What can be learned from this? To illustrate the power of the analysis, we performed simulations of copper/nickel alloys ranging from pure copper to pure nickel. Copper readily displays plasticity and grain reorientation, while adding nickel increases the stacking fault energy, and we observe increasing resistance to plastic deformation. Pure nickel is dominated by elastic deformation. Increasing both temperature and pressure enhances plastic deformation. The extent of the plastically deformed layer is revealed by analyzing the time- and depth-resolved HSV color saturation value of the laterally averaged dCOT images. Our analyses serve to provide insight into the initial development stages of the discontinuity that arises between deformed and undeformed regions.

The two-time correlation revealed by the dCOT offers an interesting prospect. While both temperature and pressure are process parameters, both alloy composition and initial microstructure are material design parameters. The dCOT analysis, when integrated into a computational materials engineering framework, can enable optimal material design for tribological and other deformation problems.