Multiscale study of the dynamic friction coefficient due to asperity plowing

A macroscopically nominal flat surface is rough at the nanoscale level and consists of nanoasperities. Therefore, the frictional properties of the macroscale-level rough surface are determined by the mechanical behaviors of nanoasperity contact pairs under shear. In this work, we first used molecular dynamics simulations to study the non-adhesive shear between single contact pairs. Subsequently, to estimate the friction coefficient of rough surfaces, we implemented the frictional behavior of a single contact pair into a Greenwood-Williamson-type statistical model. By employing the present multiscale approach, we used the size, rate, and orientation effects, which originated from nanoscale dislocation plasticity, to determine the dependence of the macroscale friction coefficient on system parameters, such as the surface roughness, separation, loading velocity, and direction. Our model predicts an unconventional dependence of the friction coefficient on the normal contact load, which has been observed in nanoscale frictional tests. Therefore, this model represents one step toward understanding some of the relevant macroscopic phenomena of surface friction at the nanoscale level.


Introduction
The empirical mathematical summary of Amonton's first (friction force is proportional to the applied normal load) and second (friction force is independent of the apparent contact area) friction laws is the relation F = μN, where F is the friction force, N is the applied normal load, and μ is the coefficient of friction (COF). Even though Amonton's law works well for dry friction problems in traditional engineering, the reason for this remains unclear for quite a long time. In Bowden and Tabor's work [1], Amonton's law was explained by the fact that rough surfaces in contact with each other consist of numerous smaller contact pairs. Consequently, the actual contact area due to these microscopic and plastically deformed contact pairs is much smaller than the apparent contact area. Both experimental and numerical studies have found that the real contact area is actually proportional (or quasi-proportional) to the normal load, that is, A = N/k, where k is a constant whose value depends on the material elasto-plastic property and surface roughness. The friction force can also be interpreted as F = τA, where τ is the shear strength of microcontacts. Therefore, one may interpret the macroscopic friction coefficient by considering the microscopic properties of the asperities as μ = τ/k.
Much attention has been dedicated to studying the linear dependence parameter k between the real contact area A and the normal load N. The GreenwoodWilliamson (GW) statistical model [2] provided a convenient approach that correlates the single asperity mechanical response to rough surface contact properties. A number of studies have attempted to relax the strong assumptions and constraints of the GW model to extend its range of applicability. Examples include using the simplified elliptic model [3] for the asperity shape, considering substrate deformation to include asperity interaction [4,5], extending the model for nearly complete contact [6], or even incorporating size-dependent plasticity [7]. The other influential theoretical approach is the Persson-type contact model [8], which implicitly includes multiple-length scales and solves the contact problem starting from the full contact condition through an analogy to the diffusion problem. Relevant works [914] that are based on the surface fractality also predict a linear or quasi-linear dependence between the real contact area and the normal load. The commonly accepted and Tabor [1] originally claimed that it is the material shear strength. However, the estimations of the COF for metallic surfaces are not consistent with the experimentally measured values. The reason for this may be because the real contact area consists of micro-sized contact pairs, and in particular, the yield strength behaves in a size-dependent manner when an intrinsic length scale is within a range of several micrometers or less. For this reason, discrete dislocation dynamics simulations [16] have been carried out to reveal  for different areas of contact considering the adhesive contact interface. It was found that the contact area and loading rate have a clear effect on the frictional strength,  , which is observed to have three regimes: adhesion controlled, plasticity controlled, or mixed.
It should be noted that the framework discussed above based on k and  only attributes the frictional process to the flattened asperities ( Fig.  1(c)). These ideally flattened asperities generally result from the contact between a rough surface and a rigid platen, and the friction corresponds to static friction. However, in a more realistic frictional model, one rough surface would slide relative to the other rough surface. At the surface asperity level, one would observe many asperity plowing pairs, as shown in Fig. 1(a). The mechanical responses of these plowing asperities, at the microscopic length scales, are main reasons for macroscopic dynamic friction. Asperity plowing ( Fig. 1(b)) has been studied in detail, for example, by performing finite-element method (FEM) simulations [17,18] and discrete dislocation dynamics [19,20] and molecular dynamics (MD) simulations [2124]. Those studies mainly focus on the asperity mechanical response under different loading conditions, such as plowing/interference depth, asperity size, and plasticity properties. Unfortunately, the connection to the surface frictional property was not considered. It can be seen in Fig. 1 combined average can be replaced by a time average over the simulation period if, in the simulation, the fraction of time that the system spends in each state satisfies the Boltzmann distribution. For our friction problem, we assume that the rough surfaces are sufficiently large such that all pair configurations are equally likely to be found. Therefore, we can replace the combined average by the time average of a simulation where all plowing asperity configurations are covered. Such a simulation can be the plowing of a single asperity, as shown in Fig. 1(b), where the deformable asperity (colored in blue) is fixed at the bottom, and the rigid asperity (colored in gray) moves rightwards under a constant velocity V. This simulation will contain, for example, all six configurations shown in Fig. 1(a) because they have the same plowing depth. Then, the average tangential force and normal force can be written as cont cont where cont t represents the contact time, that is, the time duration between the instances the two asperities form contact and lose contact. Using the above idea, Mulvihill et al. [25] studied the time-averaged single asperity responses under different interference/plowing depths using FEM simulations. However, owing to the complex FEM model, which has to include large deformations and material fracture effects, the numerical convergence is difficult and the asperity interference depth in their FEM modeling is quite small. Furthermore, there is also limited plastic deformation in the plowed asperity. In addition, the application of conventional continuum plasticity for surface asperity studies is always debatable as it is well known that for crystalline materials, plasticity becomes size dependent at such small length scales. Therefore, it is necessary to utilize simulation techniques that can capture material elastoplastic properties at a small length scale, such as MD simulations. A significant body of literature exists concerning friction-related studies using MD simulations. For example, for a nanoscratching problem, Zhang and Tanaka [26] revealed that there are generally four distinct regimes of deformation during a diamond-copper sliding system, that is, no-wear, adhering, plowing, and cutting regimes. The dominant regimes are governed by some sliding parameters such as invasion depth, sliding velocity, and surface lubrication conditions. Gunkelmann et al. [27] carried out MD simulations of nanoscratching on iron, and the dislocation microstructure was characterized using continuum field quantities from continuous dislocation dynamics theory [28]. The mechanical response of the material is then linked to the evolution of the dislocation field quantities. For grinding processes, the effect of the sliding velocity on the subsurface damage was investigated by analyzing the dislocation and phase transformation processes [29]. It was proposed that the subsurface damage thickness only slightly increased when the sliding velocity exceeded 180 m/s. With respect to the friction law at the nanometer length scale, Müser et al. [30] established a microscopic theory to explain the molecular origins of friction. MD simulations were then used to verify their theoretical predictions. Mo et al. [31] studied the dependence of the friction force on the applied load and contact area at a nanoscale level, and they demonstrated that the friction force is linearly related to the number of interacting atoms.
To the best of our knowledge, the connection between single asperity plowing responses at the atomistic scale and the rough surface frictional property is missing. Herein, we perform a multiscale study at the macro scale. For the small length scale (i.e., asperity scale), we study in detail the single asperity response during the whole plowing process by performing MD simulations. The effects of interference depth, asperity size, sliding velocity, and crystal orientation are considered. For the large length scale (surface length scale), the single asperity response is utilized to predict the COF during surface sliding through a GW-type statistical model. The remainder of this paper is organized as follows. Section 2 presents the methodology and model description. Next, the MD simulation results of single asperity plowing are discussed in Section 3. Section 4 focuses on the prediction of the COF based on statistical model. A discussion and concluding remarks are given in Section 5. The effects of crystalline orientation and numerical fitting of the asperity responses are discussed in Appendices A and B, respectively.

