# Adaptive FEM-DEM bridging coupling for third-body evolution in sheared granular media

Manon Voisin-Leprince, Jean-François Molinari, Guillaume Anciaux

► To cite this version:

Manon Voisin-Leprince, Jean-François Molinari, Guillaume Anciaux. Adaptive FEM-DEM bridging coupling for third-body evolution in sheared granular media. 2025. hal-05014905v2

HAL Id: hal-05014905

https://hal.science/hal-05014905v2

Preprint submitted on 19 Aug 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Adaptive FEM-DEM bridging coupling for third-body evolution in sheared granular media

Manon Voisin-Leprince¹, Jean-François Molinari¹, and Guillaume Anciaux¹

¹ Institute of Civil Engineering, Institute of Materials Science and Engineering, École polytechnique fédérale de Lausanne (EPFL), CH 1015 Lausanne, Switzerland

The frictional behavior of dry solids in sliding contact is significantly affected by the formation and accumulation of wear particles at the interface. These particles create an amorphous layer known as the third body, which undergoes intense deformation. We show that accurately capturing the evolution of the third body requires a precise modeling of the physical dimensions of the surrounding elastic bodies and boundary conditions. This implies that a computational tribology simulation engine must resolve both the fine-scale dynamics of wear particles and the large-scale structural dimensions of the sliding bodies. To address this challenge, we employ a coupled Finite Element Method (FEM) and Discrete Element Method (DEM) approach. The DEM is used to model the third body, where grain-scale discretization is crucial, while the FEM represents the first bodies as a continuum, reducing computational cost. The coupling between FEM and DEM is achieved using the bridging method. To allow unrestricted growth of the third body beyond the initial discrete region, we implement an adaptive coupling mechanism that enables FEM regions to transition into DEM regions. This transition occurs when a predefined criterion is met, based on the average change in neighboring particles. Our approach is validated for both amorphous materials and ordered crystalline lattices. Finally, as a proof of concept, we apply the adaptive coupling framework to simulate the evolution of a third body initially composed of stiff elliptical particles.

Keywords Computational tribology, discrete-element method, finite-element method, adaptive coupling, granular flow

## 1 Introduction

In tribology, a third body (or gouge) refers to a layer of materials located in between two contacting surfaces. Friction and wear are intimately linked to the evolution of the third body, which is influenced by various factors. In other fields, a third body may correspond to material features, that always involve surfaces and disorder. For instance in metals and alloys a third body is a surface layer where grain properties are changing, sometimes going amorphous. Numerical simulations provide useful insight in material science (Tadmor and Miller 2014). More specific to the study of gouges, Mollon (2021) employed a deformable discrete element particle method to investigate the impact of grain properties on the third-body flow regime, friction, and energy dissipation (Mollon 2019). It was demonstrated that the third body modifies the macroscopic friction through its ability to change the mesostructures and force chains. Y. Zhang et al. (2020) explored the influence of cohesion between particles on the third-body rheology. These simulations (Mollon 2021; Y. Zhang et al. 2020) revealed complex rearrangements, with particles detaching, re-attaching, and porosity changes.

Pham-Ba et al. (2022) developed a coarse-grained contact model that determines the inter-particle force parameters based on the desired macroscopic material properties. This model enables the replication of key wear mechanisms observed in atomic scale simulations (Aghababaei et al. 2016), while significantly reducing computational cost. The model was used (Pham-Ba et al. 2023) to understand the influence of macroscopic material properties (tensile strength, adhesion) on the emergence of different third-body configurations, including the formation of detached wear particles, of a shear band, or of an intermediate regime, as observed in Mollon (2019); Y. Zhang et al. (2020).

However, all these previous studies were conducted on relatively small Discrete Element Method (DEM) systems, with regions surrounding the third body, i.e. the first bodies, modeled
with a thin layer of particles to provide boundary conditions. A later attempt by Bouillanne et al. (2024) included the surrounding bodies modeled as a linear elastic medium, but a systematic investigation of how far-field boundary conditions affect the evolution of the third body remains to be carried out.

The objective of the manuscript is two-fold. First, we aim to show that an accurate representation of boundary conditions is crucial to capture the correct deformation mechanisms within the third body. Second, as the third-body thickness may grow with time, we propose a novel adaptive coupling scheme between a granular domain and surrounding elastic continuum media.

The manuscript therefore explores strategies to reduce the heavy computational cost of discrete element simulations. We use DEM to model the accumulation of wear particles in the third body, as here a grain discretization is essential, while we employ the Finite-Element Method (FEM) to efficiently represent the surrounding regions as continuum domains.

Coupling methods have a long history, in particular for coupling crystalline Molecular Dynamics (MD) with FEM. Generally, two main classes of coupling can be distinguished: edge-to-edge coupling (Rudd et al. 1998; Broughton et al. 1999; Tadmor, Ortiz, et al. 1996), and overlapping coupling Wagner et al. (2003); Xiao et al. (2004) for which the discrete and continuum domains coexist within part of the simulation domain. These approaches have been adapted, with important challenges, to amorphous granular media modeled using DEM, either with overlap regions (Voisin-Leprince et al. 2022; Chaudry et al. 2022; Wellmann et al. 2012), also called bridging domain Xiao et al. (2004); Dhia et al. (2005), or edge-to-edge approaches (Gao et al. 2018; Mollon 2018; Avci et al. 2012). In this paper, we will employ the FEM-DEM bridging method Voisin-Leprince et al. (2022). For these approaches, it is critical to match accurately the material properties between FEM and DEM (Wellmann et al. 2012). In an earlier work, we have determined the characteristic mesh size to employ in the FEM region Voisin-Leprince et al. (2022). We have extensively validated the approach for wave propagation, crack propagation, and wear debris formation (Voisin-Leprince et al. 2024), but in all cases maintained the partitioning between FEM and DEM domains fixed. Here, we extend the approach to adaptive domains, for which the discrete region could invade the continuum region. Chen et al. (2021) employed an adaptive overlapping coupling between a continuum domain modeled with the *Material Point Method* (MPM) and a discrete domain modeled with DEM, starting from only a continuum region at the onset of the simulations. Measuring high strain rate gradients allows the detection of shear bands formation, which was the criterion to substitute part of continuum domain for a set of re-equilibrated discrete particles. Similar strategies emerged earlier, but for "simpler" crystalline materials, as for instance in the Quasicontinuum Method (Miller et al. 2002).

