Thermostat Influence on the Structural Development and Material Removal during Abrasion of Nanocrystalline Ferrite

TL;DR

This study uses molecular dynamics simulations to show that the method of heat removal significantly affects the material removal and structural development during nanomachining of ferrite, highlighting the importance of realistic thermal modeling.

Contribution

It introduces an electron-phonon coupling based thermostat in MD simulations, providing more accurate temperature gradients and revealing the impact of heat removal mode on nanomachining outcomes.

Findings

Thermostat mode influences grain and chip morphology.

Average shear stress remains temperature independent up to 1 GPa.

Heat removal method affects material response beyond 0.4 GPa.

Abstract

We consider a nanomachining process of hard, abrasive particles grinding on the rough surface of a polycrystalline ferritic work piece. Using extensive large-scale molecular dynamics (MD) simulations, we show that the mode of thermostatting, i.e., the way that the heat generated through deformation and friction is removed from the system, has crucial impact on tribological and materials related phenomena. By adopting an electron-phonon coupling approach to parametrize the thermostat of the system, thus including the electronic contribution to the thermal conductivity of iron, we can reproduce the experimentally measured values that yield realistic temperature gradients in the work piece. We compare these results to those obtained by assuming the two extreme cases of only phononic heat conduction and instantaneous removal of the heat generated in the machining interface. Our discussion of the differences between these three cases reveals that although the average shear stress is virtually temperature independent up to a normal pressure of approximately 1 GPa, the grain and chip morphology as well as most relevant quantities depend heavily on the mode of thermostatting beyond a normal pressure of 0.4 GPa. These pronounced differences can be explained by the thermally activated processes that guide the reaction of the Fe lattice to the external mechanical and thermal loads caused by nanomachining.

Tables

Table 1. Table 1: Overview of the coupling times for the Langevin thermostats acting on the substrate base ( λ base subscript 𝜆 base \lambda_{\mathrm{base}} ) and on the rest of the Fe atoms ( λ bulk subscript 𝜆 bulk \lambda_{\mathrm{bulk}} ).

thermostatting scheme	${\lambda }_{\mathrm{base}}$ [ps]	${\lambda }_{\mathrm{bulk}}$ [ps]
base	0.5	—
electron-phonon coupled	0.5	3.5
full	0.5	0.5

Equations

$$\lambda =\frac{2m_e\kappa E_F}{{\Theta }_D{T}_{0}Ln{e}^2{k}_{B}Z},$$

$${\Phi }_q=-\kappa (T)\frac{dT(z)}{dz}.$$

Full text

\phone

+43 (0)2622 81600 161 \fax+43 (0)2622 81600 99 AC2T research GmbH] AC2T research GmbH, Viktor-Kaplan-Straße 2, 2700 Wiener Neustadt, Austria

AC2T research GmbH] AC2T research GmbH, Viktor-Kaplan-Straße 2, 2700 Wiener Neustadt, Austria

AC2T research GmbH] AC2T research GmbH, Viktor-Kaplan-Straße 2, 2700 Wiener Neustadt, Austria

Hamburg University of Technology] Institute of Advanced Ceramics, Hamburg University of Technology, Denickestraße 15/K, 21073 Hamburg, Germany

Vienna University of Technology] Institute of Applied Physics, Vienna University of Technology, Wiedner Hauptstraße 8–10/134, 1040 Vienna, Austria

Note: This document is the unedited Author’s version of a Submitted Work that was subsequently accepted for publication in ACS Applied Materials & Interfaces, copyright ©American Chemical Society after peer review. To access the final edited and published work see http://pubs.acs.org/doi/full/10.1021/acsami.7b01237.

Please cite this work as:

S. J. Eder, U. Cihak-Bayr, D. Bianchi, G. Feldbauer, and G. Betz, “Thermostat Influence on the Structural Development and Material Removal during Abrasion of Nanocrystalline Ferrite,” ACS Applied Materials & Interfaces, DOI: 10.1021/acsami.7b01237, 2017.

1 Introduction

The response of the near-surface microstructure of a work piece or a system component subjected to abrasive conditions can strongly influence its properties and behavior. This aspect of surface modification is of high interest to the surface finishing industry, where material is deliberately removed from the surface to obtain a desired roughness or texture, but also for tribological applications, where the – usually unintended – wearing of surfaces affects a system’s performance and service life. Depending on the relative velocities of the surfaces in mechanical contact, the friction and possible plastic deformation occurring in the interface can lead to the generation of considerable amounts of heat, which has a large impact on the development of the microstructure close to that interface.

The geometrical tolerances for components in various applications have increasingly entered into the nanometric length scale, e.g., in high-gloss finishing for optical purposes 1 or in nano/micro-electromechanical systems (NEMS/MEMS) 2. As the mechanisms governing material removal and wear, but also those responsible for changes to the crystal structure have their foundations in the atomistic nature of the component or work piece, it makes sense to consider the respective processes from a nanoscopic point of view.

