Curvature flows, scaling laws and the geometry of attrition under impacts

Impact induced attrition processes are, beyond being essential models of industrial ore processing, broadly regarded as the key to decipher the provenance of sedimentary particles. Here we establish the first link between microscopic, particle-based models and the mean field theory for these processes. Based on realistic computer simulations of particle-wall collision sequences we first identify the well-known damage and fragmentation energy phases, then we show that the former is split into the abrasion phase with infinite sample lifetime (analogous to Sternberg’s Law) at finite asymptotic mass and the cleavage phase with finite sample lifetime, decreasing as a power law of the impact velocity (analogous to Basquin’s Law). This splitting establishes the link between mean field models (curvature-driven partial differential equations) and particle-based models: only in the abrasion phase does shape evolution emerging in the latter reproduce with startling accuracy the spatio-temporal patterns (two geometric phases) predicted by the former.

Impact induced damage and fragmentation of solids is ubiquitous in nature and plays a crucial role in the evolution of our geological environment: repeated impacts shape particles (sand grains, pebbles, and volcanic rocks) in sediment transport [1][2][3][4][5][6][7][8][9] , affect the production of ash and pyroclast particles in volcanic eruptions 10,11 , and contribute to the generation of atmospheric aerosols with consequences on air pollution and on the global climate 12 . In the Solar system, the size and shape of asteroids and of the particles of planetary rings observed today are the results of a long lasting collisional evolution [13][14][15][16][17] . On planet Mars traces of fluvial evolution of landforms such as pebbles have been discovered 5,18 similar to river beds on Earth. Particle breakage is widely used by the industry in comminution processes of ores and minerals [19][20][21][22][23] , however, it can also be undesired in process and handling engineering due to the resulting degradation of product quality. In these natural processes and industrial applications particles collide both with each other and with hard walls presented by Earth surface (river beds, beaches, and rock walls) or by the components of the process equipment (conveyors, transportation tubes, and containers). A specific area where particle-wall collisions are of utmost importance is the damaging of aircrafts, especially jet engines by impacting hail particles which can cause power reduction and even flame-out of the engine 24 .
Over the past decades, detailed knowledge has been accumulated in geology 1,2,18,25 , physics [26][27][28][29][30][31][32][33] , and engineering 19,23,[34][35][36] on single impact breakage phenomena, however, a comprehensive understanding of low velocity impact sequences responsible for the gradual mass reduction and global rounding of solid particles is still lacking. In the physics literature, the existence of two distinct energy phases has been established; these are called the damage phase and the fragmentation phase. These two phases have not only been demonstrated for brittle materials 30,32,34,37 , and plastics spheres 38 , but also for liquid droplets 39 . Moreover, the same two energy phases have also been reported in the geophysics literature 8 for the collisional attrition of sedimentary particles. In the latter context, global mechanical and geometric understanding of impact induced breakage would be essential to decipher the information hidden in the size and shape of grains and pebbles [1][2][3][4]18 . While research in geology and physics concentrated on single impact phenomena, mathematical research related to the proof of the Poincaré conjecture [40][41][42] led to the study of a class of nonlinear geometric partial differential equations (PDEs) called curvature-driven flows which appear to be the adequate mean-field theory models for the global evolution of pebbles and other particles under a large number of low energy impacts 43,44 . One may target global shape evolution of particles either by extending the physics literature about single breakage to multiple breakage processes or by relying on mean field PDE models. Although the latter are invaluable tools to obtain qualitative • At sufficiently low velocities repeated impacts result in abrasion of the body and lead to a finite asymptotic residual mass; we call this the abrasion energy phase. • Above a first critical velocity, complete destruction is achieved within a finite number of repetitions; we call this the cleavage energy phase.
The third, highest energy phase, occurring above a second critical velocity corresponds to instantaneous fragmentation where cracks span the entire body and the sample rapidly falls apart into a large number of small pieces; we call this the fragmentation energy phase. The transitions between the abrasion, cleavage, and fragmentation phases occur at two well-defined critical velocities analogous to continuous phase transitions.
We establish the link between microscopic physical breakage models and mean-field PDEs in two steps. First, the splitting of the earlier identified damage phase into the abrasion and cleavage phases delineates the range of validity for the latter: the main feature of the now identified abrasion phase is that each impact removes only a small amount of relative mass. As PDE models are based on the limit where the removed relative mass in each collision approaches zero, our study shows that PDEs can be regarded as a mean field approximation of collision-induced attrition in the abrasion energy phase. Second, we identify one key feature of the PDE model in the microscopic simulation: we show that two geometric phases earlier identified in the context of the PDE model clearly emerge inside the abrasion phase in the microscopic breakage model.
Our finding is based on large scale computer simulations which revealed that the evolution of the mass and shape of the solid is governed by scaling laws in terms of the impact velocity. Most notably, in the abrasion phase the shape evolution of the sample is described by a universal scaling form with a power law dependence on the impact velocity predicting infinite sample lifetime at some finite, asymptotic mass, the latter being determined by the energy threshold for the creation of cracks. In the special limit when this threshold approaches zero, our findings reproduce Sternberg's Law 48 , predicting exponential decay (and infinite lifetime) for sedimentary particles undergoing collisional abrasion in fluvial environments. In addition to verify Sternberg's Law for mass evolution, in the energetic abrasion phase we also confirmed the existence of the two earlier observed geometric phases 9,46 , thus our simulations serve as the first direct mechanical confirmation of curvature-driven PDEs as models of impact-driven abrasion processes. In the cleavage phase we find that the sample lifetime decreases as a power law of the impact velocity analogous to the Basquin law 49,50 of sub-critical fracture.