In the present work, the adaptation criterion will necessitate to estimate the deformation from discrete particles. We will favor simple and efficient methods, such as the one presented in Falk et al. (1998), to adapt the coupling region, and refrain from using more costly alternatives such as Voronoi tessellation and/or Delaunay triangulation B. Zhang et al. (2015); Durán et al. (2010); Durán et al. (2010); Nguyen et al. (2020).

The outline of the paper is as follows. Section 2 presents the critical influence of the dimensions of the system on the evolution of the third body using pure DEM simulations. Then, we will focus on a system made of two FEM blocks surrounding a central DEM/MD body coupled as in Voisin-Leprince et al. (2022). Section 3 introduces the adaptive coupling method, in particular the detection strategy and adaptation criterion, while Section 4 presents validation tests for both the packing of crystalline and amorphous spheres. Finally, Section 5 examines how stiff elliptical bodies placed at the contact interface can influence the thickness of the third body and the evolution of friction.

## 2 Influence of regions surrounding a third-body

To assess how system dimensions affect third-body evolution, pure DEM systems of varying sizes are considered. The coarse-grained contact model developed by Pham-Ba et al. (2022) is used to obtain the interaction forces. The details of this model are provided in Appendix 2. The material parameters (see Table 1) are determined such that a particle assembly behaves as SiO$_2$, as was done in Pham-Ba et al. (2022). A log-normal distribution of particle sizes is selected, with a
<table><caption>**Table 1:** Material properties and particle size parameters.</caption><tr><td>Young's modulus, <i>E</i></td><td>73 GPa</td></tr><tr><td>Grain density, <i>ρ</i></td><td>2200 kg/m<sup>3</sup></td></tr><tr><td>Poisson's ratio, <i>ν</i></td><td>0.25</td></tr><tr><td>Critical length scale <i>d</i><sup>*</sup></td><td>0.1<i>L</i><sub>x</sub> nm</td></tr><tr><td>Restitution coefficient, <i>η</i></td><td>0.9</td></tr><tr><td>Mean grain diameter, <i>d</i><sub>0</sub></td><td>1.3 nm</td></tr><tr><td>Surface energy, <i>γ</i></td><td><i>d</i><sup>*</sup><i>σ</i><sub>n</sub><sup>2</sup>/(32<i>E</i>) N/m</td></tr><tr><td>Mean grain diameter, <i>d</i><sub>0</sub></td><td>1.3 nm</td></tr><tr><td>Maximum diameter, <i>d</i><sub>max</sub></td><td>1.2<i>d</i><sub>0</sub> nm</td></tr><tr><td>Minimum diameter, <i>d</i><sub>min</sub></td><td>0.8<i>d</i><sub>0</sub> nm</td></tr><tr><td>Variance, <i>σ</i><sub>p</sub></td><td>0.2(<i>d</i><sub>max</sub> − <i>d</i><sub>min</sub>) nm</td></tr></table>

mean grain diameter of $d_0 = 1.3$nm, a maximum grain diameter of $d_{max} = 1.2d_0$, a minimum grain diameter of $d_{min} = 0.8d_0$, and a variance of $0.2(d_{max} - d_{min})$nm. The system will be subject to hybrid boundary conditions, i.e. an imposed normal pressure and an imposed horizontal velocity (see Figure 1). In order to showcase the influence of boundary conditions, the shear stress $\sigma_{xy}$ will be measured at the top of the simulation cell, and will reflect changes in flow regime occurring in the third body.

For a given box size of height $50d_0$, a constant pressure of 100MPa is maintained on the top of the DEM while shearing at a steady velocity $v = 0.01c$ with c the wave speed. By adjusting the tensile strength parameter, we modeled three distinct third-body regimes as observed by Mollon (2019) and Pham-Ba et al. (2023): a mixed regime at a tensile strength of 0.1E (Figure 2 - top), a debris-rolling regime at a tensile strength of 0.11E (Figure 2 - middle), and a shear-band regime at a tensile strength of 0.2E (Figure 2 - bottom). Thus, varying the tensile strength (with constant boundary conditions) results in different flow regimes.

Then, for a given tensile strength of 0.1E, the applied top pressure was varied, along with simulation box size, to assess its impact on gouge evolution. Figure 3 shows the effective shear stress $\sigma_{xy}$ for three distinct pressures: 100 MPa, 200 MPa, and 500 MPa. At 100 MPa,

each system height exhibits different shear stress evolution and it is difficult to identify a clear trend. At pressures of 200 MPa, noticeable correlations emerge. The larger box sizes ($500d_0$, $600d_0$) produce higher shear stresses (~ 0.6 GPa) compared to smaller DEM heights ($50d_0$, $100d_0$). Rolling debris regimes are observed for the smallest DEM heights, while a regime closest to a mixed one appears for the largest ones. At 500 MPa, the flow regime turns to a shear band mechanisms with at times oscillations in the shear resistance indicating weakening mechanisms taking place within the shear band. All of these simulations reveal the richness and complexity of the rheology observed in a theoretically very simple model system. The rheology is clearly influenced by both boundary conditions (here pressure) and system size.

A simulation domain should therefore be selected to represent as closely as possible the actual system size. Given that the number of particles scales with the system size, a FEM-DEM coupling has the potential to substantially decrease the computational costs. In such a case, only the third body is modeled with DEM particles, while the first bodies are modeled with continuum

<figure><img src="image_2.png" /><figcaption>**Figure 1:** Schematic of domain geometry and hybrid boundary conditions, i.e. an imposed normal pressure and an imposed horizontal velocity, used to study shearing of an amorphous assembly of particles behaving as SiO$_2$.</figcaption></figure>
<figure><img src="image_3.png" /><figcaption>**Figure 2:** Initial position of the particles in the y-direction for a DEM box of height $50d_0$ at an external pressure of 100MPa. (a) shows a mixed regime at a tensile strength of $0.1E$, (b) depicts a wear debris regime at a tensile strength of $0.11E$, and (c) illustrates a shear band regime at a tensile strength of $0.2E$.</figcaption></figure>