The shift in tribology towards the nano length scale took place over the last decades and showed that commonly used macroscopic tribological models as well as laws cannot be readily applied at this scale 3, 4, 5. Particularly, nanoscale wear processes 6, 7, 8 are not fully understood because of their highly complex nature 9. Over the last years a high effort has been put into a better understanding of single-asperity contacts, which can be probed experimentally in detail using an atomic force microscope (AFM) 10. This approach is of high importance from a fundamental point of view; however, many realistic, nano-technological applications such as NEMS/MEMS are dominated by multi-asperity contacts 11, 12. In such systems nanowear poses a major limitation for the lifetime, performance, reliability, or the overall usability. Thus, it is imperative to investigate such multi-asperity contacts exhaustively. A further pressing issue regarding contacts at the nanoscale is the calculation of the real contact area. Commonly applied continuum mechanics approaches are hardly applicable at this length scale 13; therefore, various methods have been proposed to calculate the contact area based on atomistic principles, see for example refs 12, 14, 15.

As addition to experimental techniques and to obtain a more complete picture, classical molecular dynamics (MD) simulations have proven very useful to simulate nanoscale tribological systems 10. Modern computing power allows tracking the time development of fully atomistic systems consisting of several million atoms without the a priori assumption of constitutive equations governing the systems’ behavior, see refs 16, 17 and the references therein for some examples. One important aspect of non-equlilbrium MD (NEMD) simulations, i.e., where energy is continuously introduced into the system via external constraints such as imposed forces or velocities, is the thermodynamically correct removal of the generated heat. Several concepts have been proposed to do this in a manner consistent with some thermodynamic ensemble, such as the Nosé-Hoover 18, 19, the Berendsen 20, or the Langevin thermostats 21, 22, to name only a few.

The feasibility and quality of a NEMD simulation is strongly affected by the choice of the thermostat as well as its subsequent configuration. Particularly, this includes the decisions on which parts of the system the chosen thermostat acts upon and on the coupling strength between system and thermostat. In the following, three distinct scenarios will be presented. For a nanomachining simulation as considered in this work it would, in principle, make sense to apply the thermostat only to a part of the system some distance away from the tribological interface. Such a setup would minimize the interference with the physics of the processes occurring close to the surface. This implies that there is a large heat sink attached to the base of the explicitly considered near-surface region, namely the bulk. Unfortunately, typical interaction potentials modeling the behavior of metals, such as EAM 23 or Finnis-Sinclair 24, do not explicitly include electrons and are therefore not able to sufficiently reproduce the thermal conductivity of a metal. As the modeled systems grow ever larger, resulting in larger distances between contact zone and heat sink, the grossly overestimated temperature gradients additionally lead to highly unrealistic absolute temperature differences. We will henceforth refer to this scenario as “base-thermostatted”.

An alternative approach is to thermostat the entire work piece. Most frequently, the thermostat is strongly coupled to the atoms so that heat is almost instantaneously transported away. While this certainly allows one to maintain good temperature control over the work piece, one may completely neglect the thermal influence of the modeled mechanical process itself by eliminating heat conduction altogether. We will call this scenario “fully thermostatted”.

When it comes to atomistic simulations of systems featuring considerable heat transfer through metals, many authors do not mention which thermostatting method they use to control the system’s temperature. If they do, they seldom disclose how they parametrized the thermostat or justify why they chose the particular parametrization. Without any claim to completeness, we list some examples of otherwise high-quality simulations that are very likely underthermostatted 25, 26, 27, 28 due to base-thermostatting or possibly overthermostatted 29, 30, 31, 32 due to full thermostatting. Out of these examples, it seems that only Shiari et al. 30 are aware of some implications of their simplified thermostatting efforts and concede that another group reported that material removal is greatly influenced by the thermal conditions in the shear zone 33.

Gill 34 gives a good overview of NEMD as well as multiscale modeling of heat conduction in solids. However, most of this extensive review focuses on non-metallic materials. The presented methods intended for metals are either quasi-static, such as the two-temperature method mainly used for modeling the laser annealing of voids 35, 36, while others require a multiscale approach to the problem via dynamic coarse graining 37 or featuring coupling schemes to continuum 38. The latter apply their ad hoc technique to an examination of frictional heating during sliding by solving the heat equation and imposing a thermal conductivity.

In this work, we adopt the concept of electron-phonon coupling as laid out in refs 39, 40. This does not require us to implement a time-consuming multiscale approach, but rather assumes that due to their high mobility within the metal, the electrons can be used as an implicit heat sink permeating the substrate. So far, this setup is identical to what we described above as “fully thermostatted”. The crucial point of putting the electron-phonon coupling concept into practice is correctly parametrizing the thermostat so that it reflects how the electrons in the metal interact with the lattice vibrations. We will consider the two extreme cases of thermostatting only the substrate base (“base thermostatting”) and thermostatting the entire substrate with a strongly coupled thermostat (“full thermostatting”) and compare these two limiting scenarios to results obtained with a carefully parametrized thermostat that attempts to reflect the electron-phonon coupling in the metallic work piece as closely as possible. Additionally, we check if considerably extending the heat treatment period of the initial substrate configuration significantly changes the results of the abrasion simulations. We discuss how the application of the electron-phonon coupled thermostatting scheme affects some typically evaluated tribological quantities such as friction, material removal/wear, contact area, and surface topography, as well as the microstructural development of the work piece. This constitutes an important step towards making the results of atomistic non-equilibrium simulations of materials undergoing plastic deformation and subjected to high temperature gradients more reliable.

