Concurrent multi-scale modeling of granular materials: Role of coarse-graining in FEM-DEM coupling

The finite element method (FEM) is commonly used for modeling continuum media, while particle simulation methods like the so-called discrete element method (DEM) are used for discrete systems. Coupling the discrete (DEM) and continuum (FEM) methods is conventionally achieved through a direct mapping between discrete particles and finite elements. Coarse-graining (CG) is a micro–macro transition (discrete-to-continuum) method that maps discrete particle data onto smooth, differentiable fields that satisfy the continuum equations. By choosing an appropriate length scale (the coarse-graining width c ), the coarse-grained fields are then homogenized and projected onto a FEM spatial discretization. This concept is utilized here to reformulate FEM-DEM coupling methods, both surface and volume, where in the limiting case of c → 0, the classical coupling is recovered. For surface coupling, the discrete particle–surface contact forces are first mapped onto a continuous surface traction field (using CG) which is then coupled to the continuum FEM model. For volume coupling (also known as the Arlequin framework), the homogenization operators are enriched with CG functions, offering a non-local coupling approach between discrete particles, their continuum fields, and the finite elements. The CG enrichment represents a new strategy that consists of (1) a discrete-to-continuum mapping and (2) a continuum-to-continuum coupling based on “CG-enriched homogenization” (CGH). It is shown for surface coupling that the CG-enriched formulation not only leads to more accurate results, conserving symmetry, but also reduces energies generated by the coupling. For volume coupling, there is consistently less numerical dissipation with than without CG-enrichment, especially when the dynamic load contains high-frequency content. Finally, the optimal CG widths are identified for very simple test cases, with which the surface or volume coupling performs best. CGH can be potentially extended beyond the present examples, by considering other continuum fields (e.g., higher-order) and equations (e.g., multi-physics), and used to formulate other multi-scale modeling methods.


Introduction
Granular materials are ubiquitous in nature and industry, e.g., from the geological units (e.g., sands and rocks) that constitute the Earth's subsurface [1,2] to pharmaceutical powders and agricultural grains [3][4][5]. The composition of these particulate materials can be highly heterogeneous, with their bulk properties arising from interparticle collisions and microstructural evolution [6]. The emergent mesoscopic and macroscopic behavior can therefore be quite complex, even for one physical process alone (e.g., dilatancy [7], plasticity [8], anisotropy [9], and ratcheting [10] in the solid state and non-locality [11], shear thickening/thinning [12], and segregation [13] in the fluid state). The full spectrum of this complex behavior cannot be easily unified in a single numerical framework.
In a continuum description, a granular medium can be treated as a continuum body in which a constitutive law or closure relation is essential to solve the material's boundary value problem. For example, constitutive laws for densely packed granular materials, like soils and rocks, are developed within the framework of critical state soil mechanics, based on advanced plasticity theories in which microstructural evolution such as fabric anisotropy [14], particle breakage [15], etc., are taken into account. On the other hand, fast granular flow can be effectively modeled by incorporating a granular rheology (e.g., µ(I ) [16]) or by rewriting the yield functions of soil plasticity as a non-Newtonian viscosity [17,18] in the Navier-Stokes equations.
The Discrete Element Method (DEM), also known as the Discrete Particle Method, is commonly used to simulate granular materials [19]. It calculates the motion of individual particles by applying forces and torques that stem from either external body forces or particle-particle interactions. The Discrete Element Method can be enriched by coupling it to a continuum solver, resolving two major limitations: Firstly, DEM is too inefficient to simulate large quantities of granular material. Ideally, this should be resolved by using an up-scale continuum method, efficient enough to simulate large-scale, industrial machinery. However, all available constitutive models for granular fluids or solids are based on simplifying assumptions that limit their range of applicability. It is appealing to develop a multi-scale setting in which a continuum description is used where the constitutive equations are known and DEM is used where they are not known.
Secondly, DEM usually assumes that external objects that interact with particles are immobile and rigid. Although flexible discrete elements (e.g., facets and cylinders) have become available recently (for example [20][21][22]), the material behavior approximated with these elements can only be elastic, with limited accuracy. A continuum description of these external bodies under the impact of particles (leading to irrecoverable deformation) can be achieved most accurately with FEM. Particles' interaction with deformable structures may strongly influence the rheological properties of the granular bulk [11,23] and is relevant to many industrial applications, such as sand-filled geotubes [24,25], monopile installation [26,27], and clogging [28].
Concurrent multi-scale methods aim to address both challenges: "surface coupling" can be used to simulate the interaction of particles with deformable bodies, whereas "volume coupling" is well suited to simulate most of the granular material as a continuous bulk, and only use DEM in regions where the continuum model does not suffice. These two approaches can also be combined, as shown in Fig. 1.
The most common approach to couple FEM and DEM is surface coupling. As the name suggests, the exchange of information occurs at interacting surfaces of multiple domains [29,30]. In volume coupling (also known as the Fig. 1. Example of a concurrent multi-scale simulation. Ω FE 1 denotes a deformable structure, Ω DE the domain of the discrete particle simulation, and Ω FE 2 the domain of a continuum model for the granular bulk. The surface interactions between the particles and the deforming structure are handled at Γ C , while a volume-transition from a discrete to a continuum model is achieved in Ω C . Arlequin framework [31,32]), the computational models are coupled within an overlapping volume by applying kinematic constraints on their solutions and weighting the governing equations based on the partition of unity. Originally proposed by Dhia and Rateau [31] for coupling FEM models with different mesh resolutions, the Arlequin framework was extended by Wellmann and Wriggers [26] to couple FEM and DEM for concurrent multi-scale modeling of granular materials. Recently, Yue et al. [33] developed a similar hybrid discrete-continuum framework in which the DEM model is coupled with the material point method (MPM) model for simulating free-flowing granular materials. Another approach that has been popular recently is the hierarchical multi-scale framework [34][35][36] where the material's constitutive model is substituted by DEM simulations at the level of integration points.
A fundamental challenge in FEM-DEM coupling is how to map between microscopic and macroscopic scales such that information lost and numerical artifacts are minimized. In these aforementioned studies and those cited therein, the coupling between DEM and FEM is generally achieved through a direct mapping between discrete particles and finite elements (also referred to as homogenization), via FEM basis functions or volume averaging. Conventionally, discrete data, such as particle displacements and contact forces, is coupled to the continuum-scale governing equations as point sources. In this work, we take a two-step approach that consists of (1) a discreteto-continuum mapping using "coarse-graining" (CG) [37,38] and (2) a continuum-to-continuum coupling based on "CG-enriched homogenization" (CGH). CG is here not referring to coarser, bigger particles, but is a micro-macro transition technique that extracts (locally varying) continuum quantities from discrete particle data by applying a spatial smoothing length scale [37], the so-called CG width.
For surface coupling, we define a smooth traction field, and for volume coupling, a continuous particle velocity field, using CG. The "coarse-grained" traction over interacting surfaces is passed to the FEM mesh as an updated boundary condition, whereas the "coarse-grained" particle velocity is constrained to the FEM velocity over an overlapping volume. Using CG, the coupling terms are derived in a more generalized form, containing an internal length scale parameter, i.e., the coarse-graining width. In the limiting case of a zero CG width, the CG-enriched formulations reduce to the most commonly used ones in [26,39].
To the authors' knowledge, this is the first work that uses coarse-graining as the basis to formulate multi-scale coupling problems. The benefits of coarse-graining to surface coupling are demonstrated through a simple numerical example, that is a pair of particles sliding over an elastic cantilever. Through another numerical example, namely, wave propagation between discrete and continuum descriptions of an idealized granular material, we demonstrate the benefits of coarse-graining to volume coupling. For each application, an optimal length scale can be found that is large enough to minimize discretization errors and small enough to preserve the physical gradients of the resulting fields [38].
The paper is organized as follows: Sections 2 and 3 summarize the key theoretical aspects of the particle and continuum models using variables defined at the respective scales. We provide a brief overview of coarse-graining in Section 4 and show how to use it as a versatile tool to reformulate the surface and volume coupling problems in Sections 5.1 and 5.2. The algorithmic details can be found in Section 5.3. The numerical examples for verification study are introduced in Sections 6.1 and 6.2. Section 7 summarizes the concluding remarks and outlook. In the following, we use boldface letters for vectors and matrices. Greek subscripts α, β, and γ are used to denote particle-scale variables; lowercase letters i, j and k are used for the finite element discretization of the continuum description.