<figure><img src="image_4.png" /><figcaption>**Figure 3:** Shear stress $\sigma_{xy}$ in the third-body layer (TBL) for a tensile strength of $0.1E$ at varying external pressures: (a) shows $\sigma_{xy}$ for an external pressure of 100 MPa, (b) for 200 MPa, and (c) for 500 MPa. The green dashed rectangle shows a debris rolling regime, the red dashed rectangle indicates a mixed regime, and the black dashed rectangle denotes a shear band regime.</figcaption></figure>

mechanics (linear elastic bodies) and provide accurate boundary conditions.

## 3 Adaptive coupling

In Figure 4, the discrete domain is depicted with circular particles, while the continuum domain is made visible as a black grid mesh. In total, three regions are necessary to achieve adaptive coupling: the usual bridging and boundary regions are complemented with a detection region where the criterion for the discrete domain's evolution is evaluated. The bridging region is shown in red, the boundary region in purple, and the detection region in blue. In this paper, we use simple linear-elastic, small-deformation, first-order finite elements (T3). This is considered sufficient, as no large rotations are observed in any of the cases presented in the sequel.

The coupling is realized considering a strong formulation Voisin-Leprince et al. (2022).
An appropriate criterion triggering an extension of the discrete region should capture when inelastic rearrangements of the particles occur near the bridging region. This is to ensure that strong deformations remain distant from the bridging region, which is important as a proper coupling is guaranteed only under assumptions of regularity of the strain/displacement fields and of linear elastic behavior.

Thus, to quantify the field of inelastic deformation, we evaluate the average change of neighbors per particle. For each particle in the detection region, the list of nearest neighbors is identified within a given radius $R_{mean}$. Neighbor changes are detected by comparing the current list of neighbors with a reference list. The reference list, first established at the start of the simulation, is updated after each adaptation. Then, for each particle i, the spatial average of the number of neighbors changes is calculated and denoted by $n_{i,mean}$ (refer to Figure 5).

<figure><img src="image_5.png" /><figcaption>**Figure 4:** Schematic of the adaptive coupling regions. The discrete domain is represented by circular particles, while the continuum domain is depicted with a black mesh. The boundary region of the discrete domain is marked in purple, the bridging region in red, and the detection region in blue.</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>**Figure 5:** Schematic of the spatial average of the change of neighbors over a radius $R_{mean}$.</figcaption></figure>

We denote $N_{total}$, the total number of particles within the detection region. The percentage of particles exhibiting a number of neighbor changes greater than or equal to a prescribed value, i.e. for which $n_{i,mean} \ge \alpha$, is calculated. If this percentage exceeds the threshold, $\epsilon$, the criterion for expanding the discrete domain is met:

$$ \epsilon \le \left( \frac{\sum_{i=1}^{N_{total}} H(n_{i,mean} - \alpha)}{N_{total}} \right) \times 100. \quad (1) $$

Here $H$ is the Heaviside step function. The selection of $\epsilon$ depends on the scale of the detection region, as it is influenced by $N_{total}$. For future studies on larger systems, it could be advantageous to assess neighbor changes within subdomains smaller than the detection region currently used. These subdomains could, for instance, be on the scale of a finite element.

Once decided, the discrete domain needs to be expanded. The size of the expansion is equal to the size of the current bridging region, which for the setup used in this paper turns out to also be the characteristic height of a finite element. To expand the discrete region, several steps must be followed, as illustrated in Figure 6.

The first step consists in shifting the boundary, bridging, and detection regions by a distance equal to the original height of a FE element. Consequently, the second step introduces new particles in the empty region (Step 2 in Figure 6). These new particles are taken from a periodic (in the x, y, and z directions) DEM/MD sample, having the same size distribution and material properties as those in the original and pristine discrete domain. They are placed in a geometry matching the undeformed bridging region.

In a third step, the coupling is re-initialized with these updated geometries, before imposing the FEM displacement field onto the new particles, allowing to transfer the continuum deformations (Step 3 in Figure 6). Next, the displacement of FE elements in the area transitioning from the bridging region to a purely discrete region is adjusted. The FEM nodes have to be constrained to
<figure><img src="image_7.png" /><figcaption>**Figure 6:** Schematic of the adaptive coupling process. The discrete region is depicted with circular particles, while the continuum domain is shown as a black mesh. The bridging region, where the coupling occurs, is shown in red; the boundary region, where the particles are slaved to the continuum domain, is in purple; and finally, the detection region is depicted in blue. Step 0 illustrates the system prior to adaptation. Step 1 shows the shifting of geometries, Step 2 the addition of new particles, and Step 3 the system's final configuration.</figcaption></figure>

the discrete domain's displacements, which is done by solving a least-square problem as described by Equation 7 in the appendix. These elements will deform consistently with the contained particles and will therefore create boundary conditions, for the elements in the bridging region, through the produced stresses. Additionally, the elements are visually disabled.

Finally, a relaxation process is carried out using velocity damping to stabilize the particles that are now in the bridging, and detection regions. Throughout this relaxation, all the other particles are fixed, i.e. those not in the new detection and bridging regions. To capture the system's dynamics, the velocity field of the FEM is then projected onto the bridging and detection regions to become the discrete velocity of the concerned particles. The adaptive coupling approach has been implemented in the open-source LibMultiScale software Anciaux 2024, and visualized with ParaView software Ahrens et al. (2005).

# 4 Validation

## 4.1 Adaptive FEM-MD coupling

**Geometry** The adaptive coupling is first tested with a three-dimensional MD-FEM setup. A face-centered cubic aluminum crystal is assembled with primitive lattice vectors given by $a_x = \sqrt{2}a, a_y = \sqrt{3}a, a_z = \sqrt{6}a$, where $a = 4.045 Å$ is the lattice length characteristic of ground state aluminum. Each atom within the crystal has a mass of 26.981 g/mol. The used inter-atomic potential is the classical generalized EAM potential from Finnis and Sinclair Finnis et al. (1984). The MD region is characterized by a length of $L_x = 20a_x$, a height of $L_{y,MD} = 19a_y$, and a thickness of $L_z = 2a_z$. The top and bottom FEM regions have similar dimensions in $L_x$ and $L_z$ and a height $L_{y,FEM} = 20a_y$, with 8 elements in x direction, 10 elements in y direction, and 2 elements in z direction.