Results
Single impacts: transition from damage to fragmentation. To understand the evolution of solid bodies under repeated collisions with a hard wall, first we focus on single impact events and quantify the resulting mass reduction. We performed numerical measurements by means of computer simulations of a realistic discrete element model (DEM) of body-wall collisions in three dimensions (3D) 51-55 varying the impact velocity v 0 in a broad range. To represent freshly fractured rocks with sharp corners and edges in the initial state of shape evolution, rectangular samples of mildly elongated cubic shape were created with the aspect ratio 1:1.2:1.4 of their shortest c 0 , intermediate b 0 , and longest a 0 sides. This choice is justified by our recent finding that the average shape of fragments is well approximated by a cube for a large diversity of fragmentation processes 45 . The solid was discretized as a random packing of spherical particles, connected by breakable cohesive contacts [51][52][53][54][55] . Parameters of the model were set in such a way that our DEM provides a consistent qualitative and in certain cases quantitative description of the mechanical and fracture properties of the broad class of heterogeneous brittle materials which are abundant in our geological environment 30,[56][57][58][59] (see Methods and the Supplementary Information for details of the model construction and parameter settings). Initially, the discretized sample was placed close to a planar wall with a random orientation chosen uniformly on the sphere and the impact was initiated by assigning identical velocity v 0 perpendicular to the wall to all particles of the solid . As the body moved, it got into contact with the wall and deformed which could result in cracking and fragment formation. The impact lasted until complete rebound was achieved where all particles separated from the wall. In the final state of the process, particles connected by the surviving cohesive elements were identified as fragments. A snapshot of the impact process is presented in Fig. 1 www.nature.com/scientificreports/ Simulations revealed that for sufficiently low impact velocities v 0 < v a the sample solely underwent deformation around the impact site and rebounded elastically without suffering any damage. Cracks first occurred when v 0 surpassed a threshold velocity v a determined by the strength of the internal cohesive elements of the material. In this low velocity range, deformation and crack formation is restricted to the vicinity of the contact zone, while for high impact velocities cracks can span the entire sample giving rise to rapid breakup. To give a quantitative characterization of the final outcome and the degree of destruction caused by impacts, we determined the average masses M 1st and M 2nd of the largest and second largest fragments, respectively. After normalizing these values by the total mass M 0 we plotted m 1st = �M 1st /M 0 � , m 2nd = �M 2nd /M 0 � as function of the impact velocity v 0 . It can be observed in Fig. 2 that at low impact velocities we have m 2nd ≪ m 1st , i.e. the second largest fragment is orders of magnitude smaller than the largest one, showing that only small pieces are removed from the body around the impact site. This is characteristic for the damage energy phase. Fragmentation is achieved when the second largest piece becomes comparable to the largest one, which first occurs at the maximum of m 2nd defining the critical velocity v f of fragmentation. Beyond the critical fragmentation velocity v f both m 1st and m 2nd decrease monotonically. Figure 2 shows that, depending on the velocity, single impacts give rise either to damage or fragmentation of the sample with a sharp transition at the critical velocity v f . The damage-fragmentation transition   has already been studied in experiments and computer simulations of impacting spherical samples against a hard wall, using heterogeneous brittle materials 30,32,34,37 , plastics spheres 38 , and liquid droplets 39 . In these studies, the same qualitative behavior was obtained for the largest fragment masses m 1st , m 2nd as in Fig. 2, which implies that the overall outcome of the process in the high velocity range is entirely controlled by the impact velocity and its critical value v f , whereas neither the sample's shape nor materials' features have any relevant effect. The detailed analysis of the mass distribution of fragments revealed that the observed universality is caused by the underlying continuous phase transition from damage to fragmentation as the impact velocity is varied 26,39,60 .
The identification of the known damage and fragmentation phases also serves as a verification of our model.