2 Modeling and Simulation Details

All our simulations were carried out using the open-source MD code LAMMPS 41. The exact procedure of our model construction is described in an earlier work 42. It consists of a rough 60$\times$60$\times$20 nm3 $\alpha$-Fe polycrystal including randomly oriented abrasives with Gaussian size-distribution. A representative snapshot of the 3D model is shown in the top right corner of Fig. 1, along with some diagonally arranged tomographs colored according to the grain orientation. The mean equivalent diameter of the bcc Fe grains in the originally constructed 3D-periodic 60$\times$60$\times$60 nm3 cube is 12.7 nm, corresponding to an average grain volume of 1080 nm3. Due to the introduction of the surface at $z=20$ nm, out of the 128 remaining grains, the cleft ones no longer have the shape and size of their original Voronoi cells. Therefore, the final overall size distribution of the grains features 40 grains with volumes of less than 200 nm3, a size that was not present in the original size distribution at all, while 20 grains have sizes between 200 nm3 and 400 nm3, compared with only 7 in the original polycrystal. There are thus approximately 75 grains in the initial substrate that almost fully retain their original geometry. The grains are initially randomly oriented, and the surface topography has a fractal dimension of 2.111, an RMS roughness of 0.7 nm, and a lower frequency cut-off producing a typical lateral roughness feature extent of 23 nm. The atoms located in the lower 3 Å  of the simulation box are kept rigid to emulate bulk support, and the Fe–Fe interactions within the substrate are governed by a Finnis-Sinclair potential 43.

Two slightly different initial configurations were produced by subjecting them to two different types of heat treatments, abbreviated “HT1” and “HT2” henceforth, with the temperature curves shown in Fig. 2. In both heat treatments, a Langevin thermostat 21, 22 controlled the temperature of the entire substrate (abbreviated “full” in some figures) during the ramps, while during the constant temperature intervals only a layer with a height of 0.3 nm at the substrate base was thermostatted (“base”). In this way, the annealing process is efficient while allowing the crystal grains to reorganize without too much thermostat interference. This basic Langevin thermostat was parameterized with a time constant of 0.5 ps, corresponding to strong electron-phonon coupling in metals 44. For a comparison between the substrate microstructures after the two heat treatments, see the selected substrate tomographs in the center of Fig. 3. These tomographs are colored according to grain orientation as in electron backscatter diffraction, using the inverse pole figure coloring standard. The orientations were calculated using polyhedral template matching 45 as implemented in OVITO 46, and the color rendition was carried out using the MTEX toolbox 47, 48 for Matlab. The differences between HT1 and HT2 are subtle, even when magnified, and best visible in or around small near-surface grains, see the comparison for slice #1 shown in the center of Fig. 3. The comparison for slice #2 (right side) gives an idea of the largest occurring differences between the two heat treatments, where the small purple and orange grains located in the upper central part of the section have disappeared in favor of neighboring grains after HT2.

The 18 abrasives are modeled as rigid truncated octahedra with diameters ranging from 6–14 nm, randomly rotated and distributed laterally, leading to an effective surface coverage of 37%. As in ref 49, the interactions between abrasive particles and the Fe surface are modeled using a Lennard-Jones potential with the parameters $\varepsilon_{\mathrm{LJ}}=0.125$ eV and $\sigma_{\mathrm{LJ}}=0.2203$ nm. These values are similar to refs 25, 50 and yield an interaction roughly one order of magnitude weaker than the one for Fe–Fe, which means that there will be some adhesion between the abrasives and substrate as well as the wear particles. The abrasives are pulled across the surface at a sliding velocity in $x$ direction of $v_{\mathrm{(slide)}}=80$ m/s and an angle of $6.42$° with the $x$-axis, so that they re-enter the simulation box at different $y$ positions every time they pass the periodic box boundaries and therefore never follow exactly in their own grinding marks. The relative abrasive positions are locked throughout the simulations, and abrasive rotation is disabled. Furthermore, the abrasives can change their $z$ position collectively depending on the topography, similar to a grindstone, but not individually. Because of their rigidity, the abrasive particles themselves are not subject to any wear. The normal pressure $\sigma_{z}$ on the abrasives (defined as the total force in $-z$ direction acting on the abrasive particles divided by the lateral cross-section of the simulation box, 3595.4 nm2) is kept constant at values ranging from 0.1 to 0.9 GPa for a simulation time of 5 ns. During the grinding simulations, the Langevin thermostat acts only in $y$ direction, (nearly) perpendicular to the directions of normal pressure and grinding, so as not to overly interfere with these external constraints.

3 Thermostat Parametrization