The bridging (red region in Figure 7) and boundary MD (dark blue region in Figure 7) regions are chosen with a height $2a_y$ (the height of two elements), and the height of the detection region (light blue region in Figure 7) is set to $4a_y$. A void is created at the center of the MD region by removing a square region of atoms of area $a_x \times a_y/2$. Similarly to previous (amorphous) simulations, compressive forces are then applied to the FEM top surface, while the bottom surface is maintained fixed. Figure 7 shows a schematic of the simulated system. Two simulation types are compared, one non-adaptive and one adaptive both with an initial height $L_{y,MD} = 19a_y$.
For the neighbor change detection, the reference list is first assembled at the beginning of the simulation and is then updated after every adaptation, which changes the local neighborhood nearby the pad region. The choice of neighboring atoms can be computed based on distances, thanks to the regular crystalline structure.

Each neighbor list includes the **12 closest atoms** within a radius $R_2 = 2R$, with $R$ being the classical equilibrium inter-atomic distance $R = a\sqrt{2}/2$. This **accounts for** thermal vibrations that should not be assimilated as **neighbor changes**. During the simulation the current neighbors within $R_2$ distance are compared with the reference neighbor list. **If the percentage of atoms with more than three neighbor changes ($\alpha = 3$) is greater than 1%**, the evolution criterion is reached, which will trigger adaptation and a moving of the bridging domain. **Neighbor changes are measured every 100-timesteps.**

**Results** A straightforward detection of isolated dislocations within a crystal is achieved by measuring per-atom potential-energy peeks. In Figure 8, the potential energy is shown for the non-adaptive case (upper row) and for the adaptive coupling case (lower row). The three columns correspond to growing time steps, from left to right, allowing to spot positions of dislocations and to view how they interact with the bridging region. In the non-adaptive coupling case,

the dislocations soon get stuck within the bridging region. However, the adaption of the coupling region lets dislocations propagate further, as shown in Section 4.1, where the evolution of the vertical position of dislocations with time is plotted.

Removing such boundary artifacts was possible thanks to the chosen detection criterion, which triggered adaptation at the right moment. It revealed that the average change of neighbors allows to accurately capture the trace of dislocations, as illustrated in Section 4.1 where the neighbor changes correspond to the dislocations travel paths. Therefore, the presented criterion also allows atomic-precision identification of dislocation traces.

<figure><img src="image_8.png" /><figcaption>**Figure 7:** Schematic of the MD-FEM validation setup. The MD region is surrounded by two FEM regions, located at the top and bottom. The detection regions, the bridging regions, and the boundary regions are placed in the MD region.</figcaption></figure>
<figure><img src="image_9.png" /><figcaption>**Figure 8:** Potential energy in the MD, highlighting dislocation propagation. At the top, the non-adaptive coupling is represented, while at the bottom, the adaptive coupling is represented. The three columns correspond to distinct simulation times.</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>**Figure 9:** (a) Dislocation position in the y direction for adaptive and non-adaptive coupling. (b) Measure of neighbor changes for a non-adaptive MD-FEM coupling.</figcaption></figure>
## 4.2 Adaptive FEM-DEM coupling

**Geometry** The adaptive coupling is now further studied with a three-dimensional DEM-FEM setup. To model the discrete domain, the coarse-grained contact law developed by Pham-Ba et al. (2022) is employed (more details in the appendix Section 7.2). This contact law provides the spring-stiffness between particles that reproduces the desired macroscopic material properties. It was chosen so that FEM properties and the DEM macroscopic properties are matching. Similarly to Pham-Ba et al. (2022), the material properties of SiO₂ are considered, along with the particle size distribution described in Section 2.

We consider a system characterized by a length $L_x = 100d_0$, a total height $L_y = 200d_0$, and a thickness $L_z = 3d_0$. Periodic boundary conditions are chosen in the x and z directions. The FEM element characteristic sizes are selected to be $10d_0$ in the x and y directions. The thickness in the z direction is spanned by two elements. As mentioned previously, the bridging height (red region in Figure 10) is chosen as one FEM element's height, corresponding to $10d_0$. The boundary

region (dark blue region in Figure 10) height is set to $3d_0$. And the detection region height (light blue region in Figure 10) is chosen to be $7d_0$, which is enough to prevent spurious boundary effects, c.f. Voisin-Leprince et al. (2022). A constant normal pressure of 100MPa is applied to the top boundary, as well as a constant shearing velocity of $v = 0.1c$, with $c$ being the push-wave celerity. The bottom surface boundary is fixed. An initial junction of size of $30d_0$ is created in the central DEM region to facilitate the creation of a debris. Figure 10 depicts the geometry of the considered system. The monitoring of neighbor changes is made on the six closest particles within a radius $R = 2d_0$. Then, for each particle the average number of neighbor changes is computed amongst the neighbor particles within a radius $R_{mean} = 2d_0$ (see Section 3 and Figure 5). If more than 2% of the particles have an average change in neighbors greater than or equal to $\alpha = 3$, then the criterion is met. This criterion is assessed every 500-timesteps.

Three cases will be considered: 1) a non-adaptive coupling with a small DEM height $L_{z,DEM} = 50d_0$, 2) an adaptive coupling with the same initial DEM height $L_{z,DEM} = 50d_0$, and 3) an initially larger DEM region, with height $L_{z,DEM} = 70d_0$, also with an adaptive coupling. These three cases highlight the benefits provided by the coupling adaptivity.

**Results** In Figure 11, the averaged change of neighbors, for each particle, is shown within the detection regions. A significant change of neighbors is observed at the crack tip, where the cracks propagate to form the debris. This confirms that the average change of neighbors is a reliable measure of DEM inelastic deformation, which will be used to trigger discrete domain expansion in the sequel.

<figure><img src="image_11.png" /><figcaption>**Figure 10:** Schematic of the FEM-DEM validation setup. The bridging region is depicted in red, the detection region in light blue, and the boundary DEM region in dark blue.</figcaption></figure>
<figure><img src="image_12.png" /><figcaption>**Figure 11:** Average change of neighbors per particle. The DEM domain is depicted using circular/spherical particles, and the continuum is illustrated with a black mesh. Changes in neighbors are only observed within the detection regions.</figcaption></figure>