The microscopic particle model
The Discrete Element Method models granular materials as assemblies of rigid particles that interact via binary interparticle contact forces and torques. Simplified, rigid geometric objects (walls) can be added that also interact with the particles via binary contact forces and torques. External forces and torques can be applied to the particles as well. The macroscopic behavior of the particle assembly is then obtained by resolving the translational motion of individual particles according to Newton's second law and the rigid-body rotation via Euler's equations (rotational rigid body dynamics [40]). Consider particle α that has N and N w existing contacts from the neighboring particles and walls, 1 for given initial and boundary conditions, the time evolution of the velocity v α , position x α , angular velocity ω α , and orientation q α of a particle P α is given by where m α and I α are particle's mass and inertia tensor, C(q α ) is the transformation matrix that allows an efficient handling of its orientation q α as quaternion, f b α are the body forces, and f α and τ α are the total forces and torques acting at the particle position x α . Eqs. (1) and (2) are closed using so-called contact laws that define interparticle contact forces f αβ between a pair of interacting particles P α and P β . Another contact law, similar to, or different from the interparticle contact law, can be used to compute particle-wall interaction forces f w αγ . The branch vector l αβ = x c αβ − x α connects the particle center x α with the contact point x c αβ , and the cross product l αβ × f αβ contributes a torque that affects the particle rotation and angular velocity [41]. In its simplest form, a linear-elastic forcedisplacement law relates the contact force in the direction normal to the contact area and the interparticle overlap u n (resembling deformation) by a linear spring constant, namely, ∥ f n αβ ∥ = k n u n . Similar force-displacement laws are often taken in the tangential direction with an added Coulomb yield criterion, ∥f s αβ ∥ ≤ µ∥f n αβ ∥, in order to capture irrecoverable plastic deformation, using a friction coefficient µ. For more details on all forces see [41].
The above differential equations are solved numerically using the Velocity-Verlet algorithm. Using a higher-order accurate time integration scheme would not increase accuracy, because most contact laws used in the DEM have discontinuous derivatives. The open-source DEM code MercuryDPM [42] is used for the particle side of our new DEM-FEM coupled code.

The macroscopic continuum model
We will now introduce the FEM model used for modeling deformable bodies. We restrict ourselves to describing the deformation of an elastic body; however, this can be changed without impacting the coupling methods. We solve for the motion of this body using a continuum approach for given initial and boundary values. oomph-lib -an open-source Object-Oriented Multi-Physics Finite-Element Library [43] -is used for the macroscopic modeling. The governing equations are briefly introduced below, for the sake of completeness.

Governing equations
We assume that the deformable body behaves as a three-dimensional elastic solid and describe its behavior using finite strain theory. A finite-strain framework is essential to correctly handle the impact of granular materials, in particular for deformable structures that have large displacements and/or geometrical nonlinearity. The vectors ξ and X, attached to the body, are used to define the material's undeformed and deformed configurations. A strain measure, the Green strain tensor, is defined as ε = 1 2 (F T · F − 1), where F = 1 + ∇u is the deformation gradient w.r.t. ξ , u = X − ξ the displacement field, 1 the identity matrix, and (·) T the transpose.
Let us consider a solid body subjected to a surface traction t on a subset Γ t of the body's boundary ∂Ω FE , a body force density b acting on the domain Ω FE , and a displacement boundary condition, prescribed by X (BC) on another subset of the boundary Γ X . Using the stress and strain measures σ and ε, the deformation is governed by the principle of virtual displacements [44] ∫ where ρ is the density of the undeformed body, δ(·) represents the variation of (·), and ":" denotes a double contraction between two rank-two tensors. Eq. (3) must be closed with a constitutive law that relates the stress and strain tensors. This closure relation can be the elasticity tensor C, via dσ = C : dε. 2 The symmetric second Piola-Kirchhoff stress tensor σ is a measure of force per unit undeformed area and is the work conjugate to ε [45].
In this work, we consider only elastic material behavior, using Hooke's constitutive law (see Appendix A for details).

Finite element discretization
oomph-lib uses the total Lagrangian formulation for solid mechanics problems, based on the large-displacement, finite strain theory. In this formulation, all variables in Eq. (3) are defined as functions of the initial, undeformed configuration. For isoparametric finite elements, the local position vector X = (X 1 , X 2 , X 3 ) is interpolated using the element's basis functions ψ i (i = 1, . . . , N ), via where X pi is the pth coordinate of the position vector at nodal point i, ψ i is the corresponding interpolation or basis function, and N is the number of nodes per finite element. The discretization of Eq. (3) can be expressed by the local residual for an element e, 3 Assembling the local equations in an element-by-element fashion yields a system of global residuals, R(U), where U is the vector of all unknowns X pi . The exact solution U then satisfies R(U) = 0. This system of equations is discretized in time and solved at each time step. Let U k be the solution vector at time t k . Since the system Eq. (5) is nonlinear, the solution U k+1 at time t k+1 is found using Newton's method: Starting with an initial guess U k+1 = U k , we evaluate the residual R and the Jacobian J = ∂R ∂U at U k+1 and solve the linear system Then we update our guess, We repeat this iteration until ∥R∥ is sufficiently small. Note, in Eq. (5) we assumed thatÜ k+1 is known. In the time-discretized scheme, we use the Newmark-beta method [46] to interpolateÜ k+1 from the values of U in the previous time steps. Details about the construction of the Jacobian matrix and residual vector and the time integration scheme can be found in [43].

Coarse graining
The key issue in the coupling between particle and continuum models is to have an accurate homogenization that maps microscopic quantities (e.g., density, momentum, and stress) onto the macroscopic scale, while obeying the fundamental laws of physics, such as mass and momentum conservation, to ensure minimum information loss. In what follows, the theoretical basis of CG will be briefly reviewed.
The aim of coarse-graining is to define continuum fields (ρ, v, and σ ) that locally satisfy mass and momentum conservation. By applying a spatial smoothing kernel W(x) to the statistical-mechanics definition of point density ρ m (x) = ∑ N α=1 m α δ(x − x α ), the coarse-grained density is obtained as 2 Choosing the position vector in the deformed configuration as the unknown variable allows the variation of deformation δε to be expressed as δε = 1 . 3 Here, we took advantage of the fact that the 2nd Piola-Kirchhoff stress tensor is symmetric. with • denoting a convolution and δ(·) the Dirac delta function. We require that the volume integral over density is equal to the mass of all particles, ρ(x) is non-negative, and for efficiency reasons, only particles within a small enough distance of x α from x (c in Eqs. (9)-(10)) contribute to the density ρ(x). Therefore, the CG function W in a n-dimensional space: • is normalized: • has compact support: ∃c ∈ R: W(x) = 0 for all |x| > c.
Two typical CG functions are the cut-off Gaussian and the Lucy polynomial [42]: where C G and C L are the appropriate prefactors for normalization. Note, for spherical particles, W is isotropic in space and c is the only parameter that controls the support of a CG field [47].
Extending the coarse-grained definition of density to velocity and stress [48] allows us to define the continuum fields v and σ from particle-scale quantities, using any CG function W. For example, to satisfy the mass and momentum conservation laws, velocity and momentum density are defined as This coarse-grained velocity field will be constrained to the continuous velocity field interpolated from the finite elements to derive the FEM-DEM volume coupling formulation.
To account for the particles' interactions with external boundaries, Weinhart et al. [48] proposed to add an additional term to the momentum balance, the boundary force density, where x c αγ is the γ th contact point and f w αγ the interaction force between particle α and the rigid wall. Note that b BC is only non-zero in the vicinity of the boundary Γ C . For coarse-graining functions that vanish in the direction perpendicular to the boundary Γ C , b BC is equivalent to a surface traction defined on Γ C , where W Γ C is a CG function evaluated on the boundary Γ C and normalized to satisfy We will use this coarse-grained surface traction to couple the DE and FE formulations for surface coupling problems.

CG-enriched formulation for concurrent multi-scale modeling
Within the concurrent multi-scale modeling framework, the computational domain Ω is divided into multiple subdomains; in each subdomain the material behavior is described by a different sub-model. In this work, DEM and FEM are the two methods used to model the material behavior arising from a particulate or continuum nature, thus Ω = Ω FE ∪ Ω DE . This multi-scale model is set up "concurrently", with its sub-models synchronized in time and coupled either on an interacting surface Γ C = ∂Ω FE ∩ ∂Ω DE (surface coupling) or in an overlapping volume Ω C = Ω FE ∩ Ω DE (volume coupling). Differently from conventional approaches in which the finite element basis functions are used for micro-macro projection, we formulate the multi-scale problems with coarse-graining (CG), which is used to map conservative variables (i.e., mass, momentum, energy) between microscopic and macroscopic scales. In the following, we will present the CG-enriched formulations for multi-scale surface and volume coupling problems. Note, the superscripts "DE" and "FE" will be used throughout the subsections to clarify the microscopic and macroscopic nature of the computational domains, and the subscripts {α, β, γ , . . .} (particles) and {i, j, k, . . .} (nodes) to distinguish variables of the microscopic and macroscopic models.