The electronic contribution to the thermal conductivity $\kappa(T)$, which dominates in metallic systems, is not accounted for in the interaction potentials applied in classical MD. We calculated $\kappa(T)$ obtained from the state-of-the-art Finnis-Sinclair potential for Fe used in this work 43, which is based only on the phononic contribution, for temperatures ranging from 300 to 900 K using the Green-Kubo formulas, which relate the ensemble average of the auto-correlation of the heat flux $q$ to $\kappa(T)$. When comparing the results to the experimental TPRC data series 51, we found that $\kappa(T)$ calculated via MD correlates well with the experimental values, but that the MD values are, on average, 4.5 times lower than the experimental ones.

We will now describe a “smart” thermostatting scheme, based on refs 39, 40, which attempts to reflect the coupling between electrons and phonons in metals. Such an approach is justified as long as the externally imposed velocities lie well below the speed of sound in the work piece, and the mean free path for the electron-phonon interactions is short enough to be accommodated in the work piece along the direction of heat conduction. With our grinding velocity amounting to less than 2% of the speed of sound in iron, the first condition is easily met. An estimation based on the assumptions in ref 52 yields an electron-phonon mean free path of approximately 2.5 nm, which is roughly one order of magnitude smaller than our work piece thickness. The Langevin thermostat acting on the bulk of the substrate is parametrized with a coupling constant corresponding to the characteristic time that it takes the lattice vibration energy of the ionic system to be transferred to the electron gas, which, due to its high thermal conductivity, can be considered an implicit heat bath. Based on the Sommerfeld theory of metals, an estimation for this electron-phonon coupling time $\lambda$ can be expressed as 40

[TABLE]

where $T_{0}$ is the temperature of the heat bath, $\Theta_{D}$ the Debye temperature, $L$ the Lorenz number, $n$ the free-electron density, $Z$ the valence, $\kappa$ the thermal conductivity, $E_{F}$ the Fermi energy, $k_{B}$ the Boltzmann constant, and $e$ and $m_{e}$ the electron charge and mass. If we feed this expression with the data for Fe at the expected mean substrate temperature of $T\simeq 310$ K, we obtain a $\lambda$ of about 0.5 ps.

On the other hand, we can roughly estimate the thermal conductivity of our substrate by assuming that almost all the frictional energy dissipated in the interface must find its way to the heat bath eventually. We neglect the heat that flows from the sliding interface into the wear particles as well as radiation losses from the surface. The former contribution depends on the load and amounts to several percent of the total, while the latter are not even included in our MD model, but are roughly six orders of magnitude smaller than the friction power at the surface temperatures we encounter in our simulations. We thus set the one-dimensional heat flux to $\Phi_{q}=\mu\sigma_{z}v_{\mathrm{(slide)}}$, where $\mu$ is the coefficient of friction, $\sigma_{z}$ is the normal pressure, and $v_{\mathrm{(slide)}}$ is the sliding speed. While $\sigma_{z}$ and $v_{\mathrm{(slide)}}$ are kept constant as external constraints, $\mu$ turns out to maintain a value of 1 throughout the simulations discussed in this work (see Sec. 4.1 for details), so we may also treat it as a constant. The effective thermal conductivity measured from our grinding simulations can then be obtained via the Fourier law

[TABLE]

We can now search for the optimum value of $\lambda$ starting from the value 0.5 ps estimated using Eq. (1) to parametrize our Langevin thermostat for practical use. This is done by performing an abrasion simulation at constant $\Phi_{q}$ up to the point where it reaches a sufficiently linear temperature gradient $dT(z)/dz$ within the substrate close to the machining interface and checking how well that gradient corresponds to the experimental thermal conductivity of Fe 51 at the mean substrate temperature, see Fig. 4 for some associated temperature profiles at various loads and values of $\lambda$. We calculated $\kappa(T)$ for several values of $\lambda$ in this way and found that the best reproduction of the experimental thermal conductivity of Fe can be achieved with $\lambda=3.5$ ps throughout our desired range of normal pressures, c.f. the inset in Fig. 4. Simulations of ion beam mixing and damage production in Fe 53 with an implementation of the electron-phonon coupling model of ref 40 also came to the conclusion that the time constant best reflecting the experimental data was several times larger than the one calculated via Eq. (1), arriving at values between 1.6 and 2.2 ps. Other work dealing with metal ablation by picosecond laser pulses 36 arrives at an electron-phonon coupling time of 7 ps. We can thus conclude that the estimate for $\lambda$ calculated via Eq. (1) can only serve as a basis for empirically finding the electron-phonon coupling time that best reflects experimentally observable quantities. Table 1 summarizes the coupling times parametrizing the Langevin thermostats acting on the base and the bulk of the substrate to implement the three discussed thermostatting schemes.

An additional benefit of the electron-phonon coupled thermostatting scheme discussed above is that the cooling of the removed matter (e.g., via a coolant fluid) is, at least for all practical purposes, accounted for in this way, i.e., even if the wear debris are already detached from the substrate surface and would therefore no longer be thermally coupled to a base thermostat, they can cool off by themselves.