Figure 12 shows the DEM displacements for the three previously described cases. The top row shows the situation at the beginning of the simulations, before any adaptation could occur. The bottom row displays a later situation, where the adaptivity criterion was met. It is clear that without adaptivity, the cracks are blocked in the bridging regions. This was expected as smooth displacement and deformation fields are necessary conditions for the coupling to operate well, which creates fictive elastic forces opposing crack's motion towards the continuum. Another visible consequence is the appearance of erratic FEM deformations generated from coupling a DEM region containing voids with a continuum FEM unable to represent it. When adaptive coupling is employed, the cracks follow consistent paths, independently of the initial DEM size.

Figure 13 presents a measure of the third-body layer (TBL) thickness as a function of time, which was computed as follows. First, the initial volume fraction occupied by particles in the bulk DEM region is computed. The volume fractions in several horizontal slabs of particles are evaluated (each band of particles have a vertical height of $d_0$). When the change in the current volume fraction is greater than 5% of the initial bulk, such a band of particles belongs to the TBL. It allows to characterize the boundary of the TBL region, during the entire simulation. Both small and large adaptive cases show a similar TBL thickness evolution. Although the adaptive coupling cases exhibit a similar behavior, they are not perfectly identical, which results from discrete-continuum mis-matches of two kinds. First, the particles' non-linear forces may deviate from the small deformation linear-elastic constitutive law employed for finite elements, leading to subtly distinct stiffnesses and evolution. Secondly, the amorphous arrangements of DEM particles will differ in the two systems. Indeed, the insertion of new particles during the adaptation (see steps 2 and 3 in figure 6) will bring slightly distinct particle positions, leading to a stochastic influence on system evolution. Still, the non-adaptive coupling with a small DEM region is clearly nonphysical and constraints the TBL to remain small. Furthermore, a non-monotonic evolution can be observed with a decrease after 10'000 timesteps, suggesting a totally spurious re-compaction in the DEM region, as well as a debris size reduction. This difference is necessarily due to the presence of the bridging region which prevents large deformations (crack and disorder) to penetrate the continuum. Thus, an adaptive coupling is crucial to obtain a proper dynamical evolution of such three-body layered systems.

# 5 Third-body abrasive wear with lubricating stiff elliptic bodies

It is known that stiff particles at a sliding interface result in abrasive wear (Woldman et al. 2012). In this section, we explore the abrasive or lubricative role of stiff elliptical particles trapped at a sliding interface. The elliptical particles do not exhibit adhesion with the surrounding particles.

The considered system has dimensions $200d_0 \times 200d_0 \times 3d_0$. Periodic conditions in x and z directions are used. FEM elements sizes are $10d_0$ in both x and y directions, with two elements forming the thickness along z axis. The bridging and boundary regions are chosen to match the
<figure><img src="image_13.png" /><figcaption>**Figure 12:** On the left, a non-adaptive coupled system is shown with a DEM height of 50d₀. In the middle, an adaptive coupling system with the same DEM height of 50d₀ is depicted. On the right, a larger adaptive coupling system with a DEM height of 70d₀ is presented. Among these, only the middle system experiences a top and bottom coupling adaptation. On the top, the systems are illustrated before any expansion of the discrete region for the middle system, while the bottom shows the systems after the discrete region has been expanded. The horizontal dashed black line marks the largest discrete domain upper and lower limits, showing how the top and bottom coupling regions adapted.</figcaption></figure>

<figure><img src="image_14.png" /><figcaption>**Figure 13:** Third-body thickness evolution for an initial junction size of 30d₀. Three scenarios are analyzed: a small DEM (50d₀) with a non-adaptive coupling, an initially small DEM (50d₀) with an adaptive coupling, and a large DEM (70d₀) with an adaptive coupling.</figcaption></figure>
<figure><img src="image_15.png" /><figcaption>**Figure 14:** Schematic of elliptic stiff bodies in the adaptive FEM-DEM coupling. The top particles are shown in gray (particle type: 1), bottom particles in blue (particle type: 2), and elliptical stiff bodies in red (particle type: 3).</figcaption></figure>

height of a single FE element. The initial DEM height is $70d_0$ including bridging and boundary regions. The top of the system is subject to a normal pressure of 500MPa and a constant shearing velocity $v = 0.03c$, with $c$ being the push-wave celerity while the bottom of the system is held fixed. The setup of the system is depicted in Figure 14.

The discrete domain is modeled with particles using the coarse-grained DEM law of Pham-Ba et al. (2022). Material properties are selected in line with Section 2, featuring a Young’s modulus $E$, a density $\rho$, and a restitution coefficient $\eta$ typical for SiO$_2$. The Poisson’s ratio $\nu$ is fixed at 0.25, which restricts the interaction forces to the normal direction (see Eq. 13 in the appendix). The critical length scale $d^*$, i.e. the junction size where debris begins to form, is $10d_0$. The tensile strength $\sigma_n$ will be varied in order to create different gouge regimes, similarly to what was done in Pham-Ba et al. (2023). The strengths $0.14E$, $0.16E$, and $0.18E$ are considered, which also affect the surface energy because $\gamma = \frac{d^*\sigma_n^2}{32E}$. Finally, the distribution of particle size is a Gaussian described by the parameters in Table 1.

In our model, elliptic bodies are assigned a Young’s modulus of $E_2 = 5E$, making them stiffer than the surrounding regions. There is no adhesion force acting between the elliptical particles and the surrounding matter, which leads them to act as a lubricant. Nevertheless, the lubricative effect is mitigated by their abrasive effect. The elliptical particles are stiff and at times plough through the surrounding material. The elliptic bodies have a major axis of $10d_0$ and a minor axis of $5d_0$. The particles are initially equally spaced, and we will investigate how their numbers influence the overall rheology. In Figure 14, top DEM particles are displayed in gray (particle type 1), bottom DEM particles are shown in blue (particle type 2), and elliptic bodies are depicted in red (particle type 3). An initial spacing of $d_0/2$ is introduced between the top and bottom DEM particles.