Methodology and model description
The large-scale atomic/molecular massively parallel simulator [32] is employed for MD simulations. The geometry of the simulation model is illustrated in Fig. 2. A hemispherical asperity with the face-centered cubic copper is located at the center of a copper cubic substrate. The radius R of the hemisphere varies from 5a to 50a, with a representing the lattice constant of The interaction between asperities is simulated by a repulsive potential with a force of magnitude Here, K is a constant that is related to the effective stiffness of the rigid asperity indenter, and is set to 10 eV/Å 3 , R is the radius of the indenter, and r i is the distance from the i-th atom to the center of the rigid asperity ( Fig. 2 only shows its lower half). Adhesion at the contact interface plays an important role during nanoscale contact; however, in this study, we focused on dislocation plasticity; therefore, we excluded the adhesion. The repulsive potential eliminates the effect of interface adhesion and the tangential interaction at the contact interface, which is equivalent to having a frictionless interface. In the simulations, the x, y, and z-axes were initially chosen as the [001], [010], and [001] lattice directions, respectively. The mechanical properties of the asperity also depend on the crystal orientations; a detailed discussion can be found in Appendix A. Periodic boundary conditions are imposed along the x and y directions, while the z direction is non-periodic. The interactions among copper atoms in the asperities are described by the embedded atom method (EAM) potential [33], which has been widely used in the simulation of copper [3437]. Prior to the movement of the rigid asperity, the initial structures are fully relaxed at 0.01 K, with a time step of 0.0015 ps for 150 ps (in total 10 5 time steps). After the sample attains thermal equilibrium, the rigid asperity is moved along the y direction with constant velocity V ranging from 3.6 to 360 m/s (these velocities approximately correspond to 0.01a/ps to 1.0a/ps, with a being the lattice constant) with the different interference depth h (i.e., the asperity overlap). A canonical ensemble NVT (constant number of atoms, volume, and temperature) is imposed where a Nose´-Hoover thermostat is used to maintain the system at a constant average temperature of 0.01 K, so that thermal effects are effectively eliminated. This allows us to focus solely on the effect of the plowing depth and loading velocity on the asperity plasticity. In the analysis, dislocation tracking is performed using the dislocation extraction algorithm (DXA) [38], and Ovito [39] is used to visualize the defect structures.
During the asperity plowing process, the normal (along the z direction) and tangential (along the y direction) interaction forces between the asperities are recorded. The tangential force, opposite to the sliding direction, and the compression load along the z direction are defined as positive. These forces, F y and F z , are averaged over the whole plowing process, as explained in Eq. (1), and they are rewritten as: Here, t 0 and t 1 are the moments when the asperities come in contact and lose contact, respectively. For asperities of different sizes, normalized forces are also introduced as follows: where E is the Young's modulus of copper and is assumed to be 120 GPa in this study.