CG-enriched formulation for FEM-DEM surface coupling
In the current implementation of surface coupling (see Fig. 2(a)), Ω DE = R 3 \Ω FE , thus particles can move anywhere in space outside the deformable object and hit at any part of the boundary ∂Ω FE . The DEM and FEM formulations introduced in Sections 2 and 3 are used to solve the governing equations in their respective subdomains. The exchange of forces and velocities on the coupling surface is done with or without a CG-enrichment for the micro-macro transition. As shown in Figs. 2(b) and 2(c), a surface coupling force f w αγ can be coarse-grained into a continuous surface traction field t SC e (r). The triangulation approach as in [42] allows for detecting particlefinite element interactions. However, an extension of the existing algorithm is needed to handle particle interactions with arbitrary concave/convex surfaces.

Contact detection for particles-triangulated-wall interactions
The coupling surface Γ C is meshed (see Fig. 3a), and the meshed surface ∂Ω FE is replicated in DEM by a set of triangle walls {T DE γ } N w γ =1 . Note that there is a choice of how to represent the meshed surface by triangle walls (e.g., we represent each rectangular surface element by four triangles). Ideally, this should not affect the calculation of the contact forces; i.e., the contact force computed between a particle and the set of triangle walls should be the same as if the meshed coupling surface was represented by a single, non-triangulated wall.
For a flat surface, this implies that the contact force between particle α and the discretized surface should be the same irrespective of whether the contact area contains a triangle edge, vertex, or face. In a naive implementation, however, if the contact area contains a triangle edge (or vertex), then contact forces would be computed with two (or more) triangle walls; thus the particle experiences a force bump when sliding over a triangle edge (or vertex). In order to remove this effect we have implemented two modifications to the contact force calculation: 1. First, naively compute all contact forces between particle α and triangles T DE γ . As each triangle is convex, at most one contact force f w αγ is computed per triangle. See [42] for the calculation of contact forces at the triangle's face, edges, and vertices, respectively. 2. Modification 1: If any two of these contact forces have an overlapping contact area (i.e., both contact areas contain a common triangle edge), and if the common edge is convex or flat, then eliminate the smaller of the two forces. 3. Modification 2: If there are still multiple contact forces between α and the coupling surface, then we assume that these contacts share a concave edge or a vertex. Thus, we multiply each force f w αγ by a weight Note that these weights add up to 1, thus, the total contact force between α and the coupling surface will be a weighted average of the individual contact forces.
The modified contact detection algorithm is shown in Algorithm 1. The resulting forces are illustrated in Fig. 3.
Modification 1 ensures that a particle feels no bump when crossing a flat edge. Modification 2 was added to ensure that the force is not sensitive to small changes in the edge angle, i.e. there is no sudden jump in force if the convexity of an edge changes due to deformations in the solid. Note that the current implementation also allows hysteretic forces to be kept between triangles, which is essential for simulating history-dependent microscopic behavior (e.g., interparticle sliding friction).
Furthermore, the potentially large number of triangles {T DE γ } N w γ =1 in the coupling region Γ C requires an efficient treatment of the detection of particle-wall contacts. The detection of particle-particle contacts is done efficiently in MercuryDPM via the hierarchical-grid algorithm [49,50]: Particles are sorted into a set of L spatial grids with different cell sizes s 1 < s 2 < · · · < s L . Each particle P is mapped into a specific cell C ln containing all particles of diameter d ∈ [s l−1 , s l ] (with s 0 = 0) and position x j ∈ [(n j − 1 2 )s l , (n j + 1 2 )s l ], j = 1, 2, 3. A particle P ∈ C ln can only have contact with particles in the same or neighboring grid cells of the same level l, or with nearby grid Algorithm 1: Contact force calculation, with modifications 1 and 2 highlighted in blue and orange for all walls W do for all particles P close to wall W do if particle P and wall W overlap then add (P, W ) to the list of interactions I for all walls V in the same group of walls as W do if (P, V ) ∈ I and the contact areas of (P, V ) and (P, W ) overlap then if the angle between V and W is flat or concave then if the overlap of (P, V ) is larger than the overlap of (P, W ) then remove (P, V ) from the list of interactions I else remove (P, W ) from the list of interactions I end end multiply the force f and torque t by w = |f (P,V ) |/F end end cells on a different level. Thus, the grid can be used to efficiently compute all interparticle contacts. The algorithm's complexity has been proved to be O(N ) for realistic wide size distributions [51].
However, many wall types in MercuryDPM (semi-infinite planes, infinitely long cylinders, etc.) are unbounded, and thus cannot be placed in such cells. Thus, this efficient method of contact detection is not used for wall-particle interactions. Detecting particle-wall contacts, therefore, requires checking the overlaps of all particles and all walls. The cost of this computation is O(N w · N ), making simulations with many walls inefficient. Triangle walls, however, have a finite size, and thus can make use of the hierarchical grid, as shown in Algorithm 2: First, we compute a bounding box, the range of x, y, and z values in which the triangle resides, then we look for overlaps with particles in all grid cells nearby the bounding box. Thus, the complexity or wall-particle contact detection is reduced to O(N w ). This has been implemented in MercuryDPM for the purpose of this paper and can be used for any wall type for which a bounding box can be defined.

The surface coupling forces
We denote all finite elements {Ω e } N e e=1 part of whose boundary is on the meshed coupling surface ∂Ω FE as surface coupling (SCoupling) elements. The contact surface Γ e = ∂Ω e ∩ ∂Ω FE , belonging to a SCoupling element Ω e is discretized in the DEM sub-model by a set of interconnected triangle walls , with N w e the total number of triangle walls per SCoupling element. The positions and velocities of the triangle vertices are computed by interpolating values from the SCoupling elements. At every time step, a list of particles P α (α = 1, . . . , N p γ ) are found to be in contact with a given discrete triangle T DE γ . For each particle-wall pair, the interaction force f w αγ at a contact point x c αγ is computed from the overlap via a contact law. To map the interaction forces onto the finite element formulation, we coarse-grain all particle-wall contact forces with a given CG width. Thus, according to where W ∂Ω FE is a CG function evaluated on ∂Ω FE and normalized to satisfy ∫ ∂Ω FE W ∂Ω FE (x − x c αγ ) dA = 1 (in practice, it will be the quadrature formula that is normalized). The CG width c is the only user-defined parameter. We will show later that there is a systematic way to choose c that minimizes the error (in comparison with a given analytical solution) and any other numerical artifacts that violate the fundamental laws of physics.
To compute the nodal contribution to each element's residual, the coarse-grained field t SC (x) is multiplied by the finite element basis functions and integrated, that is It is worth noting that the limit of f SC ei as c approaches 0 is This means that, in the limit of c → 0, Eq. (14) becomes the point-source formulation of the surface coupling force widely used in literature [30,39,52,53]. The integral in Eq. (15) can be computed numerically, using the standard Gauss quadrature rule. In this case, the CG function is evaluated at the integration points and then normalized. Eq. (16), however, is an analytical expression obtained directly using the integration property of the Dirac delta function, ]. Note, in Eq. (16), the summation over all SCoupling element is gone. This is because with the Dirac delta (c = 0) the contact forces only contribute to the SCoupling elements that they are in contact with.
With coarse-graining, the discrete particle data is mapped to continuum fields, such that the conservation of mass and momentum are obeyed. This is exactly what is needed in surface coupling, where the discrete contact forces have to be mapped to a continuous surface traction field -discretized at the FEM mesh/integration points. The difference with the conventional approach is that it provides a smooth, mechanically meaningful definition of the traction field from a collection of discrete particle-wall contacts, even for very few contacts. Both formulations are implemented in the coupled MercuryDPM-oomph-lib code, as detailed in Section 5.3. The influence of CG functions (including the Dirac delta) on the surface coupling will be investigated in Section 6.1.

The Arlequin framework
The volume coupling method, also known as the Arlequin framework for multi-model coupling [31], uses a Lagrange multiplier or penalty-based approach to enforce kinematic constraints in an overlapping domain Ω C .
Within Ω C , it is required that the virtual work amounts to the weighted sum of the work respectively done by the sub-models, that is δW = δW FE + δW DE . The virtual work done in the DEM model δW DE is weighted by a smooth and continuous function w(x) of particle position x ∈ Ω C , which monotonically increases from zero on ∂Ω DE ∩ Ω C to one on ∂Ω FE ∩ Ω C , as shown in Fig. 4(a) and later in Fig. 18(b). The virtual work done in the FEM model δW FE is weighted by the function 1 − w(X) of the FEM position vector, such that the coupling weights on the FEM and DEM sides sum to unity at any given location. The coupled governing equations also need to be modified using the respective coupling weights, as demonstrated in Section 5.2.3.