Repeated impacts and the two sub-phases of damage: abrasion and cleavage.
In the previous subsection, confirming earlier results, we established for single impact phenomena the existence of the two main energy phases. Now we will show that, if we consider not just a single impact but impact sequences, the damage phase can be subdivided into two narrower energy phases: abrasion and cleavage. The damage phase, characterized by v 0 ≪ v f is often observed in natural and industrial processes at lower energy levels. Under such conditions, the large residue of the sample typically undergoes repeated collisions which give rise to a complex evolution of its size and shape. In the following we extend the global phase diagram of Fig. 2 refining the structure of the damage phase by characterizing qualitatively different evolution histories of residues under repeated sub-critical impacts. To simulate sequences of particle-wall collisions, in the final state of an impact event we identified the largest fragment as the residue of the body, which was further processed to obtain a completely relaxed object for the initial state of the next impact (see "Methods" for the details of the preparation of the residue). The residue was subsequently randomly rotated in three dimensions and was impacted against the wall with the same impact velocity v 0 < v f as before. The above procedure was repeated up to N max = 400 times, or until complete destruction of the body, at ≈ 60 different impact velocities, respectively. As an example, Fig. 3 illustrates the 8th impact of a residue. For each sequence, 120 different initial samples were used, while in subsequent impacts the residues were randomly rotated by uniformly choosing a direction on the sphere. These calculations revealed an astonishingly rich phase structure of the sub-critical v 0 < v f regime.
To quantify the gradual mass reduction during the collision sequence, Fig. 4a presents the average m r = �M r /M 0 � of the residual mass M r normalized by the initial mass of the sample M 0 , as a function of the impact number N for several values of v 0 . We remark that for a single impact event we have m r ≡ m 1st . At very low velocities v 0 ≪ v f , a single impact always gives rise only to a few fragments which are typically single spheres, i.e. powder in the model. As a consequence, in Fig. 4a the residual mass m r gradually decreases with increasing impact number N, however, mass reduction gets limited for high N values and a finite asymptotic residual mass emerges m r → m a r as N → ∞ . The reason is that due to the decreasing mass M r , the kinetic energy E 0 = 1 2 M r v 2 0 , imparted to the sample decreases, since the impact velocity v 0 is fixed. Consequently, beyond a certain impact number, i.e. below a certain value of M r , the emerging deformation is not sufficient to induce further cracking.
Since only small pieces are removed in single impacts, we term this velocity regime as the abrasion phase of the system characterized by the existence of a finite asymptotic residual mass m a r > 0 . It can be observed in Fig. 4a that the value of m a r decreases with increasing impact velocity v 0 . The value of m a r depends also on the energy threshold for the creation of cracks, i.e. on the strength of cohesive contacts. In the limit when this threshold approaches zero, our findings reproduce Sternberg's Law 48 , predicting exponential decay to zero mass and infinite lifetime for sedimentary particles undergoing collisional abrasion in fluvial environments.
When v 0 gets sufficiently high, the functional form of m r (N) qualitatively changes: the mass of the residue sets to a rapid decrease with N, and repeated impacts give rise to a complete destruction of the sample within a finite number of repetitions. This behavior is characterized by impact velocities in the range v c < v 0 < v f and we call this interval the cleavage phase of the impact sequence. The critical velocity v c of cleavage is defined as www.nature.com/scientificreports/ the threshold velocity above which the asymptotic residual mass is zero even at finite energy threshold for the creation of cracks. In our discrete element model, a complete destruction of the sample is reached when the largest fragment comprises solely a single particle of the discretization. For real materials this state is realized when the residual size approaches a characteristic length scale of the meso-structure, e.g. grain size of materials. The transition from abrasion to cleavage at the critical velocity v c is driven by the changing mechanism of cracking. In the abrasion phase the dominating mechanism of mass removal is chipping, i.e. crack formation parallel to the contact surface with the wall, which leads to the formation of tiny fragments 61,62 . However, in the case of cleavage, cracks penetrate the solid to significantly deeper regions so that a combination of contact damage and fracture occurs, giving rise to coarser products as well. Additionally, the elastic waves generated by the collision give rise to the gradual accumulation of damage inside the residue which, in turn, can result in fatigue crack growth as the impact sequence proceeds 63 .
Our results demonstrate that above the threshold velocity of micro-cracking v a , impact attrition phenomena have additionally two well-defined critical impact velocities v c and v f , which separate the three phases of abrasion, cleavage, and fragmentation with distinct qualitative behaviors. The phase diagram of Fig. 2 provides an overview of the distinct qualitative behaviours of impacting solids. For our model solid, the threshold velocities of abrasion and cleavage are v a /v f = 0.124 ± 0.004 and v c /v f = 0.224 ± 0.005 , with respect to the fragmentation critical velocity v f . Figure 4b demonstrates that rescaling the impact number N with a proper power α of v 0 , curves belonging to different impact velocities v 0 can be collapsed on the top of each other, yielding the scaling form where the scaling function m r (x) can be approximated by an exponential m r (x) ∼ exp (−x) (see Fig. 4b), reproducing the time evolution predicted by Sternberg's law 48 . Best collapse is achieved in Fig. 4b with the exponent α = 2.1 ± 0.15.

