FACTA UNIVERSITATIS Series: Mechanical Engineering Vol. 19, No 1, 2021, pp. 7 - 22 https://doi.org/10.22190/FUME201221012S © 2021 by University of Niš, Serbia | Creative Commons License: CC BY-NC-ND Original scientific paper A DISCRETE ELEMENT FORMALISM FOR MODELLING WEAR PARTICLE FORMATION IN CONTACT BETWEEN SLIDING METALS Evgeny V. Shilko1,2, Aleksandr S. Grigoriev1, Alexey Yu. Smolin1,2 1Institute of Strength Physics and Materials Science SB RAS, Russia 2National Research Tomsk State University, Russia Abstract. The paper describes an advanced discrete-element based mechanical model, which allows modelling contact interaction of ductile materials with taking into account fracture and surface adhesion by the cold welding mechanism. The model describes these competitive processes from a unified standpoint and uses plastic work of deformation as a criterion of both local fracture and chemical bonding of surfaces in contact spots. Using this model, we carried out a preliminary study of the formation of wear particles and wedges during the friction of rough metal surfaces and the influence of the type of forming third body (interfacial) elements on the dynamics of the friction coefficient. The qualitative difference of friction dynamics in the areas of the contact zone characterized by different degrees of mechanical confinement is shown. Key words: Dry friction, Cold welding, Wear debris, Mechanical confinement, Computer modeling, Discrete element method 1. INTRODUCTION Although adhesive wear during sliding friction of dry surfaces has been intensively studied over the past several decades, it is still a phenomenon extremely difficult to understand and predict [1-5]. The most popular approach to understanding and interpretation of key aspects of friction and wear is the concept of the third body proposed by M. Godet [6- 8]. This concept implies the formation and development of material fragments (interfacial elements) between contacting surfaces and the determining contribution of these elements to the dynamics of contact interaction including the value of the friction coefficient and wear rate. In friction pairs formed by ceramic materials (normally characterized by low surface energy), these interfacial elements are usually considered as wear particles [9, 10], while in metallic friction pairs the situation is much more complicated. Due to the high values of Received December 21, 2020 / Accepted February 01, 2021 Corresponding author: Evgeny V. Shilko Institute of Strength Physics and Materials Science, pr. Akademicheskii 2/4, Tomsk 634055, Russia E-mail: shilko@ispms.tsc.ru 8 E.V. SHILKO, A.S. GRIGORIEV, A.YU. SMOLIN surface energy and ductility of metallic materials, an important role in the formation and evolution of the third body is played by the adhesion of interfacial elements to contact surfaces and to each other. Adhesive strength can be comparable to the cohesive strength of materials themselves. Therefore, interfacial elements, being the main load carriers, undergo large plastic deformations and are hardened stronger than the surface layers, and so become strong and stable formations. In the process of friction, these formations can increase in size by accretion of smaller interfacial elements and adherence of material from the surfaces, as well as occasionally crush [11-13]. Accordingly, the characteristic dimensions of the main load-carrying interfacial elements may exceed the height of surface asperities (for both initial and current surfaces) by orders of magnitude. Therefore, they completely separate the interacting surfaces and form the actual relief of the contact zone (defined as a set of spots of local contact of surfaces with interfacial elements) [5]. The roughness characteristics of the surfaces themselves at the stage of steady-state friction have only an indirect effect on the real contact area and its change, and, consequently, on the dynamics of friction. Interfacial elements in metallic friction pairs can have not only significantly different sizes but also geometries. The latter is determined by differences in the conditions and mechanisms of the formation of third body fragments. Common names for the largest interfacial elements are prow, wedge, flake (these are generally transfer fragments, but not wear particles), and wear particles/debris themselves [14, 15]. The conditions and mechanisms for the formation of these interfacial elements (and their functional roles) are widely discussed in the literature, but there is still no unambiguous understanding of many key issues. In particular, a recent discussion by J.A. Greenwood [5] raises again the question of how wear debris forms in metallic friction pairs and how they relate to wedges/prows. An effective tool for a deeper insight into the physics (in particular, mechanics) of the formation and evolution of prows, wedges and wear particles is computer modelling [16, 17]. Macroscopic regularities of friction and wear are traditionally studied using FEM or BEM with implemented models of contact interaction and wear [18-22]. Such studies, however, do not come close to understanding the mechanisms of the formation of wear particles. The most attractive tool for performing such studies is the molecular dynamics method. In recent years, MD modelling helped to obtain important information that shed light on the conditions and mechanisms of the formation of nanoscale transfer fragments, which can subsequently contribute to the formation and accretion of wedges, prows, or wear particles [23-28]. However, due to the extreme smallness of the studied spatial and temporal scales (the contact of single nanoasperities is usually considered) and the use of the ideal crystal structure of the initial surface layers, these studies do not allow making unambiguous conclusions about the development of the studied processes at larger spatial scales and times. In recent years, various authors have proposed mesoscale numerical models that allow describing the processes of ploughing, lump/wedge formation, debris particle growth [12, 29-31]. However, these models consider single asperities, while the evolution of the contact zone (including separation of the surfaces and wear) is determined not simply by the independent evolution of single interfacial elements, but also by the mutual influence and competition between interfacial elements even located far from each other in the contact zone. Particle-based numerical method of discrete elements (DEM) is known as a promising and potentially effective tool for the numerical study of complex processes of formation of transfer fragments and small wear particles, their subsequent coalescence into large wear particles or wedges/prows, as well as competitive growth and occasional fracture of A Discrete Element Formalism for Modelling Wear Particle Formation in Contact Between Sliding... 9 such interfacial elements. It was originally developed for the simulation of granular materials, but is now widely used to simulate complex aspects of solids behavior including multiple fractures, contact interaction and friction, mass transfer, and mixing effects. In particular, this method is used to study the mesoscopic mechanisms of the formation of accretionary wedges at convergent plate boundaries in the lithosphere (an example of the contact problem on the geological scale) [32, 33]. The authors of this paper develop the original formulation of DEM, namely, the method of homogeneously deformable discrete elements [34]. This formulation makes it possible to implement models of the mechanical behavior of materials with complex rheological properties (including elastic-plastic) [35], as well as to take into account surface adhesion. Using these models, the basic modes of deformation and fracture of initial mesoscopic asperities were systematically analyzed for the first time and the influence of surface adhesion and macroscopic material properties was generalized [36, 37]. Despite large capabilities, the formalism of the method needs further development to study the evolution of transfer fragments in metallic friction pairs, including the formation of large interfacial elements and the competition of their evolutions. In particular, the model of contact interaction of discrete elements should take into account the key features of cold welding between rough (and potentially contaminated) surfaces [38-40]. The model of element-element bond-breaking should be based on a “universal” (dissipative) criterion, which can be applied both to elastic-plastic material in the initial state and to newly formed (cold-welded) fragments. This paper is devoted to the description of the advanced model of discrete element int eraction, which takes into account the features described above. We show that the model allows adequate describing the formation and evolution of large interfacial elements as well as the competition between interfacial elements in the contact zone. 2. DISCRETE ELEMENT MODEL OF CONTACT INTERACTION 2.1. The keystones of the model Within the framework of the formalism of DEM, the simulated material is represented by an ensemble of interacting elements (particles). A discrete element models a certain finite volume (region or fragment) of a material. This volume is bounded by a well- defined surface and contains a sufficient number of atoms or molecules for an integral description in terms of thermodynamic parameters. Mechanical interaction of these finite volumes (discrete elements) leads not only to their translational motion but also to the rotation, which is determined by the moments of interaction forces [41, 42]. The authors use the most widely used implementation of the discrete element method based on the approximation of equivalent disks (2D models) or balls (3D models) [43, 44]. Within the framework of this approximation, discrete elements modeling areas of the consolidated material are assumed to have an equiaxial polygonal shape and interact through flat faces. The motion of such elements is described by the classical system of Newton-Euler equations of motion for inscribed discs or balls having the same values of mass and moment of inertia as describing polyhedra. In this case, the central force of interaction of the two elements leads to their translational motion along the line connecting the centers of the approximating disks or balls, and the tangential force causes motion in the transverse plane and rotation. A pair of interacting elements can be in a 10 E.V. SHILKO, A.S. GRIGORIEV, A.YU. SMOLIN linked or unlinked state. The first case corresponds to the presence of a chemical bond between the contact surfaces of the elements, while the second case simulates the classical adhesionless contact interaction of bodies. Unlike the traditionally used models of rigid discrete elements (they apply only to brittle materials and environments), the authors consider a discrete element as deformable [45, 46]. In this work, the approximation of a uniform distribution of deformations and stresses in the volume of an element is used. This approximation is similar to the widely used approximation of linear interpolation of displacements (first-order polynomial) in the finite element method. Within this approximation, the stress-strain state of the discrete element i is characterized by the average stress i   and strain i   tensors (, = x,y,z; XYZ is the laboratory coordinate system). The components of the averaged stress tensor are calculated in terms of the forces of interaction of an element (approximated by an equivalent disk) with its neighbors: 0 0 1 [( ) ( ) ]( ) iN i i ij ij ij ij ij ij ji R S n t n V     =  =    , (1) where Ri is the radius of the equivalent disk/ball, 0 ij S is the initial value of the area of contact of the element i with the neighbor j, Vi 0 is the initial volume of the element, Ni is the number of neighbors of the element i, ijn  and ijt  are the unit vectors oriented along the normal and in the plane of contact of elements i and j, ij np ijij SF= and ij tp ijij SF= are specific values of central ( np ij F ) and tangential ( tp ij F ) forces of element-element interaction, Sij is the current value of the area of contact. The components of the average strain tensor i   are calculated by analogy. The formalism of deformable discrete elements allows modelling the mechanical behavior of materials with complex rheological properties. The specific form of relations for the forces of central and tangential interaction of elements is determined by the form of the constitutive law (elasticity, plasticity, etc.) for the simulated material and usually duplicates it. This paper deals with elastic-plastic metallic materials. The elastic response of element i to mechanical action from element j is written in a form similar to the generalized Hooke's law. In hypoelastic form (in increments), the expressions for central and tangential (shear) forces are as follows:      =          −+= ijiij i mean i i ijiij G K G G 2 3 2 12 . (2) Here, Gi and Ki are the shear modulus and bulk modulus of the material of element i, symbol  denotes the increment of the corresponding parameter per integration step, ij and ij are central (normal) strain and shear angle of element i in pair i-j (contributions of the element to the total values of normal and shear strains of the pair), i mean  is mean stress in the element i (the first invariant of the tensor i   ). The mechanical behavior of an element beyond the elastic limit (stress relaxation and accumulation of inelastic strains) is described in the framework of the associated plastic flow law using the von Mises yield criterion. Within the approximation of a homogeneously deformable element, the second invariant of the deviator of the average stress tensor is used as a yield criterion: A Discrete Element Formalism for Modelling Wear Particle Formation in Contact Between Sliding... 11 ( ) i in eq y eq  =   , (3) where y( in eq  ) is the unified hardening curve of the modelled material in terms of true stress, in eq  is the equivalent plastic strain in the element i. The plasticity of the element is numerically implemented using the radial return algorithm, developed by M.L. Wilkins. 2.2. The model of cold pressure welding of discrete elements Contact interaction and friction of ductile (in particular metallic) materials with high surface energy is typically accompanied by local effects of cold welding between the surfaces [47-52]. It is known that chemical bonding of clean and atomically smooth contacting surfaces of metals can be achieved even in the absence of significant contact pressure. However, complete welding of rough surfaces requires the application of a contact pressure comparable to or exceeding the material yield stress. In the general case of surfaces overlaid by cover layers and/or contaminants, the required contact pressure can be multiples of the yield stress [38-40]. At the same time, numerous experiments indicate that additional plastic shear deformation is capable of providing chemical bonding of surfaces of even dissimilar metallic materials under contact loads comparable to the yield stress (this is observed, in particular, when rolling bimetallic plates or twisting two pressed disks). Efficient bonding of contacting surfaces under the condition “compression+shear” is ensured, among other things, due to interpenetration and mixing of the surface layers of counterbodies, as well as dispersion of cover layer and contaminant. Thus, the contact pressure can be used as an appropriate criterion for the beginning/ development of the process of cold welding between the surfaces. Its value should exceed a certain threshold, sufficient for localization of plastic deformation on the rough surfaces and complete contact of these surfaces, as well as for dispersion of possible cover layers and contaminants. In this case, the condition for achieving complete adhesion of surfaces is the achievement of a certain threshold (minimum) value of accumulated plastic strain in the surface layers. It is clear that the threshold values of the two mentioned parameters should correlate with the material yield stress and the plastic strain before material failure. However, such a relationship is not one-valued due to a large number of other factors (roughness, contaminants, dissimilarity of contacting materials, etc.). Therefore, the values of yield stress and plastic strain to fracture (or plastic work of deformation) under tension can be considered as some reference values for characterizing the condition and criterion of cold welding between surfaces [40]. The above considerations, combined with the principles of a well-known phenomenological model by N. Bay [38-40], form the basis for a phenomenological numerical model of cold welding of contacting discrete elements, which simulate parts of interacting surfaces or third body fragments. The main theses of the model are as follows: 1) Cold welding of a pair of contacting (unlinked) discrete elements i and j proceeds at the current integration step only if the condition (–ij)>p0 weld is satisfied. Here, ij is the current specific value of the normal interaction force (it is negative when the pair is compressed), and p0 weld is the contact pressure threshold (input parameter of the model). 2) Cold welding between the surfaces of the contacting elements i and j is a progressive process and develops as the plastic straining of these elements. Accordingly, the area of the contact of the elements (Sij) generally consists of two parts: the area of the 12 E.V. SHILKO, A.S. GRIGORIEV, A.YU. SMOLIN chemically bonded (welded) part ij linkij ij link kSS = and the area of the unbonded (contact) part (1 ) ij ij link ij link S S k= − . Here, 10  ij link k is a fraction of the welded part of the contact surface. We assume linear dependence of the rate of increase in the fraction of chemically bonded contact surface on the rate of plastic straining of the pair: weld ij plastij link W W k  = , (4) where  is the increment in the value of the corresponding variable per time step of numerical integration of the equations of motion, ij plast W is the specific values of plastic work of deformation in a pair (elastic strain energy density dissipated as a result of plastic straining), Wweld is the specific value of the plastic work of deformation sufficient for complete cold welding of contacting elements. 3) In the general case, the value of Wweld is a decreasing function of the contact pressure pn (pn>p0 weld). The simplest form of dependence is linear: 0 0 1 0 1 0 ( ) ( ) weld weld weld weld n weld n weld weld p p W p W W W p p − = + − − . (5) Here, W0 weld=Wweld(p0 weld) is the plastic work sufficient for complete welding of elements at the minimum value of contact pressure at which this process is stable, and W1 weld=Wweld(p1 weld) is the plastic work sufficient for complete welding at some p1 weld > p0 weld. In the numerical implementation of the model, the value of Wweld in Eq. (4) is calculated according to Eq. (5) at each integration step with the use of |ij| as the contact pressure pn. 4) The increment in the specific work of plastic deformation in a pair i-j ( ij plast W ) is calculated taking into account the fact that the elements in a pair are deformed differently. In other words, the increment of plastic work of deformation of element i in pair i-j ( ( )i j plast W ) is not equal to the increment of plastic work of deformation of element j in this pair ( ( )j i plast W ). Therefore, ( ) ( )ij i j j i plast plast plast W W W =  + (plastic straining of both elements in a pair contribute to a tight convergence of the contacting surfaces). The increment in the specific work of plastic deformation of element i in pair i-j ( ( )i j plast W ) is calculated through the increment in its equivalent strain in this pair using a unified hardening curve ( ) in y eq   of the material of element i. The increment of the work of plastic deformation of the element j ( ( )j i plast W ) is determined in a similar way. Detailed information on the calculation of the element’s equivalent strain in a pair (in the general case, it differs from the equivalent strain calculated using element’s average strain tensor) is given in [34, 46]. The input parameters of the model can be obtained from experimental data or other theoretical constructions. These parameters are (W0 weld, p0 weld) and (W1 weld, p1 weld) for all the pairs of materials and unified hardening curves ( ) in y eq   for these materials. 2.3. The model of bond breaking in the pair of linked discrete elements Within the framework of the discrete element method, a local fracture is modeled by breaking a chemical bond between elements, that is, by the transition of a pair of elements from a linked to an unlinked state. This transition occurs out when the specified fracture criterion is fulfilled in a pair. The traditional failure criteria in the mechanics of a A Discrete Element Formalism for Modelling Wear Particle Formation in Contact Between Sliding... 13 deformable solid are “stress-based” ones. They are formulated in terms of stress tensor invariants (for example, the Mises criterion for metals and polymers). To implement such criteria within the framework of the method of deformable discrete elements, we previously developed an approach based on the determination of the local stress tensor on the contact surface of a pair of elements [34, 46, 53]. The invariants of this stress tensor are used to calculate the given “stress-based” criterion of local fracture. The bond between the elements breaks when the criterion exceeds the threshold value (bond strength of the pair). The classical “stress-based” (and similar “strain-based”) fracture criteria make it possible to adequately simulate a single fracture of the material (upon deformation from the initial unloaded state). At the same time, friction of ductile materials with high surface energy can be accompanied by multiple alternating acts of local fracture and cold welding of fragments in contact spots. The use of “stress-based” or “strain-based” fracture criteria (for example, equivalent stress or strain in the pair) for newly formed linked pairs of elements can lead to unphysical premature breakage of chemical bonds in these pairs. The reason for this is the fact that the formation of new linked pairs by the mechanism of cold welding typically occurs under high local stresses. These stresses can be close to or even exceed those required to achieve the threshold value of the fracture criterion. To solve the described fundamental problem, the authors propose to use the dissipative criterion of local fracture (breaking of a chemical bond between discrete elements) instead of “stress-based” or “strain-based”. The dissipative criterion characterizes the specific value of elastic strain energy dissipated during plastic deformation of the material. Note that in the case of conventional mechanical tests on samples of ductile material, the dissipative criterion provides the same fracture patterns as the classical Mises criterion or an analogous criterion of equivalent strain. In the developed discrete-element model, we use specific works of plastic deformation of elements i and j in the pair i-j ( ( )i j plast W and ( )j i plast W ). The method for calculating these parameters is described in the previous section when describing the numerical implementation of the cold welding model. The formulation of the dissipative fracture criterion is as follows: ( ) ( ) max{ , } i j j i plast plast fract W W W= , (6) where Wfract is the specified values (model parameter). The values of Wfract for the pair of linked elements of the same material is equal to the area under a unified loading curve ( ) in y eq   in the range from 0= in eq to = in eq s (s is inelastic equivalent strain corresponding to ultimate stress (strength) and the onset of fracture). The calculation of the parameters ( )i j plast W and ( )j i plast W is carried out in an incremental fashion (the increments of these parameters are calculated at each integration step). Note that the initial values of ( )i j plast W and ( )j i plast W are assumed to be zero for both the initially linked pairs of elements and the newly formed (cold-welded) linked pairs. Eq. (6) is a condition for the initiation (beginning) of breaking a bond in a pair of elements i and j. The process of local fracture (that is, the formation of two surfaces) is considered as distributed over time. The model of stable crack growth [34] is used as a model of such a process. Within the framework of this model, a kinetic equation is specified for the parameter ij link k (the fraction of the chemically bonded part of the contact surface): ( , ,...) ij ij ij link fract link k f k=  , where f<0 is a specified negative function and ij fract  is a 14 E.V. SHILKO, A.S. GRIGORIEV, A.YU. SMOLIN value of relative displacement of elements after the start of bond breaking. Here, we use the linear function: max,fract ij fract ij link fk −==  , (7) where fract,max is the limiting value of the parameter corresponding to a complete bond breakage between the elements ( 0=ij link k ). The parameter ij fract  increment per integration step ( ij fract  ) is calculated as a combination of the tensile and shear strain increments in the pair i-j. A rough estimate of fract,max is the length of the softening stage (from maximum resistance to zero) in the diagram of uniaxial tension of a sample of ductile material (true stress should be analyzed). 3. RESULTS OF THE ILLUSTRATIVE SIMULATION AND DISCUSSION The adequacy of the developed discrete element model and its advantages in the study of the formation and evolution of large interfacial elements (wear particles and prows/wedges) are illustrated by the example of modeling of the tangential contact of two rough surfaces of metallic materials with the same mechanical properties (Fig. 1). The simulated fragment of the contact zone includes more than 40 asperities. The modeling was carried out in a two-dimensional formulation using the approximation of plane stress or plane strain state. The linear dimensions of the model were 1000 µm along the X-axis and 220 µm along the Y-axis (the dimension in Z-axis is equal to the discrete element size 1 μm). Both surfaces have the same roughness parameters: Smi23 µm (average step of profile irregularities), Rmax5 µm (maximum profile height). Fig. 1 Structure and loading scheme of the model fragment of the zone of contact interaction of bodies with micro rough surfaces. We considered the model material with mechanical properties close to those of aluminum bronze (copper content 7%) [54, 55]. The following mechanical characteristics of the material were used: density =7800 kg/m3, Young's modulus E=118 GPa, Poisson's ratio =0.295, yield stress y=360 MPa, tensile strength UTS=560 MPa. The material is elastic-plastic with linear hardening (strain hardening modulus H=900 MPa). The value of the parameter Wfract of the dissipative fracture criterion was 141 MJ/m 3 (plastic work of deformation from y to UTS). The value of the parameter max,fract (the length of the softening stage) was assumed to be 5%. The following values of the parameters of cold welding model were used: p0 weld =y = 360 MPa, Wweld(pn) =const= A Discrete Element Formalism for Modelling Wear Particle Formation in Contact Between Sliding... 15 Wfract=141 MJ/m 3. These values of the parameters of the cold welding model can be considered as an initial ("sighting") approximation for studying the general patterns of formation of interfacial elements (third body fragments) in the contact zone. We analyzed the interaction of the counterbodies during their relative tangential displacement at a constant speed Vslide and at a given constant pressing force Fn. The force is applied to the upper surface of the upper block, its specific value is 36 MPa, or 10% of the value of the yield stress of the material. The base of the system (the lower surface of the lower block) is fixed in the Y direction and moves along the X-axis at a constant speed Vslide=5 m/s. Periodic boundary conditions are applied along the X-axis. The simulation results have shown that the type of interfacial fragments formed (and the dynamics of the change in the friction coefficient) can qualitatively change when the mechanical confinement change in the direction normal to the sliding direction (normal to the considered plane in Fig. 1). In the 2D model, two extreme cases of transverse confinement were considered: no confinement (plane stress approximation) and “rigid” confinement (plane strain approximation, out of plane strains are restricted). In the absence of transverse mechanical confinement, large and approximately equiaxial particles are formed in the contact zone. These particles are weakly bonded to the surfaces. The formation of such particles includes three main stages. At the initial stage of friction (stage I), the surfaces are “matching”. Namely, there is cutting or abrasion of the initial surface asperities and the formation of the original (small) fragments of the third body (torn off or cut off fragments of asperities). Note that even these small fragments totally separate the surfaces (Fig. 2a). Some of the fragments undergo strong compression by blocks. As a result, they are strongly hardened (in comparison with the adjacent surface layers of blocks) and adhere to the surfaces by the mechanism of cold welding. In the process of relative tangential displacement of surfaces, such fragments move by rolling. In the course of turning, they tear off layers of material from the adjacent areas of the surfaces and gradually increase in size (Figs. 2b and 2c). So, already at the initial stage, the contact spots are determined not by the current roughness of the surfaces, but by the most intensively growing fragments of the third body. Fig. 2 An enlarged fragment of the central part of the contact zone (Fig. 1) at different values of relative tangential displacement of the blocks: (a) 35 µm; (b) 80 µm. 16 E.V. SHILKO, A.S. GRIGORIEV, A.YU. SMOLIN The largest (and the most deformed and hardened) fragments of the third body grow faster than the others, since in the process of rolling they sequentially adhere to less hardened areas of both surfaces and tear off layers of material from them. The shape of these particles gradually becomes close to round. As they grow, these particles take on an increasing part of the total contact load and therefore are continuously compacted and hardened. At a certain stage (stage II), the process of their growth becomes strongly competitive. At this stage, the number of “load carrying” particles gradually decreases (Fig. 3). The end of stage II is associated with a decrease in the number of “load carrying” particles (Fig. 4) to a certain minimum level. The growth of the rest (non-bearing) particles stops, and later they coalesce with the largest ones. Fig. 3 The structure of the contact zone (stage II) at different values of relative tangential displacement of blocks: (a) 155 μm; (b) 240 μm. Fig. 4 Individual “load carrying” particles from the right part of Fig. 3b. At stage III, friction passes into a steady state. Here, the “load carrying” particles continue to grow, also due to the attachment of smaller (non-carrying) interfacial fragments (Fig. 5). Note that, the attached interfacial fragments change the shape of “load carrying” particles to significantly nonequaxial and irregular. Rotation of such an irregular “lump” A Discrete Element Formalism for Modelling Wear Particle Formation in Contact Between Sliding... 17 constrained by the surfaces is difficult. So, the “lump” stops rotation and acts as a wedge. Such wedges become the most loaded interfacial fragments and are destroyed (or lose large fragments). Episodically repeating phenomena of wedging and subsequent destruction of the largest and nonequiaxial particles control the maximum characteristic size of these fragments of the third body. Fig. 5 The structure of the contact zone (stage III) at different values of relative tangential displacement of blocks: (a) 425 μm; (b) 650 μm; (c) 870 μm. So, in the absence of transverse mechanical confinement of the contact zone, the sliding friction is actually reduced to adhesive rolling friction caused by the rotational motion of load-carrying particles and accompanied by the accretion of the material layers from surfaces on these particles. The running-in phase (stages I and II) and transition to the steady-state friction (stage III) are well manifested in the dependence of the coefficient of friction (COF) on the sliding distance (Fig. 6). At the stage of “nucleation” of the initial transfer fragments (stage I), the friction coefficient increases and can exceed 1. At the stage of competitive growth of interfacial fragments (stage II), the friction coefficient decreases to a level corresponding to the stationary regime. At the stage of transition to the steady-state friction regime (stage III), the friction coefficient varies around a certain (approximately constant) average value. The most significant deviations are associated with the formation and subsequent destruction of large “wedges” (lumps of fragments). Note that COF in Fig. 6 characterizes friction of the mesoscale area of the contact zone. It is clear that the amplitude of variations in the macroscopic COF should be many times smaller. 18 E.V. SHILKO, A.S. GRIGORIEV, A.YU. SMOLIN In the case of “rigid” mechanical confinement of the contact zone (imitated by the approximation of a plane strain state in a 2D formulation), the development of the frictional process differs from that described above not only quantitatively, but qualitatively. After the initial stage of surface matching, the peculiar surface “incrustations” of a lumpy or wedge shape are formed and gradually grow. Each of these wedges strongly adheres to one of the surfaces and moves with it. Similar to a weakly confined contact zone, the growth of wedges during friction is a competitive process. Wedges increase in size due to the attachment of smaller fragments of the “third body” (transfer fragments), as well as due to the shearing of layers of material from the opposite surface (Fig. 7). Reaching some critical size, these wedges can crush due to a high concentration of shear stresses. Therefore, similar to the case of a mechanically unconfined system, the maximum size of the wedges does not exceed a certain characteristic scale. Fig. 6 Dependence of the instantaneous COF on the relative tangential displacement of the blocks. The red line is the result of averaging the instantaneous values in a window of 100 points. Key stages I-III are designated in Roman numerals (stages are delimited by green vertical lines). The horizontal dashed line is the average value of the friction coefficient at stage III. Plane stress state approximation. Fig. 7 The structure of the contact zone at different values of relative tangential displacement of blocks: (a) 190 μm; (b) 300 μm. The frame in (b) highlights the largest and the most heavily loaded wedge. Plane strain approximation. A Discrete Element Formalism for Modelling Wear Particle Formation in Contact Between Sliding... 19 So, under the condition of “rigid” mechanical confinement of the surface layers, the sliding friction proceeds in a classical manner considered in a recent review by J.A. Greenwood [5]. Note that the friction coefficient in a mechanically confined friction pair reaches 1.2-1.4 at stage I and then decreases only insignificantly (the average value at the stage of steady-state friction is 1.0-1.1). A high value of the local coefficient of friction in mechanically confined areas of the contact zone is defined by the formation of strongly hardened wedges (instead of rolling particles), which act as an abrasive. The simulation results show that the friction can proceed in a significantly different manner in the contact zone regions characterized by different degrees of mechanical confinement of the surface layers. Particles of approximately equiaxial (and mostly round) shape are predominantly formed in low confined regions. These particles act as a kind of “rollers” and provide a relatively low value of resistance to tangential movement (low value of local COF). These interfacial fragments are apparently the wear particles that will be carried away from the contact zone due to the relatively weak adhesion to the surfaces. The formation of such particles seems to be the predominant wear mechanism in the regions of the contact zone characterized by low mechanical confinement. Wedges, the well-known interfacial elements in tribology, are predominantly formed in highly confined regions of the contact zone. These interfacial fragments have a large contact area and strong adhesion to one of the surfaces, move with it and, generally speaking, are not wear particles. The largest wedges are strongly hardened and therefore can play a role similar to that of an abrasive and produce high characteristic values of the local coefficient of friction. So, wear in highly confined regions of the contact zone is mainly due to the removal of (1) small fragments formed during crushing of the most loaded wedges and (2) small surface irregularities cut by the abrasive wedges. The relative contributions and balance of the described processes in a real three- dimensional contact zone (and, consequently, the value of the macroscopic coefficient of friction and the wear rate) should depend on the ratio of squares of low and highly confined regions. 4. CONCLUSION The proposed numerical models of the formation and breakage of chemical bonds in metallic materials are based on the use of a unified energy-based (dissipative) parameter, namely, the plastic work of deformation. The unification ensures the consistency of these models and the unity of the phenomenological description of opposing (competing) processes of fracture and cold welding under the contact interaction of ductile materials. The implementation of these models within the method of homogeneously deformable discrete elements made it possible not only to adequately model the formation and evolution of qualitatively different interfacial (third body) fragments but also to identify conditions favorable for the formation of specific types of interfacial fragments. In particular, we have shown that either equiaxial (rounded) particles weakly adhered to surfaces, or wedges/prows are predominantly formed depending on the degree of mechanical confinement of the contact region. Accordingly, the degree of mechanical confinement should largely determine the value of the friction coefficient, as well as the wear rate and the size distribution of wear particles. 20 E.V. SHILKO, A.S. GRIGORIEV, A.YU. SMOLIN The results of the preliminary study are relevant for a deeper understanding of the conditions for the implementation of various mechanisms of wear during sliding friction, as well as approaches to the control of these mechanisms. Note that the shown regularities of friction and wear of ductile materials qualitatively differ from those observed in experiments (and modeling) on the friction of brittle (ceramic) materials. The low surface energy of brittle materials provides the formation of the third body fragments (particles) by means of fracture of asperities or crushing of large third body fragments. This debris not only separates the interacting surfaces and acts as rollers and/or abrasive but fills in the surface irregularities and redistributes (equalizes) the load between them. Accordingly, surface profile matters for the tangential contact of brittle materials. This ensures the applicability of analytical and semi- analytical models of contact mechanics. As evidenced by the above-shown calculation results, such models are poorly applicable to describe the tangential contact of ductile materials, since the relief of interacting surfaces is of minor importance as it affects only the dynamics of the growth and breakage of “rollers” and wedges. This problem (including the applicability of contact mechanics models for describing the tangential contact of ductile materials) is discussed in a recent conceptual paper by the classic in this field, J.A. Greenwood [5]. The simulation results shed light on some of the important issues raised in this paper. Acknowledgement: The study was carried out with the financial support of the Russian Science Foundation (Project No. 20-19-00743). REFERENCES 1. Rabinowicz, E., 1995, Friction and Wear of Materials, Wiley-Interscience, New York, 336 p. 2. Holm, R., 2013, Electric Contacts: Theory and Application, Springer Science & Business Media, Berlin, 484 p. 3. Wen, S., Huang, P., 2017, Principles of Tribology, Wiley, New York, 608 p. 4. Schirmeisen, A., 2013, One atom after the other, Nature Nanotechnology, 8, pp. 81-82. 5. Greenwood J.A., 2020, Metal transfer and wear, Frontiers in Mechanical Engineering, 6, article 62. 6. Godet, M., 1984, The third-body approach: a mechanical view of wear, Wear, 100, pp. 437-452. 7. Bethier, Y, 1996, Maurice Godet’s third body, The Third Body Concept: Interpretation of Tribological Fenomena, C.M. Taylor, T.H.C. Childs, Y. Berthier, L. Flamand, G. Dalmaz, D. Dowson, A. Lubrecht, J.M. Georges (Eds.), Elsevier, Amserdam, pp. 21-30. 8. Denape, J., 2014, Third body concept and wear particle behavior in dry friction sliding conditions, Tribological aspects in modern aircraft industry, M. Karama, K. Delbe, J. Denape (Eds.), Trans Tech Publications, Stafa-Zurich, pp. 1-12. 9. Hokkirigawa, K., 1991, Wear mode map of ceramics, Wear, 151, pp. 219-228. 10. Savchenko, N.L., Filippov, A.V., Tarasov, S.Yu., Dmitriev, A.I., Shilko, E.V., Grigoriev, A.S., 2018, Acoustic emission characterization of sliding wear under condition of direct and inverse transformations in low temperature degradation (LTD) aged Y-TZP AND Y-TZP-AL2O3, Friction, 6, pp. 323-340. 11. de Rooij, M.B., Schipper, D.J., 2001, Analysis of material transfer from a soft workpiece to a hard tool: Part II – experimental verification, ASME Journal of Tribology, 123, pp. 474-478. 12. de Rooij, M.B., van der Linde, G., Schipper, D.J., 2013, Modelling material transfer on a single asperity scale, Wear, 307, pp. 198-208. 13. Tarasov, S.Y., Filippov, A.V., Kolubaev, E.A., Kalashnikova, T.A., 2017, Adhesion transfer in sliding a steel ball against an aluminum alloy, Tribology International, 115, pp. 191-198. 14. Hokkirigawa, K., Kato, K., 1988, An experimental and theoretical investigation of ploughing, cutting and wedge formation during abrasive wear, Tribology International, 21(1), pp. 51-57. 15. Mahato, A., Guo, Y., Sundaram, N.K., Chandrasekar, S., 2014, Surface folding in metals: a mechanism for delamination wear in sliding, Proceedings of the Royal Society A, 470, article 20140297. 16. Iordanoff, I., Berthier, Y., Descartes, S., Heshmat, H., 2002, A review of recent approaches for modeling solid third bodies, ASME Journal of Tribology, 124, pp. 725-735. A Discrete Element Formalism for Modelling Wear Particle Formation in Contact Between Sliding... 21 17. Vakis, A.I., Yastrebov, V.A., Scheibert, K., Nicola, L., Dini, D., Minfray, C., Almqvist, A., Paggi, M., Lee, S., Limbert, G., Molinari, J.-F., Anciaux, G., Aghababaei, R., Echeverri Restrepo, S., Papangelo, A., Cammarata, A., Nicolini, P., Putignano, C., Carbone, G., Stupkiewicz, S., Lengiewicz, J., Costagliola, G., Bosia, F., Guarino, R., Pugno, N.M., Müser, M.H., Ciavarella, M., 2018, Modeling and simulation in tribology across scales: An overview, Tribology International, 125, pp. 169-199. 18. Li, Q., Pohrt, R., Lyashenko. I.A., Popov, V.L., 2020, Boundary element method for nonadhesive and adhesive contacts of a coated elastic half-space, Proceedings of the Institution of Mechanical Engineers Part J: Journal of Engineering Tribology, 234(1), pp. 73-83. 19. Li, Q., Popov. V.L., 2018, Boundary element method for normal non-adhesive and adhesive contacts of power-law graded elastic materials, Computational Mechanics, 61, pp. 319-329. 20. Bazrafshan, M., de Rooij, M.B., Schipper, D.J., 2018, On the role of adhesion and roughness in stick-slip transition at the contact of two bodies: A numerical study, Tribology International, 121, pp. 381-388. 21. Woldman, M., van der Heide, E., Tinga, T., Masen, M.A., 2017, A finite element approach to modeling abrasive wear modes, Tribology Transactions, 60(4), pp. 711-718. 22. Bochkareva, S.A., Panin, S.V., Lyukshin, B.A., Lyukshin, P.A., Grishaeva, N.Yu., Matolygina, N.Yu., Aleksenko, V.O., 2020, Simulation of frictional wear with account of temperature for polymer composites, Physical Mesomechanics, 23, pp. 147-159. 23. Aghababaei, R., Warner, D.H., Molinari, J.-F., 2017, Critical length scale controls adhesive wear mechanisms, Nature Communications, 7, article 11816. 24. Molinari, J.-F., Aghababaei, R., Brink, T., Frerot, L., Milanese, E., 2018, Adhesive wear mechanisms uncovered by atomistic simulations, Friction, 6, pp. 245-259. 25. von Lautz, J., Pastewka, L., Gumbsch, P., Moseler, M., 2016, Molecular dynamic simulation of collision- induced third-body formation in hydrogen-free diamond-like carbon asperities, Tribology Letters, 63, article 26. 26. Vargonen, M., Yang, Y., Huang, L., Shi, Y., 2013, Molecular simulation of tip wear in a single asperity sliding contact, Wear, 307, pp. 150-154. 27. Yang, Y., Huang, L., Shi, Y., 2016, Adhesion suppresses atomic wear in single-asperity sliding, Wear, 352-353, pp. 31-41. 28. Dmitriev, A.I., Nikonov, A.Yu., Osterle, W., Jim, B.C., 2019, Verification of Rabinowicz’ criterion by direct molecular dynamics modeling, Facta Univesitatis-Series Mechanical Engineering, 17(2), pp. 207-215. 29. Milanese, E., Molinari, J.-F., 2020, A mechanistic model for the growth of cylindrical debris particles in the presence of adhesion, International Journal of Solids and Structures, 203, pp. 1-16. 30. Mishra, T., de Rooij, M., Shisode, M., Hazrati, J., Schipper, D.J., 2020, A material point method based ploughing model to study the effect of asperity geometry on the ploughing behaviour of an elliptical asperity, Tribology International, 142, article 106017. 31. Mamalis, A.G., Vortselas, A.K., Panagopoulos, C.N., 2013, Analytical and numerical wear modeling of metallic interfaces: A statistical asperity approach, Tribology Transactions, 56(1), pp. 121-129. 32. Burbidge, D.R., Braun, J., 2002, Numerical models of the evolution of accretionary wedges and fold-and- thrust belts using the distinct-element method, Geophysical Journal International, 148, pp. 542-561. 33. Miyakawa, A., Kinoshita, M., Hamada, Y., Otsubo, M., 2019, Thermal maturity structures in an accretionary wedge by a numerical simulation, Progress in Earth and Planetary Science, 6, article 8. 34. Shilko, E.V., Psakhie, S.G., Schmauder, S., Popov, V.L., Astafurov, S.V., Smolin, A.Yu., 2015, Overcoming the limitations of distinct element method for multiscale modeling of materials with multimodal internal structure, Computational Materials Science, 102, pp. 267-285. 35. Shilko, E.V., Smolin, A.Yu., Dimaki, A.V., Eremina, G.M., 2021, Particle-based approach for simulation of nonlinear material behavior in contact zones, in: G.-P. Ostermeyer, V.L. Popov, E.V. Shilko, O.S. Vasiljeva (Eds.), Multiscale Biomechanics and Tribology of Organic and Inorganic Systems, Springer, Berlin, pp. 67-89. 36. Dimaki, A.V., Shilko, E.V., Dudkin, I.V., Psakhie, S.G., Popov, V.L., 2020, Role of adhesion stress in controlling transition between plastic, grinding and breakaway regimes of adhesive wear, Scientific Reports, 10, article 1585. 37. Dimaki, A.V., Dudkin, I.V., Popov, V.L., Shilko, E.V., 2019, Influence of adhesion force and strain hardening coefficient of the material on the rate of adhesive wear in a dry tangential frictional contact, Russian Physics Journal, 62 (8), pp. 1398-1408. 38. Bay, N., 1979, Cold pressure welding – the mechanisms governing bonding, Journal of Engineering for Industry, 101(2), pp. 121-127. 39. Bay, N., 1983, Mechanisms producing metallic bonds in cold welding, Welding Journal, 62, pp. 137-142. 40. Zhang, W., Bay, N., 1996, A numerical model for cold welding of metals, Annals of the CIRP, 45(1), pp. 215-220. https://www.nature.com/articles/ncomms11816 https://www.nature.com/articles/ncomms11816 22 E.V. SHILKO, A.S. GRIGORIEV, A.YU. SMOLIN 41. Jing, L., Stephansson, O., 2007, Fundamentals of discrete element method for rock engineering: theory and applications, Elsevier, Amsterdam, 562 p. 42. Bicanic, N., 2017, Discrete element methods, Encyclopaedia of Computational Mechanics, Stein E, de Borst R, Hughes TJR (Eds.), Wiley, New York, pp. 1-38. 43. Potyondy, D.O., Cundall, P.A., 2004, A bonded-particle model for rock, International Journal of Rock Mechanics and Mining Sciences, 41, pp. 1329-1364. 44. Ivars, D.M., Pierce, M.E., Darcel, C., Reyes-Montes, J., Potyondy, D.O., Young, R.P., Cindall, P.A., 2011, The synthetic rock mass approach for jointed rock mass modelling, International Journal of Rock Mechanics and Mining Sciences, 48, pp. 219-244. 45. Psakhie, S., Shilko, E., Smolin, A., Astafurov, S., Ovcharenko V., 2013, Development of a formalism of movable cellular automaton method for numerical modeling of fracture of heterogeneous elastic-plastic materials, Frattura ed Integrita Strutturale, 24, pp. 26-59. 46. Psakhie, S.G., Shilko, E.V., Grigoriev, A.S., Astafurov, S.V., Dimaki, A.V., Smolin, A.Yu., 2014, A mathematical model of particle–particle interaction for discrete element based modeling of deformation and fracture of heterogeneous elastic–plastic materials, Engineering Fracture Mechanics, 130, pp. 96-115. 47. Budakian, R., Putterman, S.J., 2002, Time scales for cold welding and the origins of stick-slip friction, Physical Review B., 65, article 235429. 48. Alcantar, N.A., Park, C., Pan, J.M., Israelachvili, J.N., 2003, Adhesion and coalescence of ductile metal surfaces and nanoparticles, Acta Materialia, 51, pp. 31-47. 49. Cha, S.-H., Park, Y., Han, J.W., Kim, K., Kim, H.-S., Jang, H.-L., Cho, S., 2016, Cold welding of gold nanoparticles on mica substrate: Self-adjustment and enhanced diffusion, Scientific Reports, 6, article 32951. 50. Wagle, D.V., Baker, G.A., 2015, Cold welding: a phenomenon for spontaneous selfhealing and shape genesis at the nanoscale, Materials Horizons, 2, pp. 157-167. 51. Singh, R., Gupta, P., Yedla, N., 2019, Single-crystal Al–Cu50Zr50 metallic glass cold welds: tensile and creep behavior, Molecular Simulation, 45, 1549-1562. 52. Singh, R., Sharma, V., 2020, Molecular dynamics study of tensile behaviour for cold and linear friction welded single crystal tungsten, Journal of Molecular Graphics and Modelling, 99, article 107655. 53. Dimaki, A.V., Shilko, E.V., Popov, V.L., Psakhie, S.G., 2018, Simulation of fracture using a mesh- dependent fracture criterion in a discrete element method, Facta Univesitatis-Series Mechanical Engineering, 16(1), pp. 41-50. 54. Yin, Z., Sun, L., Yang, J., Gong, Y., Zhu, X., 2016, Mechanical behavior and deformation kinetics of gradient structured Cu-Al alloys with varying stacking fault energy, Journal of Alloys and Compounds, 687, pp. 152-160. 55. Huang, C.X., Hu, W., Yang, G., Zhang, Z. F., Wu, S.D., Wang, Q.Y., Gottstein, G., 2012, The effect of stacking fault energy on equilibrium grain size and tensile properties of nanostructured copper and copper–aluminum alloys processed by equal channel angular pressing, Materials Science and Engineering A, 556, pp. 638-647. https://www.sciencedirect.com/science/article/pii/S1365160910002169