CG-enriched homogenization
Before applying kinematic constraints on particle-or continuum-scale variables (e.g., displacement or velocity), a computational homogenization is needed to bring them onto the same scale, as illustrated in Fig. 4(b). The fact that DEM and FEM give solutions at the particle and continuum scales makes their volume coupling non-trivial: the length scale (either microscopic or macroscopic) on which the variables should be constrained is not well defined. A novel approach we take here is CG-enriched homogenization (CGH), using smooth kernel functions whose length scale can be chosen independent of particle and finite element sizes. In the following, we will demonstrate how CG is used to enrich a computational homogenization and in which limit this CG-enriched formulation reduces to the conventional one.
Let us now consider N p particles that reside within a finite element in the coupling domain Ω C e , as illustrated in Fig. 5. Using coarse-graining, we can define a homogenized velocity field v CG (x) at any position where m α is the mass, v α the velocity and x α the position of particle α, and W a coarse-graining function of width c. The definition of W can be found in [37,38]. Note, due to the limited effective range of the coarse-graining function, only particles near the evaluation point x have to be considered to compute the sum. We assume the finite element basis functions ψ i (i = 1, . . . , N ) to be orthogonal and, by convolving the velocity field v CG with ψ i (i = 1, . . . , N ), obtain the particle velocity at the nodal position x i as, This allows defining a continuous velocity field , for x ∈ Ω C . Note the similarity between Eq. (18) and the CG-enriched surface coupling force in Eq. (15).
Similar to the surface coupling force, with a vanishing coarse-graining width (lim c→0 W = δ), we recover the conventional form of the local "micro-macro projection", which is derived from a minimization of the constraint conditions [26], where each basis function ψ i is evaluated at the positions of the coupled particles x α ; here, α = 1, . . . , N C p denotes the particles whose positions reside strictly within the coupling region Ω C .
The major properties of the generalized formulation are: 1. The homogenization rule in Eq. (18) is derived from and compatible with momentum conservation; even the sub-particle scale fluctuations of particle motions can be coarse-grained using the CG formalism. 2. Particles not positioned inside a finite element can still contribute to the coupling for a nonzero CG width, that is to say, CGH is non-local. 3. We do not need to evaluate the FEM basis function at non-quadrature points which would require an expensive mapping between global and local coordinates.
Both homogenization approaches, Eqs. (18) and (19), are implemented in the coupled MercuryDPM-oomph-lib code, as detailed in Section 5.3. The benefits of using CG for volume coupling will be investigated and discussed in Section 6.2. Note that the same procedure can be utilized to homogenize other macroscopic variables, such as the stress and strain fields. However, formulating the constraints on these higher-order terms is beyond the scope of the current work.

Coupled governing equations
The governing equations 4 for the continuum and particle models are: where f b α is the body force acting at the particle position x α and f αβ are the contact forces acting at the contact point x c αβ , for β = 1, . . . , N α , denoting the contact partners of particle α = 1, . . . , N p . For frictional particle contacts, rotational degrees of freedom have to be considered as well. For the continuum description in Eq. (20), ρ is the density of the undeformed body, b is a body force density acting on the domain Ω , and t denotes a surface traction acting on a subset of the body's boundary, Γ t ⊂ ∂Ω FE . σ and ε are a work conjugate pair of stress and strain tensors, and δX and δε are virtual variations of the position vector and strain tensor [44].
We define the coupling weight function w(x) as a monotonic function of the position vector x in the DEM subdomain, and subsequently 1 − w(X) in the FEM subdomain, with Multiplying the governing equations above with their respective coupling weights 1 − w and w for δW FE and δW DE gives the weighted sum of the total virtual work δW = δW FE + δW DE , with where N e and N r are the total numbers of volume and surface elements in Ω and Γ t , respectively, and N p is the total number of discrete particles. We use the short-hand notation w α = w(x α ) and w αβ = w(x c αβ ) for the weights at the particle positions and contact points. δx α is a variation of the position of particle α.
Following [26], we require the difference between the FEM and DEM displacements in the last time step, u FE = X| t − X| t− dt and u DE = v DE dt, to be vanishingly small 5 at the macroscopic scale. To this end, we make use of the same coarse-graining approach as applied to the particle velocity: Note that both Eqs. (18) and (19) can be written in terms of a homogenization operator, We will now apply the same homogenization to the virtual particle displacements, δx α , to define a virtual particle displacement field, To enforce the kinematic constraint on Ω C , we penalize the velocity difference at the macroscopic scale. This gives rise to an additional term in the virtual work equation, δW = δW FE + δW DE + δC, with where ϵ is the penalty parameter. The volume coupling terms can now be derived from the virtual work equation; interested readers are referred to Appendix B for the details. Thus, the principle of virtual work for the coupled system results from Eq. (29) as with the FEM coupling force density b C and the DEM coupling force f C α defined as Finally, time integration is realized using the respective time steppers best suited for the particle and continuum models, namely, the Velocity-Verlet [42] for the former and the Newmark-beta [46] for the latter. At a given time t, the nodal displacements and the particle velocities and positions are solved explicitly for the synchronously updated boundary conditions and coupling forces. In this work, the penalty parameter ϵ is chosen to be sufficiently large (e.g., ϵd 2 = 4 × 10 8 Pa) such that the constraint condition is approached as much as possible.

Implementation of concurrent multi-scale methods in the coupled MercuryDPM -oomph-lib code
The generalized surface and volume coupling methods are implemented in the MercuryDPM branch https: //svn.mercurydpm.org/SourceCode/Branches/OomphlibCouplingExplicit. We plan to make a future version of this code a feature of the main MercuryDPM repository.
Our implementation produces a single executable that couples together the two open-source software packages, MercuryDPM and oomph-lib. Mainly written in C++, both codes use object-oriented programming techniques, including multiple inheritance, template classes, and function overloading. The flexibility of both codes allows implementing flexible wrappers around the existing finite element and discrete particle classes, with the additional member data and functions for a particular type of FEM-DEM coupling. The data necessary for both surface and volume coupling problems are handled from the lowest levels of both methods, namely, between the finite elements and their interacting particles. As illustrated in Fig. 6, the implementation requires three levels of inheritance: • the definitions of a MercuryDPM and an oomph-lib problem where the particle model and the continuum model are implemented, respectively, • the interface class that provides access and helper functions to the particle or continuum model, • and the definition of a coupled problem where the actual surface or volume coupling algorithm is implemented.

The surface-coupling element and algorithm
For the surface coupling, the standard finite element is wrapped into a surface-coupling element (see Fig. 7), adding data and functions to store and compute the residual contribution from surface-coupling. Fig. 7 shows the member data and functions of the SCoupling element derived from a template finite element class ELEMENT. In addition, the SCoupling element overrides the functions that add the elemental contributions to the local residual vector. A detailed description of the surface coupling algorithm is given below and summarized in Fig. 8. Note, in detecting a particle's interaction with a list of triangles connected at vertices and edges (a triangulated wall), only a single contact point is allowed between the particle and the triangulated wall. Moreover, the history of this contact point has to be kept intact when the particle moves from one connected triangle to another. This is ensured with the modifications made to the standard particle-triangle interaction calculation algorithm of MercuryDPM (see details in Section 5.1.1).