Sternberg's law and Basquin's law.
It follows from the scaling analysis that increasing the impact velocity v 0 the characteristic impact number N c of the time evolution decreases as a power law The scaling law Eq. (2) holds in both the abrasion and cleavage phases v a < v 0 < v f of impact attrition. For cleavage, the characteristic impact number N c can be interpreted as the lifetime of the sample. Since the peak stress, emerging at the contact zone during impact, increases as a power of the impact velocity v 0 64 , it follows that the expression (2) of residual lifetime is analogous to the Basquin law of sub-critical fracture phenomena 49,50,[65][66][67] . The Basquin law of fatigue life is a fundamental law of sub-critical fracture. It expresses that under a constant or varying sub-critical load, where the stress amplitude falls below the fracture strength of materials, failure occurs in a finite time which decreases as a power law of the externally applied stress amplitude 49 . Our results demonstrate that the Basquin law holds also for sub-critical impact phenomena.
In the abrasion phase N c characterizes the rate of convergence to the asymptotic residual mass m a r . Additionally, the impact velocity also determines the value of m a r , which tends to zero when approaching the critical point Shape evolution: geometric phases inside the abrasion energy phase. Mean field models. In case of polyhedral initial samples, in the abrasion phase we expect that at the beginning of the impact sequence sharp corners and edges are gradually removed, giving rise to an evolution towards an asymptotic rounded shape. In the cleavage phase, due to the breaking of coarser pieces this evolution is more erratic and eventually results in an ultimate destruction. Due to the small size of fragments, we expect that geometric aspects of the abrasion phase may be well reflected in the solutions of averaged, mean field geometric PDE models of attrition 43,44 . The simplest, two-dimensional version of these PDE models may be written as where V denotes the speed by which a surface points moves inward along the surface normal, κ is the scalar curvature and Eq. (4) is often referred to 40,41 as the curve shortening flow, or as a geometric heat equation, referring to a class of geometrically defined, parabolic partial differential equations which have been, ever since the groundbreaking work of Mullins 68,69 , broadly used in surface evolution models 43,70 . The constant c can be regarded as scaling of time and plays no role if evolution is plotted as a function of the normalized residual mass m r . From the mathematical perspective, if we restrict ourselves to classical solutions of Eq. (4) then, instead of considering polyhedra as initial data, we rather consider smooth approximating sequences from which we can pick initial conditions arbitrarily close to a polyhedron. In the approximating sequence all first and second derivatives are continuous.
We remark that Eq. (4) is written in a compact, invariant notation, details about this and other notations are given in Section 1 of the Supplementary Information. The 3D version of this impact-induced attrition model, called the Gauss curvature flow was first introduced by Firey 43 and its convergence to the sphere was ultimately proven by Andrews 71 . One advantage of the compact notation of Eq. (4) is that the 3D version can be described by an analogous formula, for details see Supplementary Information. Next we will show that these expectations are well founded and PDE models serve indeed as good approximations of impact-induced attrition processes, however, only in the abrasion phase.
Shape descriptors. To give a quantitative characterization of the rounding process, we picked three dimensionless descriptors of the overall shape of the residue which not only provide efficient monitoring of the geometric evolution but also admit meaningful comparison with earlier results: axis ratios, circularity (isoperimetric ratio) and intact surface ratio. www.nature.com/scientificreports/ 1. Axis ratios c/a and b/a are traditional geological descriptors 9 characterizing the shape of the residue 25,72 where a > b > c refer to the axes of the bounding box of the residue, aligned with the edges of the initial (cuboid) sample. 2. Isoperimetric ratio or circularity of a planar object is given as R = 4πA/P 2 , where A, P refer to area and perimeter, respectively. It has been observed 8 that circularity of the largest projection of sedimentary particles shows universal features in fluvial abrasion and its evolution is entirely determined by the mass loss during impact induced attrition processes. In our DEM, A and P of the residue were obtained as the area and perimeter of the convex hull of the point cloud of the largest projection of the spherical particles of the relaxed body. For more details on shape descriptors see Subsection 1.3 of the Supplementary Information. 3. Intact surface ratio S/S 0 , expressing the intact fraction of the initial surface, was selected following an idea of Richard Hamilton 46 who, in one of the papers dedicated to the study of curvature-driven flows (leading ultimately towards to his seminal contribution to the proof of the Poincaré -conjecture) describes a curious nonlinear phenomenon about intact surface ratio in the Gauss curvature flow which is the 3D version of (4): he predicted that S/S 0 will drop to zero after a finite time, marking the end of the first geometric phase for cuboids.Hamilton's result inspired further, detailed research on other curvature-driven PDEs 73 which found that whether or not flat sides are preserved depends on delicate features of these models. (For more details see Subsection 1.4 of the Supplementary Information.) In the initial state of DEM samples S 0 is determined as the number of particles covering the external body surface, then the surviving intact surface S is obtained by tracing the particles removed from the initial surface S 0 in subsequent impacts.
The evolution of axis ratios c/a and b/a and the evolution of the isoperimetric ratio R has been computed in the PDE model 9 for the very same cuboid initial conditions as in our DEM study. For the evolution of intact surface ratio S/S 0 in the PDE model we have an analytical result 46 . We will now establish the link between PDE models and microscopic computations by comparing these evolutions. The most striking qualitative feature of the PDE model is the spontaneous emergence of two geometric phases and our computations reveal that these phases are perfectly captured in the microscopic DEM approach. To make the comparison between plots for shape descriptors meaningful, next we seek the corresponding scaling laws.