It should be noted that the theoretically calculated maximum surface temperature occurring due to dry sliding 54, 55, 56 can rise considerably higher than the surface temperatures we observe in our simulations. This is because firstly our simulations last only 5 ns, while macroscopic equilibrium temperatures are reached after much longer periods of sliding. Secondly, the Langevin thermostat representing the heat sink of the work piece is placed only 20 nm from the friction interface while being kept at a constant and rather low temperature of 300 K, so that by fixing the allowed temperature gradient in the near-surface region, we also limit the maximum surface temperature.

4 Results and Discussion

4.1 Time, load, and temperature dependence of tribological quantities

In Fig. 5, we present the time development of the wear depth, the arithmetic mean height, and the RMS roughness of the surface topography for electron-phonon coupled thermostatting. The wear depth is defined as the combined volume of all Fe atoms traveling faster than 90% of the abrasives’ speed, so that they can be considered attached to the grinding particles 57, divided by the cross-section of the simulation box $A_{\mathrm{nom}}=3595$ nm2. Compared to our previous work on monocrystalline ferrite substrates and considerably more regular lateral distributions of abrasive particles 49, 58, 59, 60, the transition between the running-in and the steady-state regime is much less pronounced, see the kink in most curves in Fig. 5 (a) at $t\approx 0.3$ ns. As, by our definition, the substrate consists of all atoms moving in $+x$ direction at less than 10% of the grinding speed, the substrate topography does not include the shear zone 57. However, during the initial impact of the abrasives on the substrate and the brief period following it, which is characterized by strong oscillations of the contact forces due to the substrate elasticity, a shear zone is formed that can momentarily be up to two orders of magnitude larger than its equilibrium size and may extend up to 10 nm into the substrate. This explains the large deviations from the general trend during the first 0.2 ns in panels (b) and (c), where the topographic quantities are shown.

The time evolution of the wear depth exhibits a steady increase for the high pressure simulations, which is reflected in a steady decrease of the mean substrate height $z_{\mathrm{subst}}$. The low pressure simulations hardly show any reduction of $z_{\mathrm{subst}}$ and a very small increase of the wear depth $h_{w}$, restricted to the first 0.5 ns. This suggests that up to 0.3 GPa the system does not wear significantly, but that the reduction of the surface roughness, as depicted in Fig. 5 (c), is merely a result of material transfer from the peaks into the valleys, i.e., abrasion of the asperity tips and subsequent re-crystallization onto the substrate. The wear particles are formed at the very beginning for low normal pressures, and although individual particles may change their sizes over time, the combined volume of all the wear particles remains fairly constant for the rest of the grinding simulation. For the high pressure cases of 0.9 GPa and 0.75 GPa, the time evolution of $S_{q}$ is unstable, and the final roughness approximately equals the initial value.

In Fig. 6, we show exemplary $xy$ projections of the real atomic contact area, colored according to contact temperature after 5 ns of grinding. As explained in detail in refs 57, 61, this contact area is obtained by establishing which substrate atoms are in contact with the abrasives or the associated wear particles via a distance criterion, and then multiplying the number of these contact atoms with a constant per-atom contact area. The contact temperature for each contact atom is calculated via the kinetic energy averaged over a spherical control volume with a radius of 1 nm after subtracting the drift velocity of this control volume, cf. ref 42. The three examples shown in Fig. 6, from left to right, are for normal pressures of 0.3, 0.6, and 0.9 GPa, respectively. The top row shows results for base-thermostatted substrates, the center one for electron-phonon coupled thermostatting, and the bottom row fully thermostatted substrates. The contact area increases with load independently of the chosen thermostatting procedure. In detail, the shape and size of individual patches of the contact area are different depending on the chosen thermostatting. However, its more general characteristics remain unchanged since these are primarily determined by the asperities of the substrate as well as the size and distribution of the abrasive particles. Clearly, the contact temperature dependence on the normal pressure is considerable when applying only a base thermostat. For the electron-phonon coupling and the fully thermostatted substrate, increasing the normal pressure does not influence the contact temperature much. The main difference between these two variants is the temperature of the hottest spots within the contact area, which is about 100 K higher if electron-phonon coupling is applied. Additionally, the local temperature distribution within the contact area patches is more inhomogeneous for electron-phonon coupling.

It seems that the choice of the heat treatment procedure has little influence on the normal pressure dependence of some tribological, topographical, and thermodynamic quantities shown in Fig. 7. The only noticeable difference occurs for the final RMS roughness at the highest load, see Fig. 7 (e), but here the error bars are also very large and overlap to such extent that it must be questioned if the difference between HT1 and HT2 is statistically relevant.