The volume-coupling element and algorithm
For the volume coupling, it can be seen from Eq. (29) that the governing equations need to be weighted in the dynamically adapted coupling zone. Therefore, both the finite element and discrete particle classes are extended with member functions to evaluate the local coupling weights at the integration points, particle positions, and contact points. The weighting function is assumed to be in the FEM solution space (and thus moves if the elastic body deforms); thus we first set its nodal values and then interpolate for a given location, using the FEM basis functions Fig. 9 shows the member data and functions of the VCoupling element derived from a template finite element class ELEMENT.
The VCoupling element overrides the functions that add the elemental coupling forces to the residual vector. A description of the volume coupling algorithm is given below and summarized in Fig. 10: 1. At every time step, the extended bounding box 7 of each finite element is updated. To detect whether a given (extended) bounding box overlap with discrete particles, the same contact detection algorithm, namely the hierarchical grid of MercuryDPM, is used.  2. When homogenizing a continuum velocity field from discrete particles to finite element nodes using Eq. (19), one needs to obtain the local coordinates of the particle positions in the coupled element. With CGH (see Eq. (18)), this global-to-local coordinate transformation is avoided. However, evaluating the CG function at the integration points still requires a local-to-global coordinate transformation which can be easily achieved with the FEM basis functions. 3. From the difference between the coarse-grained particle velocity and nodal velocity fields at the previous time step t − dt, the coupling forces at the microscopic and macroscopic scales are computed using Eqs. (32) and (31). 4. The last modification to the governing equations is with the coupling weights: they are computed and assigned to the coupled elements and particles, depending on the distances from the locations where they are evaluated (e.g., integration points, particle centers, and contact points) to the boundaries of the coupling zone. 5. With the updated boundary conditions and coupling forces, the velocity and position vectors in the particle and continuum subdomains are updated by MercuryDPM and oomph-lib, respectively. 6. The simulation repeats steps 1 to 5 as long as t ≤ t max is true.

Surface coupling: particle-cantilever interaction
In this section, the performance of the conventional and CG-enriched formulations are compared for surface coupling. The benefits of coarse-graining are exemplified by a simple verification case, namely, a pair of frictionless particles sliding over an elastic cantilever. This numerical example is also well suited to investigate the smoothness of the contact detection algorithm for particle-wall interactions. The material parameters of the FEM and DEM submodels are listed in Table 1. In order to compare with simple analytical solutions, both sub-models are undamped, purely elastic, without any dissipative terms (e.g., viscous and plastic) in the governing equations. It should be noted that the coupling framework works for any DEM contact law, e.g., for elastic, plastic, cohesive, non-cohesive, wet and dry materials, because it is only the resulting forces/tractions that are mapped to the FEM mesh. Here, we use the same time step dt for the FEM and the DEM sub-model, although it can be much larger for the FEM model than for the DEM model.
The simplest test case that exploits the full capability of surface coupling is a particle sliding over an elastic cantilever. During this dynamic process, both the smoothness of particle-convex/concave wall interactions and the   surface coupling algorithm can be verified. Therefore, an elastic cantilever is implemented, using a 2 × 2 × 30 mesh of eight-node linear cubic elements, with an element size ∆X = ∆Y = ∆Z = 0.4 m. The left end of the cantilever is fixed, while the other boundaries are free to move and interact with the particles. To test whether the triangulation and the coupling scheme cause any asymmetry, two particles of diameter d = 0.5∆X = 0.2 m are positioned on the upper surface of the cantilever, at the center of the two leftmost elements (in the y direction) of the cantilever, as shown in Fig. 11. A linear-elastic contact model without friction is applied. The particles are subjected to gravitational forces while the elastic cantilever is not. The initial particle-beam overlap is chosen such that the particle-wall interaction force −f w αγ cancels the gravitational body force f b α . At time t = 0 s, the two particles are released with an initial velocity of v y = 0.1 m/s. As time proceeds, the particles move along the cantilever in the longitudinal (y) direction while bending it in the negative perpendicular (z) direction. Table 2 gives the CG widths selected for comparing the simulation results obtained using the conventional and CG-enriched formulations. In the limit of c = 0, Eq. (16) is used; for c > 0, we use Eq. (15), with the Lucy polynomial function Eq. (10) as the coarse-graining kernel. The CG width c is varied from 0-15 times the particle diameter d, or 0-7.5 times the finite element size ∆X , in order to identify the characteristic length scale that leads to a converged, reasonable macroscopic solution.

Role of coarse-graining in surface coupling
The benefit of a CG-enriched formulation for FEM-DEM surface coupling is illustrated by comparing the trajectories of the particles when they slide on and bend the elastic cantilever. A coarse-graining (CG) width of , respectively, show the trajectories of the two particles in the x-y and y-z plane. The analytical solution for the particle trajectory in the y-z plane is obtained by assuming a point load equal to the weight of the particles, exerted at various locations on the elastic cantilever, in static equilibrium. As the cantilever reacts to the load exerted by the particles with a certain amount of inertia, the mismatch between the analytical and numerical solutions increases, as shown in Fig. 12(a). Although both surface coupling formulations give numerical solutions in good agreement with the analytical solution, without CG enrichment (c = 0d), the particles oscillate quite significantly while bending the elastic cantilever. This numerical artifact is greatly reduced with CG and a CG width of c = 10d for the current example where the FEM mesh resolution is not high in comparison with the particle size.
Theoretically, both particles should only move along the longitudinal (y) axis with zero displacements along the lateral (x) axis. However, a point force contribution that is coupled only locally to one finite element may cause unexpected oscillations if there are jumps in the particle-wall interaction force, for example, when the particles move across the joint edges/vertices of the interface triangles. Fig. 12(a) shows the erroneous particle trajectories for c = 0. Because of the nonzero lateral velocity, the particles start to deviate from x = 0.2 and −0.2 m as soon as the cantilever starts to bend. The non-smooth particle velocity (e.g., in the x direction) eventually causes the particles to collide at y = 5.4 m. With a CG width c = 10d, not only is the symmetry of the particle trajectories in the x-y plane recovered, but also the oscillation in both particle trajectories in the z direction becomes much smaller, as shown in Fig. 12(b).

Dependence on the coarse-graining width
From the previous subsection, it is evident that with coarse-graining the surface coupling forces are mapped more accurately to the finite elements, and the particles' paths are more smooth as they slide along the surface of the cantilever. Here, the effects of varying the coarse-graining width are investigated. Fig. 13 shows the particle trajectories obtained with various CG widths c, including the c = 0d and c = 10d cases shown in Fig. 12. Note that the difference between the numerical and analytical solutions for the particle trajectories in the y-z plane z − z ref is plotted against y in Fig. 13(b), in order to better visualize the effect of CG width on the oscillatory behavior caused by the surface coupling forces. Fig. 13(a) shows that the asymmetry of the particle trajectories and the deviation from a constant x-component is reduced even for small CG widths (c = 4d). Furthermore, the oscillatory behavior of the particle trajectory in the z-direction is reduced, though not removed completely, as shown in Fig. 13(b). As c increases, both the lateral motion of the particles and the oscillatory behavior in the z direction decrease dramatically. It can be seen from Fig. 13(b) that, as the CG width increases, not only is the oscillation caused by the coupling reduced, but also the Fig. 13. Trajectories of the particles P 1 and P 2 sliding over the elastic cantilever and the discrepancy from the analytical solution. The coupling forces are projected from the particle-triangle interactions using a range of coarse-graining widths.
numerical result closer to the analytical solution. However, the discrepancy becomes larger as the particles accelerate due to gravity, and the cantilever gains inertia which violates the assumption for the analytical solution.
To obtain an objective measure of the numerical errors caused by the coupling methods, we compute the L 2 error norm, err (t) = ∥x sim (t) − x ref (t)∥ 2 , between the numerical and analytical solutions for the positions of the two particles. From the particle trajectories in Fig. 13, the evolution of the L 2 error while the particles slide on the cantilever can be calculated. The error evolution over time (data not shown here) is then averaged over time, i.e., err = ∫ T 0 err (t) dt/T , until the time T where the particles fall off the beam. The time-averaged error err is plotted for different coarse-graining widths in Fig. 15. The result clearly shows that the L 2 error is significantly reduced with an increasing coarse-graining width, from 0.12 m at c = 0d to less than 0.01 m at c = 10d. Increasing the CG width from c = 10d to c = 15d does not reduce the error anymore, and for even higher CG widths the FEM solver does not converge.
The decrease of the total error is almost quadratic and remains constant after c = 10d. This suggests that a CG width of around c = 10d is optimal, i.e., for this particular problem c = 10d is large enough to minimize discretization errors and small enough to preserve the physical gradients of the resulting fields. The suboptimal choice of c = 15d in comparison with c = 10d can be seen in Fig. 14, where the evolution of the L 2 error norm over the y coordinate of the particles is plotted. The two curves deviate when the particles move to y = 10 m away from the fixed end of the cantilever, resulting in slightly larger err for c = 15d, as shown in Fig. 15.