Results and analysis of a single plowed asperity
In this section, we focus on the single asperity mechanical response under different conditions, such as the interference depth, asperity size, and plowing velocity. The mechanical response of the asperity is then analyzed based on the evolution of the dislocation microstructure and atomic deformation.

Effects of the interference on asperity plowing process
The effect of the interference depth on the mechanical response of the asperity was studied during the single asperity plowing process. The deformable asperity is orientated as x[100], y[010], z [001]. The interference depth h is normalized by the asperity radius R: α = h/R, and the dimensionless overlap α is in the range from 0.05 to 0.8. In Figs. 3(a) and 3(b), the normal and tangential forces are plotted as a function of the plowing distance S y normalized by the asperity size. The plowing velocity is set to 36 m/s. For the asperity of radius 30a, the normal force, which is shown in Fig. 3(a), increases initially with an increasing plowing distance, and then decreases; eventually, it becomes zero when asperities lose contact. In general, a larger interference depth results in a higher peak normal force. For α = 0.5, the dislocation microstructures and asperity deformations at different plowing distances are illustrated in Fig. 3(c). It can be seen that during the plowing process, complex and dense dislocation microstructures are formed, and the asperity is plastically deformed. Tangential forces under different interference are shown in Fig. 3(b), which exhibit similar features except for the one at a small interference depth. This can be seen in Fig. 4(a) for α = 0.05; the tangential force exhibits anti-symmetry across the entire plowing distance. The main reason is that the contact interface is essentially frictionless, and there is not enough plasticity generated in the deformable asperity. The deformation process is mostly elastic, which has an antisymmetric feature for the tangential force. The dislocation microstructures at peak forces are also shown in Fig. 4(b). There are dislocations that are nucleated; however, these dislocations cannot glide into the asperities owing to the image force from the surface. The dislocations escape through the surface soon after they are generated; therefore, deformable asperity is essentially elastic.