Scaling laws.
Increasing v 0 accelerates mass removal and thus shape evolution. Figure 6a, c demonstrates that both axis ratios c/a (N, v 0 ) , b/a (N, v 0 ) remain initially constant, display sudden growth between the characteristic impact numbers N r and N s and subsequently saturate. Both the overall shape of these functions and their saturation values remain the same in the entire abrasion phase, however both N r and N s decrease with increasing v 0 . We found that rescaling these curves with v γ 0 , they collapse onto master curves (see Fig. 6b, d) implying the scaling structure where �(x) and �(x) denote the scaling functions. This also implies that N r and N s both have the same power law dependence where the exponent γ was obtained numerically γ = 3.0 ± 0.07 . The saturation values �c/a� ≈ 0.865 and �b/a� ≈ 0.925 show that the asymptotic stable shape of the object is slightly anisotropic which may be a consequence of the finite number of the non-breakable discrete elements in the simulation. Our simulations revealed that under the condition of isotropic impacts, the origin of the universal scaling forms is that the shape of the evolving object is controlled by the total relative mass µ(N) = 1 − m r lost in N repeated collisions. Recently, it has been suggested 8 that µ(N) is also controlling the evolution of the circularity R, so henceforth we use this representation for all shape descriptors.
Geometric phases. The PDE model (4) predicts for the evolution of cuboid blocs with moderate initial axis ratios the emergence of two geometric phases: in phase 1 axis ratios c/a, b/a remain approximately constant while roundness increases steeply and saturates close to 1. In phase 2 the opposite happens: axis ratios increase steeply and saturate close to 1 while roundness remains constant. The conceptual plot of this evolution (as a function of the relative abraded mass µ = 1 − m r ) is shown in Fig. 7b1, accompanied by conceptual contours of the specimen, projected along the shortest (c) axis (b2) and representative snapshots of DEM simulations (b3). Figure 7c1 shows the same plot for b/a and R, obtained from the numerical computation 9 of the PDE (4). Figure 7(c2) presents the hand-drawn sketch of Hamilton 46 of his analytical result on the same PDE: intact surface area S/S 0 survives for a finite time and this marks geometric phase 1.
In Fig. 7a we compare the DEM computations to the aforementioned analytical predictions. In Fig. 7a1 we show evolutions of the average axis ratio b/a and roundness R in the abrasion energy phase v a < v 0 < v c . Note that curves of different impact velocities all fall on the top of each other in agreement with the scaling collapse predicted in the previous section. It is apparent that we have good qualitative agreement with Fig. 7c1: b/a remains constant at the initial value �b/a� = 1.2/1.4 until µ * ≈ 0.34 while R increases sharply and the opposite can be observed for µ > 0.34 . Based on this observations we can clearly record the presence of the two geometric phases in the abrasion energy phase for the evolutions of the axis ratios and the roundness. www.nature.com/scientificreports/ In Fig. 7a2 we show the evolution of the intact surface ratio S/S 0 , also in the abrasion energy phase v a < v 0 < v c . We can observe that this shape descriptor drops to zero at the same relative abraded mass value ( µ * ≈ 0.34 ) which separates the two phases for the evolution of axis ratios and roundness. This is in agreement with the prediction of Hamilton 46 who claimed that intact surface area will survive for a finite time. It is easy to see that as long as intact surface area exists, the corresponding axis ratio of the cuboid (computed from the bounding box) will remain constant so here again we see a perfect match between the DEM computations and the prediction based on the PDE. The transition point µ * between the two geometric phases in the microscopic DEM and macroscopic mean field PDE descriptions of shape evolution have a very good agreement.
This confirms our claim that in the abrasion energy phase v a < v 0 < v c the PDE model offers adequate description of the shape evolution. In sharp contrast, Fig. 7a3 illustrates the evolution of the axis ratio b/a and roundness R in the cleavage energy phase v c < v 0 < v f , both displaying a non-smooth behavior: here we do not expect any mean-field PDE model to provide an adequate description.