On the energy conservation properties of CG-enriched surface coupling
The previous results show that the particle motion oscillates around a trend that is close to the analytical solution, and the oscillations and numerical errors decrease with an increasing CG width. In order to investigate on the mechanisms that cause this numerical artifact, the motion of the elastic cantilever is analyzed. Here, we focus on the energetic aspects, including the kinetic, potential, and total energies of the cantilever and their differences that vary as a function of the CG width, and compare them to those of the particles/contacts. Fig. 16 shows the evolution of the kinetic energies E D E k and E F E k in the particle and continuum sub-models over time, for CG widths of c = 0d, 2d, 4d, 6d, and 10d. As seen above in Fig. 13(b), the oscillation is gradually reduced with an increasing CG width. From Fig. 16(a), it is clear that excessive kinetic energy is stored in the FEM sub-model and the amount of energy induced decreases with the CG width.
Finally, the time evolution of the total energy E t in the sub-models and their sum are plotted in Fig. 17. The potential energies of the particles and cantilever are chosen to be zero at the initial time step. Note that only the strain energy is used to calculate the potential energy, because gravity is not applied to the cantilever.
We observe that a large amount of induced energy remains mostly in the FEM sub-model and this effect becomes more significant as the CG width decreases (see Fig. 17(a)). The total energy in the DEM sub-model is almost    unaffected by coarse-graining, as shown in Fig. 17(b), except that the oscillations are significantly reduced. The total energy in each sub-model is conserved after the sub-models are uncoupled at t > 35 s, when the particles fall off the cantilever and keep accelerating under gravity.
Because both sub-models are purely elastic, the total energy of the coupled system should be strictly conserved irrespective of the coupling method. This is evidently not the case, as shown in Fig. 17(c), where the total energy of the coupled system in logarithmic scale is plotted over time. The initial total energy at t = 0 s arises from the nonzero particle-wall contact forces that maintain the force equilibrium in the z direction. Without coarse-graining, E t increases sharply only two seconds after the particles start to slide. As the CG width increases from 2d to 10d, the induced energy is decreased by several orders of magnitude. With c = 10d, the total energy is kept at the same order of magnitude between 0.9 and 2 J while the sub-models remain actively coupled, although there is a slight drop in E t after the sub-models uncouple at t ≈ 35 s. Even with c = 6d, the highest total energy that the coupled model has is less than 40 J, which is 125 times lower than the total energy obtained without coarse-graining. It seems that the optimal CG width is between 6d and 10d. However, the optimal CG width is not known before performing the coupled simulations and investigating on the energy conservation, and could also be problem-and material-dependent.

Volume coupling: wave propagation between discrete and continuum domains
In this section, the performance of the conventional and CG-enriched formulations for FEM-DEM volume coupling are compared against a full DEM reference case. The benefits of coarse-graining are exemplified by a simple numerical example, that is the propagation of elastic waves through a granular medium. This numerical example not only serves as a basis for verification (by comparing wave speeds in the granular material and the elastic medium), but also allows for a systematic investigation of the conservation properties of the coupling formulations.
Let us consider a beam of 20 m in length made from granular material. We model the first 12 m of the beam with DEM, and the last 12 m with FEM; thus, there is a 4 m long overlapping volume, as shown in Fig. 18(a). The goal is to propagate elastic waves between the discrete to the continuum domain, and investigate the influence of CG on propagation and conservation behavior. Data from a full DEM model with the same particle configuration and initial and boundary conditions are used to verify the volume coupling schemes with/without CG enrichment.
The FEM domain is discretized by a 2 × 2 × 30 mesh, using eight-node linear cubic elements with an element size ∆X = ∆Y = ∆Z = 0.4 m. The DEM particles are positioned in a cubic packing with zero overlaps and thereby zero initial interaction forces. The diameter of the particles is an integer fraction of the finite element size, d = ∆X/n, thus the number of particles per finite element is n 3 , with n = 1, 2, 3, or 4. Each pair of particles is linked with a linear spring k n , in both compression and tension directions, perpendicular to the contact plane. Both the DEM contact model and the FEM constitutive equations are linearly elastic, i.e. conserve energy. To match the (macroscopic) density and wave speed of the two submodels, the parameters of the micro-model are related to the macro-model parameters via ρ g = 6ρ/π , k n = E · d and ν = 0, due to the cubic structure (see [55]). The material parameters used for the sub-models are listed in Table 3.  Table 4 Initial conditions for the simulations in each subsection.

Sec. Waveform
Frequencies ω (kHz) All discrete particles and finite elements are free to move, except for those at the left and right boundaries. The right boundary is fixed, whereas the particles in the leftmost layer (perpendicular to the y axis) are given an initial prescribed velocity (0, v y , 0). The wave motions used in Sections 6.2.2-6.2.4 are summarized in Table 4. For t > t 0 , the boundary particles are either free to move (Section 6.2.2) or fixed in space (Sections 6.2.3-6.2.4) with zero velocity v y = 0. The highest frequency ( f = 1.95 kHz) was chosen such that a single wavefront contains approximately eight finite elements. For t > t 0 , the pulse wave propagates from the particles into the elastic beam until it is reflected from the fixed right end of the FEM domain.
The coupling domain Ω C between the DEM and the FEM sub-model is updated dynamically. As shown in Fig. 18(b), at a given time step, the coupling weight w applied to the DEM sub-model monotonically decreases from 1.0 on the left to 0.0 on the right, and the corresponding weights on the FEM side 1 − w are assigned to the FEM nodes. The profile of w in the coupling region follows a half-cycle cosine in the y direction and remains uniform in the x − z plane.

Setup of the test cases
In Sections 6.2.2 and 6.2.3, the CG-enriched volume coupling scheme is first verified by comparing the simulation results obtained using the two formulations described by Eqs. (18) and (19), with the full DEM case as reference. Details about the test cases for a range of CG widths, namely 0∆X -2.5∆X , are given in Table 5. The verification is done with the simplest configuration for the overlapping domain, that is n 3 particles per finite element in a lattice structure, see the case for n = 1 in Fig. 18. For Section 6.2.4, we increase the number of particles per finite element in order to investigate how CG-enriched volume coupling performs in "normal" conditions where the particles are smaller than the finite elements. Relevant parameters for the three configurations 2P, 3P, and 4P, are given in Table 6. The material parameters for the simulations reported in this section are summarized in Table 3. The penalty parameter ϵ, see Eq. (29), is chosen to be sufficiently large in order to satisfy the kinematic constraint as much as possible.

Role of coarse-graining in volume coupling
To verify the volume coupling formulation, we consider the simplest particle-element configuration, that is one particle per element in the coupling zone (see Fig. 18(b)). A sharp increase of the particle velocity, as defined by  Table 6 Penalty parameter for different finite particle-element configurations for volume coupling.  Fig. 19. Wave propagation from the granular sub-model to the continuum sub-model, using CG-enriched (c = 0.0∆X -2.5∆X ) and non CG-enriched coupling. Note, the coupling zone takes one third of each computational domain and is located in the middle from 6.67 to 13.33 m (see Fig. 18). The results from a full DEM model with the same geometry and initial and boundary conditions are added for verification.
the first row of Table 4, is given to the particles within the leftmost layer. From t > t 0 , the particles in the leftmost layer are not fixed in space. In this way, mass, momentum, and energy are strictly conserved; only momentum changes when the wave is reflected from the fixed boundary on the right. The data from a full DEM model with the same particle configuration and initial and boundary conditions are presented to verify the coupled FEM-DEM model. To understand the influence of CG on the conservation properties of the coupling scheme, the CG width c is increased from 0.5 to 2.5 times the finite element size ∆X, denoted as 0.5∆X -2.5∆X in Table 5. Figs. 19(a) and (b) show two snapshots of an elastic wave propagating through the coupling zone. The agreement between the full DEM and FEM-DEM data shows that volume coupling formulations, without and with CGH, have been correctly implemented and can effectively propagate an elastic wave between the particle and continuum models. Nevertheless, a slight deviation from the DEM waveform can be seen, after the wave passes through the coupling domain. The velocity profile given by the DEM sub-model shows an excellent agreement with that of the full model, irrespective of the CG widths chosen. This can be clearly seen in Fig. 19(a) as all markers for a given position overlap. However, as the wave enters the coupling zone, the waveform in the FEM domain starts to deviate from the DEM waveform, with a decreasing amplitude and a widening pulse width. This is likely due to numerical dissipation that occurs when the two systems are actively coupled and the fact that the elements have a larger length scale than the particles. However, with a CG width c = 1.5∆X the discrepancy between the waveforms in the FEM and the DEM domain is minimized. As the CG width increases beyond 1.5∆X , the decrease of the amplitude inside the coupling zone becomes larger. This is to be expected because of the spatial smoothing effect of CG. For this particular problem, it seems that the optimal CG width is 1.5∆X . After the wave travels through the coupling zone, the decay in amplitude seems to be slightly weaker for all CG-enriched cases. To further understand the effect of coarse-graining on the performance of the coupled model, the time evolution of mass, momentum, and energy in the respective sub-model and the totals as the weighted sums (see Eqs. (25) and (26)) are plotted in Fig. 20. From Fig. 20(a), for both formulations, irrespective of the CG width, the total mass is strictly conserved for all cases. As the wave enters the coupling zone (indicated by the deviation from the full DEM reference solution at t ≈ 0.03s), the y component of the linear momenta in the DEM sub-model starts to decrease while the counterpart in the FEM sub-model starts to increase, as shown in Fig. 20(b). Nevertheless, the total momentum remains nearly constant, as shown in Fig. 20(c).
The influence of CG and CG width is also visible in the time evolution of the total energy. An improvement on the energy conservation, compared to the non-CG case, when the elastic wave travels within the coupling zone can be seen in Fig. 20(d), although numerical dissipation is not avoided completely. Fig. 21 shows how the ratio of the dissipated to the total energy varies with the CG width. Compared with the non coarse-grained formulation, the dissipation is reduced from approximately 33.5% for c = 0∆X to 29.5% for c = 1.5∆X . A further increase of the CG width (e.g., c = 2.0∆X and 2.5∆X ) does not lead to a further decrease of the numerical dissipation. In fact, increasing c beyond 2.5∆X leads to numerical instability and a converged solution cannot be obtained. This is because when the particle velocity is homogenized over too many finite elements, the difference between u DE and u FE is excessively large, especially for those elements which do not overlap with the particles. As a result, the coupling terms that arise by penalizing u DE − u FE lead to non-converging iterations in finding the solution to approach u DE = u FE . As shown in Fig. 21, for a sharp signal, the CG width should not be larger than twice the element size.