Finally, the evolution criterion is reached when at least 2% of the particles have change in neighbors greater than or equal to $\alpha = 3$. This is evaluated every 1500 timesteps.
<figure><img src="image_16.png" /><figcaption>**Figure 15:** (a) shows the evolution of shear stress $\sigma_{xy}$, while (b) represents the evolution of third-body layer (TBL) thickness. Several tensile strengths are considered: 0.14E, 0.16E, and 0.18E, along with either two or four elliptic stiff bodies.</figcaption></figure>

<figure><img src="image_17.png" /><figcaption>**Figure 16:** Final state of the coupled FEM-DEM system with elliptical stiff bodies placed at the interface. Continuum regions are displayed with a black grid. Top DEM particles are in gray while bottom DEM particles are in blue. The domain was extended in the x-direction using the periodic boundary conditions to enhance visualization. (a) and (b) feature a tensile strength of 14E, (c) and (d) exhibit a tensile strength of 16E, while (e) and (f) present a tensile strength of 18E. Configurations with two lubricant bodies are depicted on the left, and those with four lubricant bodies are shown on the right.</figcaption></figure>

We analyze the influence of elliptic particles on the overall friction and on the formation and evolution of a third-body. The friction is directly related to the macroscopic shear stress $\sigma_{xy}$
measured at the top FEM-boundary. The third-body thickness evolution is measured as was done in Section 2.

Figure 15(a) shows macroscopic $\sigma_{xy}$ when varying the tensile strength and the number of stiff lubricating bodies. Clearly, increasing the number of lubricating particles reduces the overall friction. The tensile strength is also varied, which leads to surprising effects. With four elliptical bodies, a higher tensile strength leads to more friction. This is expected as the surface energy (and effective adhesion) grows as $\gamma \propto \sigma_n^2$. However, surprisingly, when there are only two elliptical particles, the friction appears independent of the tensile strength/adhesion. A detailed investigation should be conducted in future work to identify the precise mechanisms leading to this unexpected behavior. Here, we only observe that the number of stiff friction-less particles leads to junction sizes that vary non-monotonously with particle number and strength. For the case of two elliptical particles, the contact junction sizes (between gray and blue regions) appear to be a decreasing function of strength.

In Figure 15(b), the third-body thickness evolution is shown. For low tensile strengths, $\sigma_n = 0.14E, 0.16E$, no clear effect of the number of stiff particles on the gouge thickness can be measured, as it remains rather constant in these cases. From Figure 16 it becomes visible that for low tensile strengths a mixed regime was always created, with particles attaching, detaching, and agglomerating without forming a clear wear debris. Also, the interface gap stays rather small, which is consistent with observing small third-body thicknesses.

The only case with a clear increase of the third-body thickness is for the largest tensile strength (0.18E) and two stiff elliptical particles. A salient feature of this case is that it meets the conditions for crack initiation, propagation and therefore debris formation, see Figure 16(e). Remarkably, the same strength but with four stiff particles does not lead to third-body growth.

**Our model, controls only the number of stiff lubricating particles and material strength, which leads to a complex and rich evolution of the third-body and frictional behavior.** Our adaptive method has the potential to unleash more unexpected tribological behavior, as it allows a seamless transfer of fields between discrete and continuum domains. This is seen in Figure 17, which shows the velocity fields for the highest tensile strength 0.18E for two and four stiff particles, and is displaying good continuity through the discrete-continuum coupling interfaces.

# 6 Conclusion

We developed an adaptive Finite Element Method-Discrete Element Method (FEM-DEM) coupling approach that enables the discrete region to dynamically expand in response to third-body deformation. While this method was developed to address tribology systems at engineering and geology-lab scales, it is also suitable for molecular assemblies, in both crystalline and amorphous cases.

This coupling employs the bridging method, initially designed for regular lattices and later extended to amorphous materials. Changes in the neighborhood per particle within a detection zone placed ahead of the bridging region trigger the advancement of the discrete domain followed by coupling relaxation and interpolation steps to equilibrate the numerical approach.

The adaptive coupling was tested on amorphous and crystalline materials. Wear debris simulations showed that without domain adaption cracks are arrested in the bridging region, leading to invalid mechanical behavior. This demonstrated the necessity and accuracy of the proposed adaptive coupling approach. We also investigated a tribological system with stiff elliptic particles at the contact interface. The simulations revealed the influence of hard wear particles as well as the material's tensile strength on the evolution of the third-body thickness and friction.

# Acknowledgements

G.A., J.-F.M., and M.V.-L. acknowledge financial support from the Swiss National Science Foundation (Grant 200021_197152, "Wear across scales").

# Supplementary Material

The scripts to reproduce the study are available in the following git repository.
<figure><img src="image_18.png" /><figcaption>**Figure 17:** Velocity field measured during the final timestep of the coupled FEM-DEM simulation with stiff elliptical particles inserted at the interface. Black iso-lines for the velocity field are also plotted. Horizontal black markers indicate the boundary between continuum and discrete regions. To enhance visualization, x-direction periodic replicates are used to extend the domain. (a) tensile strength $\sigma_n = 18E$ and two stiff elliptical particles. (b) $\sigma_n = 18E$ and four stiff elliptical particles.</figcaption></figure>
# 7 Appendices

## 7.1 Coupling method Voisin-Legrince et al. (2022); Voisin-Legrince et al. (2024)

The FEM-DEM bridging coupling involves an overlapping region $\Omega$ where the DEM and FEM are combined. Here, we provide a summary of the method (for more details, refer to Voisin-Legrince et al. 2022). The first step in the coupling process involves defining an energy weighting that controls the influence of each domain. The global energy is then defined as a Hamiltonian:

$$H = \int_{\Omega} \alpha(X)E^{C}(X) + (1 - \alpha(X))E^{D}(X) \ dX. \qquad (2)$$

Here $E^C$ and $E^D$ are the energy densities in the continuum and discrete regions, respectively. These are weighted using a scaling function $\alpha(X)$ which defines the influence of the continuum domain for a given spatial location $X$. However, the energy weighting does not ensure kinematic consistency between the coupled models. To address this, Lagrange constraints are employed to link the two models. For a strong coupling, the constraint applied to the particles is expressed as:

$$g = N^T u - d = 0. \qquad (3)$$