It is more striking that the pressure dependence of the shear stress $\sigma_{x}$ does not appear to be influenced by the thermostatting procedure at all. This is in agreement with other literature 62, where it is reported that the substrate temperature, varied from 0 K to 500 K, has little effect on what is there called the “scratching force”. We therefore carried out additional calculations at a considerably higher normal pressure of $\sigma_{z}=2.5$ GPa with base- and electron-phonon coupled thermostatting. They revealed that under these severe conditions a difference between the obtained shear stresses emerges. While the electron-phonon coupled results remain on the straight line suggested in Fig. 7 (a), the results for the base thermostat lie approximately 25% lower. This might be explained by a pronounced softening effect of the polycrystal due to the high temperatures in the surface zone, where the base-thermostatted substrate ground at $\sigma_{z}=2.5$ GPa reaches $z$-averaged temperatures of 1160 K, while the $T$-maxima in the contact spots even rise up to 1350 K. Compared to the much lower surface temperatures of only 450 K for the electron-phonon coupled substrate, the hotter and thus softer base-thermostatted substrate opposes much less resistance to shear, equivalent to reduced work hardening of macroscopic samples at high temperatures.

As expected, the wear depth in Fig. 7 (b) and the arithmetic mean height in Fig. 7 (d) exhibit opposing trends, both showing almost identical behavior for all thermostatting procedures up to 0.4 GPa and a splitting up into three distinct branches for greater normal pressures. The fully thermostatted substrates are worn less than the electron-phonon coupled ones, which in turn wear less than the base-thermostatted ones. The final RMS roughness of the surface in Fig. 7 (e) already shows some differences at a normal pressure of 0.4 GPa, but this difference remains small up to 0.6 GPa. At the highest normal pressure, however, the “best” achievable roughness of the base-thermostatted substrates is almost double that of the ones obtained with the other two thermostatting procedures. Note that the initial RMS roughness is 0.7 nm, which means that grinding at normal pressures of 0.75 GPa or more with a base-thermostatted substrate does not produce a smoother surface. Whereas the pressure dependence of the final roughness is almost linear for full and electron-phonon coupled thermostatting, which is in good agreement with ref 49 (where the substrates were so thin that the thermostatting procedure would not have had much influence on the result), the base-thermostatted systems show a strongly superlinear pressure-dependence. This might be attributed to the following: in the base-thermostatted systems, heat transport from the chips towards the heat sink is only possible via phononic heat conduction through the contact regions with the substrate. The chips can therefore reach temperatures close to 1000 K, making them soft. This, in turn, leads to their occasional slumping over to the sides of the abrasives, leaving high scratch ridges in the likewise softened near-surface substrate. A linear increase of the RMS roughness with the normal pressure implies that the substrate compliance does not change within the substrate volume that interacts with the abrasives, even at higher loads. Only for the base-thermostatted substrate the temperature gradient is steep, and Fig. 7 (e) shows that this temperature increase is only relevant above 0.6 GPa, as there is not enough friction energy to generate sufficiently high temperatures at lower pressures.

Figure 7 (c) shows the normalized contact area $A_{c}/A_{\mathrm{nom}}$ over the normal pressure $\sigma_{z}$, where $A_{\mathrm{nom}}=3595$ nm2 is the lateral cross-section of the simulation box. $A_{c}/A_{\mathrm{nom}}$ can be considered the degree to which the two counterbodies are in contact. Since this contact is not necessarily flat (in contrast to $A_{\mathrm{nom}}$), the dimensionless ratio can in principle exceed 1. For the fully thermostatted systems, the normal pressure dependence is basically linear, while for the base-thermostatted systems it exhibits a distinct kink around $\sigma_{z}\approx 0.5$ GPa, separating the normal pressure dependence into two regimes of smaller and larger slope for lower and higher normal pressures, respectively. For the systems with electron-phonon coupled thermostatting, the data lie in between, but markedly closer to those obtained for the fully thermostatted systems, with an apparent change in slope between $\sigma_{z}\approx 0.6-0.75$ GPa.

As can be seen in Fig. 7 (f), it becomes evident that in the base-thermostatted substrates, the average contact temperature increases linearly with $\sigma_{z}$ independent of the heat treatment procedure, with a temperature increase of $\approx 45$ K per 100 MPa. For full and electron-phonon coupled thermostatting the contact temperature remains at a constant 350 K and 410 K, respectively. The temperature difference between the heat bath and the contact zone may therefore only be a function of the abrasive particle geometry, the sliding velocity, and the thermostat coupling constant.

In Fig. 8 we analyze how the shear stress $\sigma_{x}$ and the final wear depth $h_{w}^{\mathrm{(end)}}$ vary with the normalized contact area $A_{c}/A_{\mathrm{nom}}$. In panel (a), the slope corresponds to the effective shear strength $\tau$ of the system. For the fully thermostatted substrates, this quantity is almost a constant, while for the other two thermostatting procedures we observe some more pronounced flattening of the curve, which reflects the onset of thermal substrate softening. Panel (b) shows that the final wear depth can essentially be considered a superlinear function of the contact area, almost independent of the contact temperature and therefore also the thermostatting procedure. We can therefore formulate that the contact area increases due to thermal softening of the substrate, while the separation of atoms from the substrate is eased at higher temperatures. As this correlation of wear depth and contact area depends on the local temperature of the near-surface region immediately beneath the contact spot, the same superlinear trend is observed irrespective of any thermostatting procedure.

4.2 Thermostat influence on microstructural development