From a sharp impulse to continuous waveforms
Although the previously considered boundary condition for the DEM model (i.e., free end on the left in Fig. 18(b)) ensures exact mass, momentum, and energy conservation, it is not convenient for prescribing arbitrary waveforms for the wave excitation. Differently from the previous section, the particles within the leftmost layer are fixed in space. These particles are given a time-varying velocity for t ≤ t 0 as defined by the waveforms described in Table 4. Note, this boundary condition has been widely used for the DEM modeling of wave propagation in dry and saturated granular media [56][57][58]. In the following, various waveforms will be utilized in the coupled FEM-DEM simulations, with the same particle-element configuration as in the previous section. For the sake of brevity, only the time evolution of the total linear momenta and total energies will be shown and discussed.
As shown in Figs. 22(a)-22(c), the benefit of CG can be observed at a higher input frequency (e.g., ω = 1.95 kHz). With CGH, the energy is better conserved as the CG width increases, similar to the trend shown in Fig. 20(d). Nevertheless, the optimal CG width seems to be between 1.5∆X and 2.0∆X; a further increase of the CG width (e.g., c = 2.5∆X ) leads to a slightly larger ratio of energy dissipation, compared with the ratio for c = 2.0∆X . It can seen from Fig. 22(d) that the numerical dissipation decreases with a lower input frequency (from 34% at ω = 1.95 kHz to 31% at ω = 0.49 kHz). Nevertheless, a decreasing numerical dissipation can still be observed as the CG width increases to 1.5∆X . One of the reasons for a decreasing energy dissipation ratio is the increasing wavelength that suppresses the particle-scale response, as ω decreases from 1.95 to 0.49 kHz. In the following sections, an input frequency of ω = 1.95 kHz is chosen because the decrease of numerical dissipation is more significant compared with other input frequencies, although the absolute values of ∆E t /E t are relatively large. With these input frequencies, the sub-continuum scale responses are ensured to be sufficiently large, i.e., at the length scales smaller than the finite element size.

Many particles in each volume-coupled finite element
In this section, we consider a more general situation for FEM-DEM volume coupling where multiple particles are coupled to each finite element. In such a case, the smallest wavelengths that can be propagated in the particle model are much smaller than those in the continuum model. This raises the question: can the sub-continuum scale responses, coming from the particles, be up-scaled into the continuum medium, and if so how much information is lost due to the coupling? In this section, we further investigate the performance of CG-enriched volume coupling by varying both the number of particles coupled to each finite element and the CG width. Four particle-element configurations are considered by reducing the particle diameter d such that each finite element in the coupling zone contains 1, 2 3 , 3 3 , or 4 3 discrete particles, respectively. Energy conservation under a single cycle of wave excitation. First, waveforms consist of a single pulse are considered. Here, we select and show only the simulation data obtained with an input frequency of ω = 1.95 kHz. Details on the waveforms can be found in Table 4. For the sake of brevity, only the weighted sum of the total energies in the sub-models are shown in Fig. 23. At ω = 1.95 kHz, a decreasing trend of energy dissipation with respect to the CG width is observed for each configuration, as shown in Fig. 23(b). The curves closely match with those in Figs. 21 and 22(d), although a minimum dissipation cannot be found at c = 1.5∆X for these waveforms.
Moreover, it appears that the number of particles per finite element also plays an important role in determining the energy dissipation. In the case of 4P (see Table 6), the percentage of dissipated energy is as high as 55%. With coarse-graining, the ratios ∆E t /E t can be reduced to 50%, 43%, and 35% for the particle-element configurations 4P, 3P, and 2P, respectively. One of the reasons for ∆E t /E t to increase is that with an increasing number of particles the non-affine motion becomes more considerable, and the non-affine motion that contributes to the fluctuating energy in the particle system is not incorporated into the kinematic constraint, without using CG.
Energy conservation under multiple cycles of wave excitation. Elastic waves generally contain many cycles and have more complex waveforms than a single-cycle cosine, and a wide range of frequency contents are usually present in the wave signals. Here, the waveforms are further generalized to include multiple cycles, as summarized in Table 4. The waveform utilized here remains the same as in Section 6.2.4 but three and five complete cycles are  considered, with an input frequency of ω = 1.95 kHz. The particle-element configurations are also kept the same as before, namely, each finite element is coupled with 2 3 , 3 3 , and 4 3 discrete particles, with the same boundary conditions. Again, only the weighted sum of the total energies will be discussed here. Figs. 24(a) and 25(a) show respectively the time evolution of the total energy in the coupled system under three and five cycles of wave excitation. As can be seen in both figures, the total energy in all cases monotonically decreases as the waves travel through the coupling zone. In the case of five-cycle wave excitation, the duration of wave propagation in the coupling domain is longer. An interesting observation from Figs. 24(a) and 25(a) is that the rate of numerical dissipation is lower when the wave travels in the coupling domain, resulting in decreased energy dissipation, as the number of repeated cycles increases. By comparing the curves in Figs. 24(b) and 25(b), one can see that the amount of energy dissipated under the three-cycle wave excitation exceeds that under the five-cycle wave excitation by 3%-8%, depending on the input frequency. The avoided numerical dissipation is largest for the case of 4P considering a five-cycle wave excitation, namely from 36% at c = 1.5∆X to 26% at c = 2.5∆X . By comparing the variation of numerical dissipation with Although the numerical dissipation is not removed via coarse-graining, a key message that can be interpreted from the previous results is that the effect of coarse-graining on FEM-DEM volume coupling is most significant when (1) many particles are present in each coupled finite element and (2) the mechanical load contains many cycles and large spatial and temporal gradients in the stress field (related to the high frequencies and small wavenumbers from the dynamic load). These two circumstances are usually met in the practical applications of FEM-DEM volume coupling, such as pile installation [26] and selective laser sintering [59] where non-homogeneous, dynamically varying strain fields are present.

Conclusion and outlooks
Based on a micro-macro transition from particles to continuum, known as "coarse-graining" [38], we have reformulated and generalized the surface coupling and volume coupling methods for the concurrent multi-scale modeling of granular materials. From momentum conservation, we derive micro-macro projection or homogenization rules that are shown to be the generalized forms of those reported in literature. For surface coupling, coarse-graining (CG) smoothly distributes the particle-element point interaction forces to several close-by integration points from their contact points. As a result, the surface coupling forces can be projected beyond the finite element on which the interaction is physically exerted. For volume coupling, we use CG to obtain a continuous velocity field from the discrete particles, then couple this field to the continuum FEM model, instead of coupling the particle velocities directly to the FEM displacement field These generalized coupling terms contain only one user-defined parameter, the coarse-graining smoothing width, setting a length scale w for the micromechanics-based macroscopic fields. To the authors' knowledge, this is the first work that utilizes CG to formulate diverse multi-scale (two-way) coupling methods, rather than simply using it to post-process particle-scale data (one-way). The limiting case, w = 0, yields the classical coupling formulations for both surface and volume coupling. A nonzero coarse-graining width allows the discrete particles and finite elements that do not overlap in space to be coupled dynamically.
In addition to the theoretical derivation of the coupling terms for surface and volume coupling, we discuss their implementation in the coupled MercuryDPM-oomph-lib code. Mainly written in C++, both codes use objectoriented programming techniques, including multiple inheritances, template classes, and function overloading. The flexibility of both codes allows implementing versatile wrappers around the existing finite element and discrete particle classes in oomph-lib and MercuryDPM. Therefore, the code framework is very general and offers the possibility to mix surface-and volume-coupling in one unified monolithic coupling framework. The surface and volume coupling methods are implemented in the MercuryDPM branch https://svn.mercurydpm.org/SourceCode/ Branches/OomphlibCouplingExplicit; see https://www.mercurydpm.org for installation instructions. We will make a future version of this code a feature of the main MercuryDPM repository in the near future.