Here $u$ and $d$ denote the displacement of the nodes and particles, respectively, while $N$ represents the standard finite-element shape functions evaluated at the positions of each particle. This means that each particle is constrained to an interpolated position obtained from the finite element discretization. Both the energy weighing and the Lagrange constraint lead to a Lagrangian $H_L$:

$$H_L = H + \lambda^T g, \qquad (4)$$

where $H$ is the weighted Hamiltonian and $\lambda$ the Lagrange multiplier vector. This modified Hamiltonian leads to new equations of motion:

$$\left\{ \begin{array}{l} \bar{M}\ddot{u} = \bar{F} + \lambda \cdot \frac{\partial g}{\partial u} \\ \bar{m}\ddot{d} = \bar{f} + \lambda \cdot \frac{\partial g}{\partial d} \end{array} \right. . \qquad (5)$$

In these equations, $u$, $F$, and $M$ represent the displacements, forces, and masses of the continuum domain, while $d$, $f$, and $m$ represent the displacements, forces, and masses of the discrete domain. Through energy weighing, the masses and forces are altered, giving rise to $\bar{m}$, $\bar{f}$ and $\bar{M}$, $\bar{F}$. Considering the Verlet time integration scheme Allen et al. (1989) this results in solving:

$$\dot{g} = \Lambda\lambda, \qquad (6)$$

where $\Lambda$ is the constraint matrix, which is determined by the constraint function $g$. Therefore, the velocities of the nodes and particles in the overlapping region are corrected using the Lagrange multiplier vector $\lambda$.

In addition to the coupling region, a boundary/pad region, sets the boundary conditions for the discrete domain. Here, the displacement of particles $d$ are constrained by the displacement of the FEM node $u$,

$$d = N^T u, \qquad (7)$$

with $N$ the standard finite-element shape function.
## 7.2 Contact law from Pham-Ba et al. 2022

The force $F$ acting between a given pair of particles is the sum of a normal component $F_N$, a tangential component $F_T$, and velocity damping components $c_N v_N$ and $c_T v_T$, where $v_N$ and $v_T$ are the relative velocities in the normal and tangential directions, and $c_N$ and $c_T$ are the damping factors. We note $n_N$ and $n_T$ the units vectors in the normal and tangential directions:

$$F = -(F_N + c_N v_N)n_N - (F_T + c_T v_T)n_T \quad (8)$$

The normal component depends on the normal overlap $\delta_N$ between the particles and is expressed as:

$$F_N = \begin{cases} k_N \delta_N & \text{if } \delta_N \le \delta_e, \\ -\frac{k_N \delta_e}{\delta_f - \delta_e} (\delta_N - \delta_f) & \text{if } \delta_e < \delta_N \le \delta_f, \\ 0 & \text{if } \delta_N > \delta_f, \end{cases} \quad (9)$$

with $k_N$ the Hookean stiffness, $\delta_e$ the maximum elastic distance, and $\delta_f$ the fracture distance. The tangential force depends on the sliding distance $\delta_T$ and is defined as:

$$F_T = \min(k_T \delta_T, F'_{m,T}), \quad (10)$$

$$F'_{m,T} = \min \left( \frac{\delta_f - \delta_N}{\delta_f}, 1 \right) F_{m,T}, \quad (11)$$

with $F_{m,T}$ the maximum tangential force, and $k_T$ the tangential stiffness. The inter-particle force parameters between particles are determined based on desired material properties using the following relations:

$$m_{\text{eff}} = \frac{m_i m_j}{m_i + m_j}, \qquad r_{\text{eff}} = \min(r_i, r_j) \quad (12)$$

$$A_N = \sqrt{2r_{\text{eff}}^2} \frac{1}{1 - 2\nu}, \qquad A_T = \sqrt{2r_{\text{eff}}^2} \frac{1 - 4\nu}{(1 - 2\nu)(1 + \nu)}, \quad (13)$$

$$k_N = \frac{A_N E}{r_i + r_j}, \qquad k_T = \frac{A_T E}{r_i + r_j}, \quad (14)$$

$$c_N = \frac{2(1 - \eta)}{\pi} \sqrt{k_N m_{\text{eff}}}, \qquad c_T = \frac{2(1 - \eta)}{\pi} \sqrt{k_T m_{\text{eff}}}, \quad (15)$$

$$\delta_e = \frac{(r_i + r_j)\sigma_N}{E}, \qquad \delta_f = \frac{4\gamma}{\sigma_N}, \quad (16)$$

$$F_{m,T} = A_T \sigma_T, \quad (17)$$

with $r_i$ and $r_j$ the radii of the interacting particles, $m_i$ and $m_j$ their masses, $E$ the Young's modulus, $\sigma_N$ the tensile strength, $\sigma_T$ the shear strength, $\gamma$ the surface energy, and $\eta$ the restitution coefficient. These relations are valid for radii of interacting particles, which respect:

$$r_i + r_j \le d_c, \text{ with the critical diameter } d_c = \frac{4\gamma E}{\sigma_N^2}. \quad (18)$$

Additionally, if a particle is smaller than the neighbor size $\delta_f$, it will interact with particles beyond its closest neighbors, extending the interaction range. Therefore, to avoid an important increase in computational cost, the minimum particle's size should be:

$$d_{\min} = \frac{4\gamma}{\sigma_N} \quad (19)$$
## References

Aghababaei, R., D. H. Warner, and J.-F. Molinari (2016). "Critical length scale controls adhesive wear mechanisms". *Nature communications* 7.1, p. 11816

Ahrens, J., B. Geveci, and C. Law (2005). *An End-User Tool for Large Data Visualization, Visualization Handbook*. Elsevier

Allen, M. P. and D. J. Tildesley (1989). *Computer Simulation of Liquids*. Oxford University Press

Anciaux, G. (2024). *LibMultiScale*. https://gitlab.com/libmultiscale/libmultiscale

Avci, B. and P. Wriggers (2012). "A DEM-FEM Coupling Approach for the Direct Numerical Simulation of 3D Particulate Flows". *Journal of Applied Mechanics*. DOI: 10.1115/1.4005093

Bouillanne, O., G. Mollon, A. Saulot, S. Descartes, N. Serres, G. Chassaing, and K. Demmou (2024). *Wear in Progress: How Third Body Flow Controls Surface Damage*. preprint. DOI: 10.21203/rs.3.rs-3890447/v1

Broughton, J. Q., F. F. Abraham, N. Bernstein, and E. Kaxiras (1999). "Concurrent coupling of length scales: Methodology and application". *Physical Review B*. DOI: 10.1103/PhysRevB.60.2391

Chaudry, M. A., C. Woitzik, A. Düster, and P. Wriggers (2022). "A multiscale DEM-FEM coupled approach for the investigation of granules as crash-absorber in ship building". *Computational Particle Mechanics*. DOI: 10.1007/s40571-021-00401-5

Chen, P. Y., M. Chantharayukhonthorn, Y. Yue, E. Grinspun, and K. Kamrin (2021). "Hybrid discrete-continuum modeling of shear localization in granular media". *Journal of the Mechanics and Physics of Solids*. DOI: 10.1016/j.jmps.2021.104404

Dhia, H. B. and G. Rateau (2005). "The Arlequin method as a flexible engineering design tool". *International Journal for Numerical Methods in Engineering*. DOI: 10.1002/nme.1229

Durán, O., N. P. Kruyt, and S. Luding (2010a). "Analysis of three-dimensional micro-mechanical strain formulations for granular materials: Evaluation of accuracy". *International Journal of Solids and Structures*. DOI: 10.1016/j.ijsolstr.2009.09.035

Durán, O., N. P. Kruyt, and S. Luding (2010b). "Micro-mechanical analysis of deformation characteristics of three-dimensional granular materials". *International Journal of Solids and Structures*. DOI: 10.1016/j.ijsolstr.2010.04.014

Falk, M. L. and J. S. Langer (1998). "Dynamics of viscoplastic deformation in amorphous solids". *Phys. Rev. E* 57 (6), pp. 7192-7205. DOI: 10.1103/PhysRevE.57.7192

Finnis, M. W. and J. E. Sinclair (1984). "A simple empirical N-body potential for transition metals". *Philosophical Magazine A*. DOI: 10.1080/01418618408244210

Gao, K., B. J. Euser, E. Rougier, R. A. Guyer, Z. Lei, E. E. Knight, J. Carmeliet, and P. A. Johnson (2018). "Modeling of Stick-Slip Behavior in Sheared Granular Fault Gouge Using the Combined Finite-Discrete Element Method". *Journal of Geophysical Research: Solid Earth*. DOI: 10.1029/2018JB015668

Miller, R. E. and E. Tadmor (2002). "The Quasicontinuum Method: Overview, applications and current directions". *Journal of Computer-Aided Materials Design*. DOI: 10.1023/A:1026098010127

Mollon, G. (2018). "A unified numerical framework for rigid and compliant granular materials". *Computational Particle Mechanics*. DOI: 10.1007/s40571-018-0187-6

Mollon, G. (2019). "Solid Flow Regimes Within Dry Sliding Contacts". *Tribology Letters*. DOI: 10.1007/s11249-019-1233-0

Mollon, G. (2021). "The soft discrete element method". *Granular Matter*. DOI: 10.1007/s10035-021-01172-9

Nguyen, N. H., H. H. Bui, and G. D. Nguyen (2020). "An approach to calculating large strain accumulation for discrete element simulations of granular media". *International Journal for Numerical and Analytical Methods in Geomechanics*. DOI: 10.1002/nag.3076

Pham-Ba, S. and J.-F. Molinari (2022). "Adhesive wear with a coarse-grained discrete element model". *Computer Methods in Applied Mechanics and Engineering*. DOI: 10.1016/j.cma.2022.115124

Pham-Ba, S. and J.-F. Molinari (2023). "Role of Minimum Adhesive Wear Particle Size in Third-Body Layer Properties". *Tribology Letters*. DOI: 10.1007/s11249-023-01772-x
Rudd, R. E. and J. Q. Broughton (1998). "Coarse-grained molecular dynamics and the atomic limit of finite elements". *Physical Review B*. doi: 10.1103/PhysRevB.58.R5893

Tadmor, E. B., M. Ortiz, and R. Phillips (1996). "Quasicontinuum analysis of defects in solids". *Philosophical Magazine A*. doi: 10.1080/01418619608243000

Tadmor, E. B. and R. E. Miller (2014). *Modeling Materials: Continuum, Atomistic and Multiscale Techniques*. Cambridge: Cambridge University Press

Voisin-Leprince, M., J. Garcia-Suarez, G. Anciaux, and J.-F. Molinari (2024). "Two-scale concurrent simulations for crack propagation using FEM-DEM bridging coupling". *Computational Particle Mechanics*. doi: 10.1007/s40571-024-00788-x

Voisin-Leprince, M., J. Garcia-Suarez, G. Anciaux, and J.-F. Molinari (2022). "Finite element method–discrete element method bridging coupling for the modeling of gouge". *International Journal for Numerical Methods in Engineering*. doi: 10.1002/nme.7171

Wagner, G. J. and W. K. Liu (2003). "Coupling of atomistic and continuum simulations using a bridging scale decomposition". *Journal of Computational Physics*. doi: 10.1016/S0021-9991(03)00273-0

Wellmann, C. and P. Wriggers (2012). "A two-scale model of granular materials". *Computer Methods in Applied Mechanics and Engineering*. doi: 10.1016/j.cma.2010.12.023

Woldman, M., E. van der Heide, D. Schipper, T. Tinga, and M. Masen (2012). "Investigating the influence of sand particle properties on abrasive wear behaviour". English. *Wear* 294-295, pp. 419–426. doi: 10.1016/j.wear.2012.07.017

Xiao, S. P. and T. Belytschko (2004). "A bridging domain method for coupling continua with molecular dynamics". *Computer Methods in Applied Mechanics and Engineering*. doi: 10.1016/j.cma.2003.12.053

Zhang, B. and R. A. Regueiro (2015). "On large deformation granular strain measures for generating stress–strain relations based upon three-dimensional discrete element simulations". *International Journal of Solids and Structures*. doi: 10.1016/j.ijsolstr.2015.04.012

Zhang, Y., G. Mollon, and S. Descartes (2020). "Significance of third body rheology in friction at a dry sliding interface observed by a multibody meshfree model: Influence of cohesion between particles". *Tribology International*. doi: 10.1016/j.triboint.2020.106188