Effect of the asperity size
We now discuss the effect of the asperity size on its mechanical behavior under plowing for a plowing depth α = 0.5. As shown in Figs. 5(a) and 5(b), normal and tangential forces during the plowing process increase as the asperity size increases. This is expected because more atoms are in contact with the indenter for larger asperity. The evolution of the forces shows the same trend, namely, as the plowing distance increases, the two forces first increase and then decrease until the asperities are no longer in contact. The normalized forces nz F and ny F , which are defined in Eq. (4), are also shown in Figs. 5(c) and 5(d). It can be seen  that the normal and tangential forces generally exhibit size-dependent behavior; a larger asperity results in a larger normalized force. Furthermore, this effect is more pronounced for the tangential force.
Typical asperity deformation and dislocation microstructures for various asperity sizes after the plowing process are shown in Fig. 6(b). It can be seen that for a very small asperity (R = 5a), after deformation, the asperity is dislocation free. For larger asperities (e.g., R = 30a), complex dislocation microstructures (mainly Shockley partial dislocations) are left in the asperity. The dislocation microstructures in bigger asperities appear to fill the entire asperity. The total length of dislocations can be calculated using DXA [38]. With this, the is the initial volume of the asperity, as shown in Fig. 6(a). For asperities with radius R ≥ 20a, the figure shows that the dislocation density is approximately constant.
Following the general idea schematically shown in Fig. 1, we also calculated the time-averaged asperity strength during plowing: the normal and tangential forces during the plowing process, shown in Fig. 5, are time-averaged following Eq. (3) and then normalized using Eq. (4). Figure 7 shows that the time-averaged asperity strength, including both normal and tangential directions, becomes almost size independent with the increase in the asperity size. The physical interpretation of this phenomenon is that for the same dimensionless plowing depth , the average strength of a "statistical asperity ensemble" (i.e., the average strength of asperities with all possible contact configurations) becomes size independent. If one defines the friction coefficient as the ratio of the tangential strength to the normal strength, in our simulation system, the friction coefficient exhibits size-dependent behaviors when the asperity size is smaller than 7.2 nm (20a). Hereafter, we decide to only consider large asperities because for them, the effect of size is not important, and this enables the study of other features/aspects, such as surface roughness and loading velocity in a cleaner manner.

Effect of the plowing velocity
Subsequently, we focus on the effect of the plowing velocity on the mechanical responses of the asperity. The normal and tangential forces for different velocities are shown in Fig. 8. Two clear features are observed. First, with the increase in the plowing velocity, the forces increase rapidly. For instance, the peak force in the normal direction increases from 0.69 to 90.7 μN when the plowing velocity changes from 7.2 to 360 m/s. Second, for a low plowing velocity (e.g., V = 7.2 m/s), the peak force appears earlier rather than in the middle of the whole process.
The effect of the plowing velocity on the dislocation microstructure evolution is shown in Fig. 9 where the plowing distance y S is normalized by the total plowing displacement S of the overall plowing process. Two clear features can be observed:   (1) With increasing plowing velocity, the dislocation microstructure becomes denser, as confirmed by the evolution of the dislocation density in Fig. 9(a). In single crystals, two main mechanisms are responsible for the stress relaxation in the system through plastic deformation: dislocation nucleation and dislocation propagation; the dominant mechanism depends on the external loading rate [40]. When the external loading velocity is small, the dislocations get nucleated and then propagate, having enough time to relax the stress of the system. As a consequence, the dislocation microstructure is more localized because the activated slip systems are sufficient to relax the stress. However, when the loading rate/velocity is high, nucleated dislocations either do not have enough time to propagate or the plasticity generated by dislocation propagation is not sufficient to relax the high internal stress. In the latter case, further nucleation on other slip systems is triggered, which eventually results in a denser dislocation microstructure [40].
(2) The material pile up at the right side of the asperity grows with increasing plowing velocity. This is caused by massive dislocation nucleation in the asperity and subsequent dislocation annihilation at the surface. The shear stress distribution in the asperity caused by plowing has a high value front ahead of the rigid asperity ( Fig.  3(b) in Ref. [19] for the two dimensional (2D) plane strain case). Under high plowing velocity, because the stress cannot be effectively relaxed, the high-stress region would continuously nucleate dislocations, which in turn continuously annihilate at the free surface and leave a surface step. These surface steps accumulate to form a material pile up, which would eventually turn into wear debris. This phenomenon may also explain the increasing wear rate with increasing sliding velocity [41] for metal surfaces.
The time-averaged normalized forces for different plowing velocities are shown in Fig. 10. It can be seen that the normalized forces can be clearly divided into two regimes distinguished by the different plowing velocities. The distinction between the two regimes is made through a Fig. 9 (a) Dislocation density as a function of plowing distance and (b) dislocation microstructures in the asperity at different plowing velocities when α = 0.5. Color scheme of dislocation type follows that of Fig. 6(b).  are the time-averaged-normalized tangential forces corresponding to a given velocity V and an extremely low plowing velocity V 0 , respectively. In this study, V 0 is taken as 3.6 m/s. The mechanical response is considered to be low-velocity dependent when 20%   . As shown in Fig. 10(a), when the asperity interference is α = 0.5, for the low plowing velocity (e.g., V < ~7.2 m/s), asperity responses exhibit low-velocity dependence. When the plowing velocity exceeds a certain value, which is appropriately 12 m/s, the normalized forces start to show a strong velocity dependence. The critical/transition plowing velocity also depends on the interference depth, as shown in Fig. 10(b), where α = 0.2, and the transition velocity is approximately 50 m/s.