CG-enriched surface coupling
The benefits of coarse-graining for the surface-coupling scheme are exemplified by modeling particle-cantilever interactions. The dynamic process of particles sliding over an elastic cantilever is simulated via surface coupling. Compared with the conventional approach that uses the Dirac delta function to map contact forces, the new, non-local, CG-enriched SCoupling formulation: (i) removes the nonphysical oscillations in the particle motion when the particles slide across and bend the elastic cantilever, (ii) recovers the symmetry of particle trajectories in the horizontal plane when two particles are released from both sides of the line of symmetry with an equal offset, (iii) gives a smaller deviation from the analytical solution and removes the excess energy, at an optimal CG width, that is 10d for the current numerical example.
Furthermore, in modeling the interactions between discrete particles and elastic bodies, the particle-convex/concave wall detection algorithm ensures that the interaction is smooth not only in magnitude but also in the contact normal direction, i.e., one unique contact force is defined.

CG-enriched volume coupling
The benefits of coarse-graining for the volume-coupling scheme are demonstrated by modeling an elastic wave propagating from discrete particles to a FEM-discretized continuum medium. In this example, an elastic beam is considered, with a coupling domain (as big as one-third of its volume) overlapping with the discrete particle domain. The particle-scale parameters are calibrated such that the homogenized macroscopic properties of the granular assembly match those of the elastic beam. A continuous field of coupling weight that monotonically increases from 0 to 1 is used for the FEM sub-model (1 to 0 for the DEM sub-model) within the coupling zone. Compared with the conventional approach that uses FEM basis functions to obtain a homogenized field, the new, non-local, CG-enriched VCoupling formulation (i) helps better satisfy the kinematic constraints in the coupling zone, which results in a somewhat better agreement between the coupled and the full DEM solutions, (ii) avoids computationally expensive global-to-local coordinate transformation, and (iii) improves energy conservation when waves propagate through the coupling zone, particularly in the case of complex waveforms with high-frequency/small-wavelength contents.
A huge benefit of the CG-enriched formulation is that nonphysical phenomena caused by the volume coupling, such as numerical reflection and dissipation, can be kept at a minimum. The new formulation also has the advantage to separate the sub-continuum scale response from the particles while enforcing the continuum-scale solutions to match with the coarse-grained fields.

Outlooks
The main contribution of this work is the generalization of FEM-DEM surface and volume coupling methods using coarse-graining, and the improved accuracy and reduced numerical energy loss/gain at optimal CG widths, in simple yet effective test cases. Through coarse-graining, continuum fields at various length scales (e.g., particle, mesostructure, continuum) can be obtained, and the resulting solution of the coupled model can be different from one and another. However, within a certain range of CG widths [37], the solution remains unique (see Fig. 15), indicating a characteristic length scale. It is only within this range that the information loss, e.g., due to the micromacro mapping or coupling, is kept at minimum. One could also think of the CG width as a numerical parameter that has to be optimized with respect to some constraint, such as an error tolerance, or conservation properties. One could also think of the CG width as a numerical parameter that has to be optimized with respect to some constraints, such as some error measures and conservation properties. The identification of suitable length scale (CG widths) has a practical relevance and is by itself an interesting topic, especially for general granular systems with strong polydispersity and strain localization. The interplay between particle sizes, the CG width, and the FEM discretization length scale in coupled particle-continuum systems should be better understood.
Another advantage of coarse-graining is the possibility to define higher-order terms such as stress and strain tensor fields. It would be interesting to investigate whether a surface/volume coupling formulation based on stress/strain fields would improve further the accuracy and performance of the multi-scale model. Therefore, ongoing work focuses on identifying appropriate coarse-graining parameters for various deformable materials or structures coupled with many particles.
One limitation of the current surface coupling implementation is that the finite element surfaces are discretized with rigid and flat triangles in DEM. This introduces certain roughness to the deformable bodies' surface and can be avoided by using mathematically well-defined geometric objects (e.g., curved walls [42]), mapped from the finite element. Other shape representations, including signed distance functions [60], level set methods [42,61,62], and Minkowski sums [63] could also be used. Secondly, currently, only the traction due to particle-wall interaction forces are mapped onto the finite elements; the kinetic traction due to velocity fluctuations is not explicitly considered. We speculate that the contribution of the kinetic traction might become important when the interaction becomes dynamic Another important extension of the presented coupling framework is to incorporate elasto-plasticity into the material's constitutive description of the continuum sub-model. For surface coupling, this would allow describing the strong impact of granular materials on structures, causing irrecoverable deformation in the structures. For volume coupling, implementing a continuum-level elasto-plastic constitutive law would effectively reduce the coupling domain size, leaving most inelastic material behavior (up to moderate particle-scale discontinuities) to the continuum sub-model and the rest to DEM. Examples of relevant physical processes include the sintering of polymer powders [59], the fragmentation of rocks and concretes [64], etc.
Currently, only fast, small-amplitude, transient loads are considered, such as wave propagation and particles rolling/sliding on a cantilever. Further work is needed to investigate the performance of the proposed coupling framework under steady-state mechanical loads, large deformation, and the phase transition of granular material, e.g., between solid and solid.

Declaration of competing interest
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

Code and data availability
We have shared the link to the source code hosted on a SVN repository (https://svn.mercurydpm.org/SourceCode/ Branches/OomphlibCouplingExplicit). The coupled code is openly accessible to the public. Part of the numerical data presented in Sections 6.1 and 6.2 are available at https://github.com/chyalexcheng/Oomph-Mercury.git. Nevertheless, all simulations can be reproduced by compiling and running the coupled code, using the parameters given in Tables 1-6. where λ = Eν (1+ν)(1−2ν) and µ = E 2(1+ν) are the bulk and shear moduli, with E and ν the Young's modulus and Poisson's ratio of the material. In oomph-lib, Hooke's law is generalized for finite strains, that is the 2nd Piola-Kirchhoff stress σ i j is linearly related to the Green's strain tensor γ i j = 1/2 (G i j − g i j ) via Here g i j and G i j are the metric tensors associated with the undeformed and deformed configurations; g i j = δ i j when there is no growth of the element (e.g., due to thermal expansion); G i j = X il ∂ψ l ∂ξ j is the derivative of global position with respect to Lagrangian coordinates ξ j where X il is the ith coordinate of the nodal position at node l. This constitutive law reduces to the classical version of Hooke's law for infinitesimal strains when G i j → g i j .
Appendix B. Derivation of volume coupling terms on the particle and the continuum sub-model The additional virtual work δC done by the kinematic constraint in the volume coupling domain Ω C is given by where u DE and u FE are the displacement fields coarse-grained from the particles and interpolated via the finite element basis functions ψ i (i = 1, . . . , N ), respectively, δx and δX denote the variation of position vector at the microscopic and macroscopic scales, and ε is the penalty parameter. Substituting u FE = ∑ N j=1 u FE j ψ j , u DE = ∑ N j=1 u DE j ψ j , δX = ∑ N i=1 δX i ψ i , and δx = ∑ N i=1 δx i ψ i , the discretized forms of δC become: where the homogenization operator Π (see Eq. (27)) is used to relate the virtual displacements from the finite elements to the particles via δx i = ∑ N p α=1 Π iα δx α .
Substituting Eqs. (37) and (39) into the virtual work equation 8 yields where the macroscopic and the microscopic variables are defined for the governing equations of FEM and DEM in Eqs. (25) and (26), and the coupling weights by Eqs. (22)- (24). Note that the additional coupling terms due to the kinematic constraint closely resemble that of a body force (density). Further, since the homogenization operators involve a normalization ∑ N p α=1 Π iα = 1, the coupling forces applied to the FEM nodes is equal and opposite to the coupling forces applied to the particles,

Appendix C. Supplementary materials
Supplementary material related to this article can be found online at https://doi.org/10.1016/j.cma.2022.115651.