Discussion
Impact induced attrition processes cover a broad variety of phenomena ranging from the gentle removal of fine powder from the surface of rock pieces by low velocity impacts to the immediate disruption of objects in energetic collisions. Understanding gradual mass removal due to a sequence of impact events is crucial in sedimentology since pebbles can be considered as witnesses of the geological conditions of their creation. Universal scaling laws of lifetime, size, and shape of evolving particles are indispensable to decode the information imprinted in pebbles [1][2][3][4]18 . In the initial state of this evolution process freshly fragmented rocks are generated 25 by dynamic breakup of rock masses due to high velocity impacts. While the theory of single impact of solid particles with a hard wall is well understood at the level of particle-based models, impact sequences have been so far only www.nature.com/scientificreports/ modeled by mean field theory which necessarily included gross simplifications of the breaking process. Here we offered the first link between particle-based models and mean field theory for collision sequences.
The main methodological novelty of our study is that we use the discrete element method to realistically simulate the entire physical process of all the individual impacts of long sequences without any additional assumption. Although at high computational costs (by simulating ≈ 5 × 10 5 collisions with samples consisting of ≈ 12,000 discrete elements), this approach enabled us to unveil the rich phase structure of impact induced attrition processes. Based on experimental observations, a descriptive classification of single impact breakage has been proposed in 19 , where low, intermediate, and high velocity ranges were distinguished according to the amount and structure of the resulting damage of the body. Here we demonstrated that in multiple impact processes these regimes are separated by universal phase transitions. In addition to the already known damage and fragmentation phases (separated by the critical impact velocity v f ) we identified the abrasion and cleavage phases inside the damage phase (separated by the critical velocity v c ). Abrasion results in finite asymptotic mass (analogous to Sternberg's Law 48 ) while cleavage results in a complete destruction after a finite number of impacts, with sample lifetime decreasing as a power law of the impact velocity (analogously to Basquin's law).
By identifying the abrasion energy phase we were able to provide the link between microscopic, particlebased models and mean-field curvature-driven equations. We showed that the latter can be regarded as adequate approximations of the former, however, only in the abrasion phase. Our simulations revealed an astonishing universality of the evolution of rounding of the residue. Both the axis ratios and the circularity of the largest projection proved to be entirely determined by the attrition mass: evolutions at different impact velocities v 0 can be collapsed onto a single curve by rescaling the number of impacts with a proper power α (also called the lifetime exponent) of v 0 . This universality confirmed earlier conjectures and observations 8,9 on the existence of two geometric phases and also helped to identify a scaling law of the dynamics: the characteristic event number of the onset of shrinking of initially angular objects proved to decrease as a power law of the impact velocity. We were also able to verify a curious effect of geometric nonlinearity, first predicted by Hamilton 46 : in case of polyhedral initial shapes, a finite amount of the initial surface area survived abrasion for a finite amount of time.
Our findings also fit into the broader picture of efforts to approximate PDE models by microscopic, particlebased simulations. In the context of curvature-driven surface evolution, closest to our current topic, Monte-Carlo  www.nature.com/scientificreports/ simulations of the Kardar-Parisi-Zhang (KPZ) equation proved to be a powerful tool to understand the global dynamics 74,75 . However, in contrast to our approach, discrete KPZ models do not use a mechanics-based DEM Kernel and most often they are aimed at surface growth in an orthogonal [xyz] frame. It is important to emphasize that the excellent qualitative and quantitative agreement (e.g. for the transition point between the two geometric phases) of the microscopic DEM and macroscopic PDE descriptions of shape evolution were obtained without any parameter tuning of DEM simulations. This confirms the high degree of robustness of the results for the broad class of heterogeneous brittle materials. For the initial state of shape evolution we considered mildly anisotropic cuboids, since it has proven to be the generic average shape of freshly fractured rocks 45 . Cuboids with other axis ratios would only change the time scale of shape evolution and shift the transition point µ * between the geometric phases. Inside the energy phases of abrasion and cleavage, the temporal evolution of mass and shape is controlled by the impact velocity which we could cast into scaling laws. The value of the scaling exponent of lifetime (cleavage) α falls close to 2, while the exponent γ controlling the shape evolution (abrasion) has a higher value γ ≈ 3 . Based on fracture mechanics, approximate analytical expressions have been derived for the threshold velocities of the onset of abrasion v a and fragmentation v f 62 . These calculations showed that the critical velocities separating the energy phases of impact attrition phenomena depend on material properties as well as on the mass and linear extension of the sample 62 . Based on the analogy to continuous phase transitions and on the good quantitative agreement between the PDE and DEM approaches for the transition point µ * , we conjecture that the critical exponents α , β , and γ are universal, they depend neither on mechanical, nor on geometrical features of the system.