Statistical model for slides of macroscopic rough surfaces
In the previous section, we focused on the mechanical response of a single asperity under different loading conditions. In this section, we aim to utilize the single asperity response to predict the rough surface response. Macroscopic rough surfaces consist of various asperities of different heights, and the surface property, such as COF, is the integrand of all asperity pairs. The GW-type model [2] is introduced here to study the COF. As shown in Fig. 11, we consider rigid wavy rough surface slides over a deformable rough surface whose asperity height follows a Gaussian distribution: Here, σ is the standard deviation of the asperity heights. The distance between the asperity peak on the rigid wavy surface and the mean asperity height of the deformation surface is d, which is determined by the normal force that puts two surfaces into contact before sliding. A larger normal force results in a smaller d.
To obtain the macroscopic contact forces (i.e., the normal and tangential forces), we assume that all asperities deform independently. The plowing/interference depth for an asperity of height z is zd when z > d; therefore, for the macroscopic contact with N 0 asperities on each surface and the total normalized (normalized by ER 2 as shown in Eq. (4)) tangential and normal forces between the contact surfaces during sliding are: By substituting the single asperity response in the above equations, the COF can then be calculated by respectively. Here, we note that our study is multiscale in the sense that the mechanical response on the asperity level is integrated to obtain the mechanical response of the surface level. It can be seen that the statistical model requires a single asperity response to be a function of the asperity height z (alternatively α, since α = (z d)/R); therefore, we have carried out a large number of simulations to fit the functional forms for single asperity responses. Details can be found in Appendix B. Because the single asperity responses also strongly depend on the crystal orientation, we have also studied the orientation effect. Details can be found in Appendix A. Overall, the functional forms are fitted for three different loading velocities (when the crystal orientations were the same), for different crystal orientations (under the same loading velocity). This enables us to study the effect of loading (sliding) velocity and crystal orientation on the friction coefficient using a statistical model. As shown in Fig. 12(a), both the normal and tangential forces depend on the surface separation distance d: a smaller d results in larger forces. It can also be seen that both forces depend on the surface roughness (which is defined as σ/R): a rougher surface will lead to larger forces. Furthermore, in the inset, we show the sensitivity of the forces with respect to d. It can be clearly seen that the tangential force is more sensitive to d than the normal force. In Fig. 12(b), the COFs of different surface roughness are plotted against the change in d between contact surfaces. Under a specific sliding velocity of 36 m/s, the increase in surface roughness results in a larger COF. It is interesting to note that decreasing d will increase the COF, which is consistent with the sensitivity shown in Fig. 12(a) inset. This appears to contradict Amonton's third law and macroscopic experimental measurements, where the COF is believed to be a constant that only depends on the material property and surface roughness. However, it has been observed experimentally [42] that the nanoscale friction coefficient resulting from plasticity depends on the normal force; a larger normal force (i.e., smaller d) results in a higher COF.
It can be seen in Fig. 13(a) that the COF evaluated by the statistical model exhibits velocity dependence, which is consistent with the experimental observation: the COF generally increases with an increasing sliding velocity [43]. At the same time, it is seen that the velocity dependence in our model varies with d. When d is large (i.e., the normal contact force is small), most of the asperity plowing pairs will have small interference. As shown in Fig. 10(b), for small interference, the critical transition velocity is high; therefore, 36 and 7.2 m/s exhibit weak velocity dependence. Similarly, when d is small, there will be more asperity pairs with high values of interference, where the critical transition velocity is low, as shown in Fig. 10(a); therefore, the COF exhibits strong velocity dependence. It should be noted that the velocity-dependent COF here mainly originates from rate-dependent dislocation plasticity, rather than other mechanisms in other frictional systems. However, the friction coefficient in our system can also be described well using the commonly used rate and state friction framework, whose general form is [44,45]. In this framework, , which is difficult to observe in the main figure of (a). The ratios show the unbiased sensitivity. the velocity dependence is described by parameter (ab), which needs to be measured experimentally for the specific testing system. In our model, the parameter (ab) is determined by d. Figure 13(b) shows the dependence of the COF on the crystal orientation. The results obtained confirm the estimation made based on the single asperity response in Appendix A; the friction coefficient is smaller along the atom close-packed direction than of those in other directions.