The grinding process changes the microstructure as well as the defect structure regardless of thermostatting procedure or normal pressure. These changes are more pronounced at higher normal pressures, and their extent depends critically on the thermostatting procedure. We therefore chose the abrasion simulation carried out at the highest normal pressure (0.9 GPa) to best illustrate the possibly drastic differences between the three thermostatting variants, see Fig. 9. For an explanation on how the tomographic slices correspond to the 3D model, recall Fig. 1. The analysis scheme using computational substrate tomography is explained in detail in ref 42. In the following, individual panels of the substrate tomographs in Figs. 3, 9, and 10 will be referred to by their row number $m$ and their column number $n$, abbreviated as R$m$/C$n$, cf. the bottom left legend in Fig. 9. The defect structure evolving due to dislocations as well as any changes of the grain boundaries can be observed via the centro-symmetry (CS) parameter 63 panels in Fig. 9 (a,c,e). The differentiation between grain growth and recrystallization or lattice rotation is possible via the orientation panels in Fig. 9 (b,d,f) with grain coloring according to inverse pole figure (IPF) standard. For an interpretation of the microstructural changes as a result of thermally activated processes, we refer to the temperature plots resulting from low (0.3 GPa) and high pressure simulations (0.9 GPa) in Fig. 10. In general, we observe grain growth, abrasion of entire grains, and full or partial rotation of the lattice structure. Part of the friction energy is stored in dislocations or can lead to the formation of new grain boundaries.

The base-thermostatted simulation results in the most dramatic changes to the initial polycrystalline structure, compare Fig. 3 (a) with Fig. 9 (b). We observe massive grain growth, and only grain boundaries pinned to the rigid base via the boundary conditions do not experience any changes. The grain structure at the surface is strongly disturbed. As the wear particles cannot cool down by themselves in the base-thermostatted simulations, their temperatures rise beyond 900 K, see Fig. 10 (b). Although the surface can cool down somewhat between two passing abrasives, the temperature in the near-surface zone remains high in the base-thermostatted simulation, never dropping below 600 K. The grain boundary mobility in this zone is therefore much higher than in the differently thermostatted simulations. Due to the high temperatures, the zone featuring modified grain orientations extends to depths of 10–15 nm. The latter can be visualized using IPF standard coloring, which reveals massive changes between initial and final grain orientations and a slight tendency towards a (111) orientation, though the orientation is too weak to refer to it as a preferred texture evolving during grinding.

After grinding at 0.9 GPa, the final near-surface grain size of the electron-phonon coupled substrate is smaller than that of the base-thermostatted substrate. Only beneath the wear particles, the base-thermostatted substrate forms numerous tiny grains down to half the substrate thickness, whereas electron-phonon coupling results in smaller wear particles and hence a reduced grain-refined area, compare R2/C1 in Figs. 9 (a,c).

The IPF orientation analysis features pronounced color shading within the larger grains, which indicates elastic deformations within the grains that are also referred to as residual stresses of type II, or microstresses 64. The fact that these elastic deformations are stable in more grains than in the base-thermostatted substrate can be explained by the much lower temperatures in the electron-phonon coupled substrate. Besides, the surface can cool down to almost 300 K between two passing abrasives, and the wear particles themselves are cool and do not heat up the surface as in the base-thermostatted process. At the end of the simulation there is no tendency towards a preferred (111) orientation for the electron-phonon coupled substrate.

We observe grain growth compared to the initial microstructure for the fully thermostatted simulation as well. The CS parameter analysis of the fully thermostatted simulation in Fig. 9 (e) shows that small grains have been abraded and grain growth has taken place, but only in the first nm beneath the surface. As the temperature of the wear particles and the contact zone only reaches 350 K, any temperature rise due to frictional heat is roughly limited to the upper 3 nm of the substrate, recall Fig. 4. Since grain growth is primarily a thermally activated process, the grain boundary mobility is only high enough in this limited zone of the first nanometers beneath the surface (located at $z\approx 20$ nm). Consequently, the grain sizes at the surface of the fully thermostatted substrate are smaller than in the electron-phonon coupled substrate and much smaller than in the base-thermostatted substrate after 5 ns of grinding, see Fig. 9 (e,f). The grain orientation analysis in panel (f) reveals that grains grow at the expense of small grains at the surface, and that larger grains located in the middle of the substrate are stable with respect to their orientation and only grow towards the surface.

The fully and the base-thermostatted scenarios reflect the two extremes of microstructural changes. As the evolution of the microstructure is primarily driven by the evolving temperature (when comparing equal pressure processes), thermostatting strongly changes the temperature levels and gradients at different depths of the substrate. Electron-phonon coupled thermostatting produces temperatures between the two bounds of the fully and base-thermostatted variants, compare Fig. 10 (c,d) with (a,b) or (e,f). Electron-phonon coupling leads to higher maximum temperature levels in the contact than full thermostatting, but also to elevated temperature fields confined to substrate regions directly beneath the contact area when compared to base thermostatting, cf. Fig. 4. Furthermore, the maximum temperature gradients $dT/dz$ (occurring within the shear zone around $z\approx 20$ nm) are considerably steeper for base thermostatting than for the two other thermostatting modes.

The main microstructural feature that changes during grinding in all simulations is the grain size. Grain growth occurs in part by grain boundary movement and in part by lattice rotation. Occasionally, tiny new grains form in both non-base-thermostatted variants, but they are not sufficiently stable to withstand a passing abrasive.

For the different thermostatting procedures, the resulting temperature gradient within the substrate and the shear zone determines the reaction of the microstructure at the respective normal pressure. The external mechanical loads are fairly equal if the pressure remains fixed. For a comparison of the evolving temperature gradients at low normal pressure (0.3 GPa) for base, electron-phonon coupling, and full thermostatting, see Fig. 10 (a,c,e). As the $T$-gradient in the shear zone is very steep for base thermostatting, even this low normal pressure machining process changes the microstructure markedly. The degree to which the microstructure is modified is thus directly defined by the level and extent of the three-dimensional temperature field.

5 Conclusion and Outlook

In this work we have analyzed the sliding interface between a nanocrystalline ferritic work piece and hard, abrasive particles during a nanomachining process simulated with molecular dynamics simulations. We have shown that choosing a thermostatting procedure that can mimic the electronic contribution to heat conductivity is crucial when large amounts of energy are introduced to metallic systems in a highly localized fashion. In our application, the resulting near-surface microstructural development, surface topography, and quantities such as the wear depth depend strongly on the chosen thermostatting procedure for normal pressures exceeding 0.4 GPa, while the total shear stress is almost unaffected up to a normal pressure of approximately 1 GPa. The correlation between plastification of the substrate and varying normal pressures as well as the evolving dislocation structure will be discussed in detail in a forthcoming publication by some of the authors.

As beneficial as the thermostatting approach described in this work is for simulating surface finishing and wear processes, there are several aspects of it that could be improved further. As the target temperature of the thermostat is constant throughout the simulations, it is currently impossible to reach the practically obtained equilibrium temperatures for dry friction with simulation times restricted to several nanoseconds. A simple workaround could be to set the thermostat temperature close to the theoretically predicted value at that depth in the work piece, but a proper treatment would have to include a coupling of the base thermostat to a far-field solution of the temperature depth profile. This could be obtained using a one-dimensional Green’s function approach 54, 65, assuming an averaged heat flux at the sliding interface and a semi-infinite work piece, thus placing the bulk temperature boundary condition “far away” from the heat source. Properly done, this would allow the system to reach a realistic surface temperature. A second aspect that is not trivial to resolve is that the Langevin thermostat somewhat unrealistically damps all modes of heat conduction equally, which might be improved by employing a more sophisticated thermostat. Finally, since nanoscopic heat conduction depends on the crystal system and orientation, it might be an oversimplification to assume a constant macroscopic heat conductivity for a 20 nm thick surface layer, even if it consists of initially randomly oriented grains. A proper treatment of this issue would include a real-time grain orientation analysis and subsequent adaptive thermostat parametrization, so any benefit gained from such an approach would have to be carefully balanced with the considerable additional computational cost.

Although the thermal and microstructural trends in our results match intuition, the highly desirable experimental validation is unlikely to be exhaustive. All of the thermostatting scenarios in our work except the electron-phonon coupled ones are difficult to reproduce in an experiment, as it is not trivial to suppress electronic heat conduction 66, 67. While one could simply conduct nanomachining trials at higher temperatures, any influence of excessive thermal gradients would still be lost. That said, it may be possible to verify the time-development of the near-surface structural changes occurring in a nanocrystalline sample due to abrasive nanomachining. One aspiring suggestion would be the design of a tribo-experiment to be analyzed in situ using synchrotron X-ray diffraction, similar to what has been done to study the growth of bulk grains in deformed metals 68. Before going to such an effort, it may however be sensible to reconsider the choice of work piece material, e.g., an fcc structure with a low tendency for oxidation while still being technologically relevant.

There exists a second path for, at least qualitative, validation of our simulation results. Aided by ever-increasing computational power, we have recently succeeded in simulating the tribological response of a considerably larger work piece consisting of approximately 25 million atoms. Due to the initial grain size diameters of 30 nm, tribological loading actually leads to dislocation pile-up and near-surface grain refinement, both of which are observed in macroscopic experiments. This bodes well for the future, as it means that phenomena relevant to applied interfaces might be reproduced at a slightly smaller scale, thus uncovering the underlying mechanisms. Ultimately, MD simulation may well evolve into a crucial component of engineering design tools for manufacturing, friction, or wear applications.

{acknowledgement}

This work was partly funded by the Austrian COMET-Program (Project K2, XTribology, no. 849109) and carried out at the “Excellence Centre of Tribology”. SJE, UC-B, and DB acknowledge the support of the Province of Niederösterreich (Project “SaPPS”, WST3-T-8/028-2014). GF acknowledges the financial support from the German Research Foundation (DFG) via SFB 986 “M3”, project A4. The authors wish to thank M. Moseler, M. Müser, J. Belak, A. Vernes, and G. Vorlaufer for helpful discussions about thermal conductivity, electron-phonon coupling, and Green’s functions.