Methods
Discrete element model. We performed computer simulations of the repeated sub-critical impact of solid bodies against a hard wall in the framework of a discrete element model of heterogeneous brittle materials which has been successfully applied before to investigate fracture and fragmentation under various types of loading conditions [51][52][53][54] . In the model the sample is represented as a random packing, consisting, on the average, of 12,000 spherical particles with a uniformly distributed diameter d in a narrow interval d around the average d with �d/ �d� = 0.05 53 . The initial packing is generated by sedimenting the randomly sized particles in a rectangular container which provides a high quality representation of the disordered isotropic micro-structure of rocks. Cohesive interaction is realized by beam elements which connect the particles along the edges of Delaunay triangles constructed from the initial particle positions. In three dimensions (3D) the total deformation of a beam is calculated as the superposition of elongation, torsion, as well as, bending and shearing 76 . Cracks are formed when overstressed beams break according to a physical breaking rule. The breaking condition takes into account the stretching and shearing of particles contacts. The interaction of contacting particles which are not connected by beams is described by the Hertz contact law 76 . In the model the random packing of particles is the only source of disorder which determines the physical properties of beams such as length, cross section, elastic moduli, and moments, as well. At the broken beams along the surface of the spheres cracks are generated inside the solid and as a result of the successive beam breaking fragments are formed. The time evolution of the fragmenting solid is obtained by solving the equations of motion of the individual particles with proper initial and boundary conditions. The model has been validated by comparing (i) the stress field generated in body wall collisions to finite element calculations 30,56 , and (ii) the crack structure and fragment mass distributions to the experimental findings [57][58][59] . Further details of the model construction and of the parameter setting can be found in the Supplementary Information and in 53,55 . Simulations of body-wall collisions. We used the model to simulate the impact of a rectangular body with a hard wall. The initial state of the simulations was prepared by placing the cubic sample close to a hard wall with a random orientation assigning the same initial velocity to the particles pointing perpendicular to the wall. As the body moves, its particles overlap with the wall and experience an elastic restoring force according to the Hertz contact law 64,76 giving rise to deformation and cracking of the body.
Preparation of the residue for repeated impacts. In the final state of an impact process particles of the fragments are not completely relaxed in the sense that the fragments can be deformed and would gradually relax by dissipating energy due to the internal friction of the material captured as a viscous damping force between contacting particles. In order to reduce the computational time, we identify the particles of the residue inside the original sample and replace it by its relaxed counterpart. For each collision during the sequence, the residue is randomly rotated in the initial state to avoid any directional dependence. At each impact number, averages over 120 samples, i.e. over 120 impact directions, were performed, which ensured the high quality of the results even for small residual sizes where the surviving particle clusters have an irregular shape.

Data availability
Mass data and shape descriptors of residues obtained from DEM simulations and PDE computations are freely available in the OSF Data repository at https:// osf. io/ g2ftd/.

Code availability
The numerical code used for data evaluation in this paper is available from the corresponding author upon reasonable request.