Conclusions
In this study, we used a multiscale model for obtaining the macro-frictional behavior of rough surfaces from the plowing response of a single asperity at the nanoscale.
MD simulations are performed to investigate the non-adhesive plowing of a single asperity. The effects of interference depth, asperity size, relative plowing velocity, and crystal orientation are discussed as follows: (1) the friction forces (frictional strength) are size dependent, and the friction coefficient increases with increasing asperity size. However, it becomes insensitive to the asperity size when the asperity is larger than a critical size (radius R > ~7.2 nm); (2) we find a critical plowing velocity below which the asperity responses exhibit weak velocity dependence; (3) when the plowing direction is parallel to the crystal close-packed plane, it is easy to plow the asperity.
Using the statistical model, we studied the frictional behavior of rough surfaces. The salient conclusions are as follows: 1) Smaller surface separation d yields a higher COF, which contradicts Amonton's law, but is consistent with nanoscale experimental observations.
2) The COF increases with the increasing surface roughness. When the surface is rougher, the COF is more sensitive to d.
3) A higher sliding velocity results in a higher COF. At the same time, the dependence on the sliding velocity strongly relies on the surface separation distance d.
4) The crystal orientation has a clear influence on the COF; when the atomic close-packed surface is parallel to the sliding direction, the COF is smaller.

Appendix
The data shown in the figures along with the python plot scripts are available at https:// gitlab.com/computational-materials-science/publi c/publication-data-and-code/2020_Hu-et-al__Mult iscaleStudyOfTheDynamicFrictionCoefficient

Appendix A: Effect of crystalline orientation
Owing to the anisotropic mechanical behavior of crystalline materials at the atomic scale, the frictional properties of the material also strongly depend on the crystal orientation [4648]. In this appendix, asperities with four different crystal orientations, as listed in Table A1, were studied to reveal the underlying orientation dependence. The asperity size is chosen as 30a (i.e., R = ~10.8 nm) and the plowing velocity is 36 m/s along the y direction. The dimensionless interference depth α ranges from 0.1 to 0.8 (only two typical overlaps are shown in Fig. A1). With increasing overlap, the mechanical response of the asperities exhibits significant orientation dependence. For asperity contact with a large overlap (i.e., α = 0.5), the tangential forces in the asperities with z   As shown in Fig. A2, the distributions of strain yz  as well as the associated dislocation microstructures at a specific interference depth α = 0.5 when Sy/S is approximately 0. where dislocations are easier to get nucleated; therefore, these two orientations are expected to have smaller friction coefficients (because it is easier to deform in the tangential direction) when compared to those of other orientations. However, the dislocation structures are also quite different in these two orientations, as shown in Fig. A2(b). Because most of the activated dislocations in Fig. A2

Appendix B: Fitting functions of single asperity mechanical responses
Macroscopic contacts contain numerous micro and nanoasperity pairs at small length scales. When the two surfaces in contact slide relative to each other, asperity pairs of different sizes with different interference depths will move under the same sliding velocities. Thus, in order to obtain the surface response, it is necessary to quantitatively pre-acquire the mechanical responses of asperity pairs. In this appendix, we aim to provide fitting formulas for the asperity mechanical response-based analysis in the above subsections.
As shown in Figs. B1(a) and B1(b), for the asperity of crystal orientation O-1, the time-averaged normalized average forces obtained by MD simulations are fitted in the power-law form as a function of the interference depth when the plowing velocity is 36 m/s. Similarly, the forces between asperities are also fitted for plowing velocities of 7.2 and 108 m/s. The power-law fitting parameters at different plowing velocities are summarized in Table B1. Similarly, the time-averaged normalized average forces in the asperities of the other three orientations at a plowing velocity of 36 m/s are also fitted in the same way, and the fitting parameters are summarized in Table B2.  The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder.