Multi-scale modelling of granular materials: numerical framework and study on micro-structural features

A multi-scale model for the analysis of granular systems is proposed, which combines the principles of a coupled FEM–DEM approach with a novel servo-control methodology for the implementation of appropriate micro-scale boundary conditions. A mesh convergence study is performed, whereby the results of a quasi-static biaxial compression test are compared with those obtained by direct numerical simulations. The comparison demonstrates the capability of the multi-scale method to realistically capture the macro-scale response, even for macroscopic domains characterized by a relatively coarse mesh; this makes it possible to accurately analyse large-scale granular systems in a computationally efficient manner. The multi-scale framework is applied to study in a systematic manner the role of individual micro-structural characteristics on the effective macro-scale response. The effect of particle contact friction, particle rotation, and initial fabric anisotropy on the overall response is considered, as measured in terms of the evolution of the effective stress, the volumetric deformation, the average coordination number and the induced anisotropy. The trends observed are in accordance with notions from physics, and observations from experiments and other DEM simulations presented in the literature. Hence, it is concluded that the present framework provides an adequate tool for exploring the effect of micro-structural characteristics on the macroscopic response of large-scale granular structures.


Introduction
The intrinsic influence of the discrete micro-structure of granular materials on their effective material properties and structural response is nowadays well recognized. The morphology, material evolution and mechanical interactions at the particle scale all contribute to the observed macroscopic non-linear failure and deformation behaviour. Multi-scale approaches provide an ideal tool for the modelling of granular systems, as they allow to directly incorporate the complex behaviour of the discrete micro-structure into the response of large-scale structural problems. This is typically done by coupling the discrete element method (DEM), which accurately represents the complex particle behaviour at the micro scale [1][2][3][4][5][6][7][8][9][10], to the finite element method (FEM), B E. Bosco e.bosco@tue.nl 1 Department of the Built Environment, Eindhoven University of Technology, P.O. Box 513, 5600 MB Eindhoven, The Netherlands which enables to efficiently solve boundary value problems at the macro scale. As a general principle, each integration point in the macro-scale FEM model is connected to a corresponding DEM micro-scale model via the application of adequate homogenization relations. In specific, a macroscopic deformation measure is imposed on the granular micro-structure through the definition of appropriate boundary conditions [11,12]. The DEM model is solved in turn, providing the particle contact forces in the granular assembly. These forces are subsequently translated into a macroscopic stress measure, which is required to solve the boundary value problem at the structural level. Several examples of coupled FEM-DEM approaches for granular materials have been presented in the literature. In [13,14] a quasi-static multi-scale method is formulated within the framework of small deformations, whereby the role of the particle microstructure on the effective frictional failure response of macroscopic samples is analysed, with a special focus on the initiation of strain localization. In [15] a small-strain multi-scale framework is proposed that elegantly computes the mechanical response for various monotonic and cyclic loading prob-lems, whereby drained as well as undrained conditions are considered. In [16] this framework is applied for developing multi-scale insights into classical geomechanical problems, such as retaining wall and footing problems. Coupled FEM-DEM approaches are commonly validated by analysing the macroscopic structural response in experimental tests typical for granular media, such as a biaxial compression test [11,15,[17][18][19], a slope stability test [20], or a (cyclic) shear test [15].
In the current communication a novel multi-scale framework is presented for granular materials, which employs the formulation and implementation of the micro-scale boundary conditions recently published in [12]. This formulation is based on the first-order homogenization approach originally proposed in [11], which includes important aspects that are usually ignored in other homogenization methods for particle systems, namely (i) the Hill-Mandel micro-heterogeneity condition that enforces consistency of energy at the microand macro scales, (ii) the influence of particle rotations in the formulation of micro-to-macro scale transitions, and (iii) a rigorous generalization of the multi-scale relations within the theory of finite deformations. The implementation of the micro-scale boundary conditions is performed with a servocontrol algorithm that uses a feedback principle similar to that of algorithms applied in control theory of dynamic systems. The servo-control algorithm has several attractive features compared to other methods used for implementing microscale boundary conditions. Firstly, from the computational viewpoint the algorithm is relatively simple to implement. Secondly, it can be implemented at the level of the interface communicating information between the macro-scale FEM and micro-scale DEM models, whereby modifications of the FEM and DEM source codes are not needed. The algorithm can therefore be easily combined with commercial software, whose source codes generally are not available to the user. Thirdly, in contrast to the often-used penalty method, the servo-control methodology preserves the physical meaning of the homogenized stress measure derived from the granular assembly. Further, the limit case at which the micro-scale boundary conditions are met exactly is rigorously retrieved from the formulation, see [12] for more details.
The first aim of this communication is to demonstrate how the servo-control algorithm for the micro-scale boundary conditions can be conveniently incorporated in a multiscale FEM-DEM framework. Accordingly, the governing equations of the multi-scale framework are formulated, and their numerical implementation is validated by comparing the computational results obtained for a quasi-static biaxial compression test to those calculated by direct numerical simulations. The convergence behaviour of the numerical results under mesh refinement is analysed, and the heterogeneity of the mechanical response across the specimen height is explored. The second aim of this communication is to show how the FEM-DEM framework can be used for analysing the influence of micro-structural characteristics on the macroscopic response of a granular system. Using the biaxial compression test, the microscopic properties selected for the variation study are the particle contact friction, the particle rotation and the initial fabric anisotropy. The influence of these properties on the overall, macroscopic response is analysed by computing the evolution of the effective stress, the volumetric deformation, the average coordination number, and the induced fabric anisotropy. This study is essential for gaining confidence in the quality of the multi-scale formulation; nonetheless, most other works on coupled FEM-DEM modelling do not consider such a study, but refer to a specific example simulation for the validation of the proposed method.
The paper is organized as follows. Section 2 presents the numerical homogenization framework for particle aggregates by defining the macro-scale and micro-scale models and the scale transition relations. Section 3 discusses numerical implementation aspects. The explicit time integration scheme adopted for the macro-scale problem is outlined, and details are provided on the dynamic relaxation procedure applied for satisfying the equilibrium conditions, and on the servo-control algorithm used for defining the boundary conditions at the micro scale. The section ends with a presentation of the coupled FEM-DEM solution algorithm. In Sect. 4, the performance of the proposed multi-scale framework is analysed for a biaxial compression test by comparing the computational results to those obtained by direct numerical simulations. A mesh convergence study is performed, and the role of several micro-structural parameters on the macroscopic response is studied. Some concluding remarks are provided in Sect. 5. In this communication the following notation will be used. The cross product and dyadic product of two vectors are denoted as a × b = e i jk a i b j e k and a ⊗ b = a i b j e i ⊗ e j , respectively. Here e i jk is the permutation symbol, e i , e j and e k are unit vectors in a Cartesian vector basis, and Einstein's summation convention is used on repeated tensor indices. The inner products between two vectors and two secondorder tensors are given by a · b = a i b i and A : B = A i j B i j , respectively. The action of a second-order tensor on a vector is indicated as A · b = A i j b j e i . The symbol ∇ indicates the gradient operator with respect to the reference configuration, and |.| refers to the absolute value of a variable. Occasionally, field variables referring to the macroscopic scale are indicated by an overbar, for instanceF, in order to avoid misinterpretation.
The present study focuses on two-dimensional particle aggregates. Accordingly, the dimensions related to volume, area, stress and mass density are consistently presented in their reduced form as length 2 , length, force/length and mass/length 2 , respectively.

Multi-scale framework for particle aggregates
This section treats the main principles of a multi-scale homogenization strategy for granular structures. These principles ensue from transforming relevant theorems used in classical first-order homogenization theories [21][22][23][24] from a continuous setting to a discrete setting. For more details on this aspect the reader is referred to [12,25,26].

Macro-scale problem
Consider a two-dimensional macroscopic domain with an initial, undeformed volume and boundary ∂ , characterized by a heterogeneous, granular micro-structure. The macroscopic domain is subjected to loadings and constraints at the boundaries under which the separation of scales principle holds, i.e., the characteristic length scale of the micro-structure is much smaller than the typical length across which the macroscopic deformation varies. Under this assumption, the macroscopic domain may be considered as a Boltzmann continuum, governed by the classical equilibrium equations: with the boundary conditions In Eq. (1),P is the first Piola-Kirchhoff stress,b is the body force per unit volume, andF = ∇(x) is the macroscopic deformation gradient, which is a function of the current position x. In relations (2) and (3) the external tractiont * is imposed on the boundary ∂ t characterized by the normal direction N, and the displacementū * is prescribed on the boundary ∂ u .
In order to solve the boundary value problem defined by Eqs. (1)-(3), a constitutive relation between the stress and the deformation is required. Instead of assuming a phenomenological constitutive equation, a multi-scale procedure is adopted that retrieves the constitutive response numerically from a computational analysis of the granular domain at the micro scale. The main features of the multi-scale scheme are illustrated in Fig. 1. At each material point the macroscopic deformation gradientF is calculated, and subsequently imposed on the corresponding micro-structural domain via appropriate micro-scale boundary conditions. After solving the response of the granular medium at the micro scale, from the particle contact forces at the boundaries the effective macroscopic stressP of the particle medium is computed, which is returned to the macro scale to solve the macroscopic equilibrium expressed by Eq. (1).

Micro-scale problem
The micro-scale geometry is represented by a twodimensional square domain composed of P + Q rigid particles. These particles are partitioned into P inner particles P p with p = 1, . . . , P, defining the initial interior domain V , and Q boundary particles P q with q = 1, . . . , Q, defining the undeformed boundary ∂ V . The boundary particles can be further split into corner particles P c with c = 1, . . . , 4 and the remaining edge particles P e with e = 1, . . . , E = Q − 4. The reference configuration of the centroids of the inner and boundary particles is denoted by the position vectors X p ∈ P p and X q ∈ P q , respectively. The undeformed microscale domain is schematically shown in Fig. 2a. Fig. 2 a Two-dimensional particle aggregate of undeformed volume V and boundary ∂ V . Light blue and blue colors refer to inner P p and boundary P q particles, respectively; b particle contact forces f c p acting on inner particle p ∈ P p in its current position x p ; c boundary forces a q , boundary moments m q , and particle contact forces f c q acting on boundary particle q ∈ P q in its current position x q . (Color figure online) The effective response of the granular assembly is derived by using classical homogenization principles, in the transition from a continuous to a discrete description. In this perspective, the finite area vector A q and the boundary forces a q are defined at the centroids of the boundary particles P q as Here, N is the normal vector associated to the undeformed boundary ∂ V of the micro-scale domain, t is the boundary traction in the reference configuration, and ds denotes an infinitesimal part of the boundary surface. The initial area vector A q is computed by accounting for the different radii of the boundary particles [12,25]: in which R q+1 , R q and R q−1 are the radii of adjacent boundary particles q + 1, q and q − 1, respectively. In addition, e 3 is the unit vector in the out-of-plane direction of the twodimensional particle structure, as indicated in Fig. 2. It is emphasized that in Eq. (5) the boundary particles must be numbered in the counterclockwise direction to arrive at an area vector pointing in the outward normal direction of the boundary. The kinematics of a rigid particle i within the granular assembly departs from the linearisation of the macroscopic deformation map. The current position x i of the centroid of particle i can be expressed as The first term in the right hand side of Eq. (6) reflects the contribution on the micro-scale kinematics by the macroscopic homogeneous deformation gradientF, and the second term indicates the local fluctuation w i of the micro-scale position field with respect to the applied homogeneous deformation. The macroscopic deformation is imposed via the frame of boundary particles P q , which induces contact forces f p and f q on the inner and boundary particles, and boundary forces a q and moments m q on the boundary particles, see Fig. 2b, c. The force and moment equilibria of the overall granular micro-structure can be expressed as Q q=1 a q = 0 and Q q=1 (x q × a q + m q ) = 0 for q = 1 . . . , Q, (7) with x q the current position vector of the boundary particles. Additionally, local equilibrium conditions must be formulated for each of the inner particles P p and for each of the boundary particles P q : Equation (8) describes force and moment equilibrium of N c p contact forces f c p at discrete contact points x c p on the surface of the interior particle p, with respect to its current configuration x p , see Fig. 2b. Analogously, relation (9) expresses force and moment equilibrium of N c q contact forces f c q at contact points x c q on the surface of the boundary particle q, in relation to its current configuration x q , see Fig. 2c. Note that the combination of expressions (8) and (9) results in expression (7).
In order to solve the micro-scale problem, a contact law describing the particle interactions is finally required. In the present work, a stick-slip contact law is adopted that relates the contact forces f c i to the corresponding contact displacements u c i as [1] f n = k n u n and f s = where the superscript c and subscript i on the contact force and contact displacement have been dropped for the sake of clarity. The normal particle contact force f n (with tension considered as positive and compression as negative) depends on the normal particle overlap u n between two particles in contact through a multiplication by the normal contact stiffness k n . When the contact is fully sticking, the tangential particle contact force f s is proportional to the relative tangential displacement u s at the particle contact via a multiplication by the tangential contact stiffness k s . This elastic constitutive relation holds up to a limit value at which frictional sliding is initiated, as defined by the normal force multiplied by the particle contact friction coefficient μ.

Macro-to-micro: kinematics and boundary conditions
The kinematical averaging relation is an essential ingredient for establishing the micro-to-macro coupling, by requiring the macro-scale deformation gradient to be equal to the volume average of the local, micro-scale deformation gradients. In a discrete setting, this is equivalent to the expression [11,12,26] Equation (11) is enforced to obtain appropriate boundary conditions in terms of displacements and rotations of the boundary particles of the granular micro-structure. Different types of boundary conditions may be selected to satisfy the constraint given by Eq. (11), see e.g. [11,12], among which the periodic boundary conditions adopted in the present study. Periodic boundary conditions already provide a realistic effective response for micro-structural volumes of relatively small to moderate size, which is commonly bounded by the upper and lower estimates obtained from, respectively, displacement and traction boundary conditions [12,22,23]. Under periodic boundary conditions both the displacements and rotations of the boundary particles P q are subjected to periodicity requirements: with θ q the magnitude of the boundary particle center rotation θ q = θ q e 3 , where e 3 is the unit vector in the out-of-plane direction of the 2D particle structure. Superscripts + and − refer to corresponding particles on opposite boundaries + and − of the granular assembly. Note that Eq. (12) 1 directly follows from Eq. (6), by using the periodicity requirement that the local fluctuations at two opposite boundaries must be equal, w + q = w − q . From the viewpoint of equilibrium, the forces and moments on opposite periodic boundaries need to be anti-periodic, thus satisfying the relations with m q the magnitude of the boundary moment, i.e., m q = m q e 3 .

Micro-to-macro: macroscopic stress and Hill-Mandel condition
In the micro-to-macro scale transition the macro-scale first Piola-Kirchhoff stress tensor is required to be equal to the surface average of the micro-scale forces a q acting on the boundary ∂ V of the particle aggregate: This expression, together with Eqs. (11) and (12), satisfies the condition on energy consistency between the macro and micro scales, known as the Hill-Mandel micro-heterogeneity condition [27]. For a discrete particle system, the Hill-Mandel condition specifies into [12,25] Equation (15) essentially states that the volume average of the virtual work applied at the boundaries of the granular micro-structure equals the virtual work of a macroscopic material point.
Finally, the numerical results of the multi-scale simulations will be presented in terms of components of the macro-scale Cauchy stress tensorσ . This stress measure is obtained from the first Piola-Kirchhoff stressP computed through (14) by using the common transformation rule: 3 Numerical implementation

Finite element formulation
In order to determine the solution of the macro-scale problem, a finite element formulation based on the theory of large deformations is employed. To this end, a total Lagrange scheme is adopted, for which the initial, undeformed macroscale domain is discretized into the domain h by using n e finite elements of volume e ∈ h . In the finite element formulation the strong form of the quasi-static macroscopic equilibrium (1) is transformed into the weak form, which, after integration by parts and using Gauss theorem, results in The first and second integrals in Eq. (17) represent the internal work and external work contributions, computed with respect to the test functions δF = ∇(δū) and δū, respectively. In the spatial discretization procedure, the weak form (17) is approximated by formulating the continuous displacement fieldū in terms of finite element interpolation functions, which leads to a system of non-linear algebraic equations Here, the vector U contains the nodal values of the macroscopic displacement fieldū. The internal and external force vectors are given by where N and B are matrices incorporating the interpolation functions and their spatial derivatives, respectively.

Dynamic relaxation
In the present work an explicit time marching scheme based on dynamic relaxation is adopted. The purpose of the dynamic relaxation method is to reach static equilibrium from the equations of motion in a relatively fast and numerically robust fashion, by effectively dissipating the kinetic energy of the modelled system. This requires the computation of the effective macroscopic stress tensorP, calculated from expression (14) via the boundary forces acting on the granular micro-structure, but circumvents the additional computation of the (computationally expensive) constitutive tangent matrix typically required in implicit time marching schemes. Correspondingly, the macro-scale balance equation, originally given by relation (18), takes the form whereM is the mass matrix andC is the damping matrix, which here is taken proportional to the mass matrix as C = αM, with α a viscous damping coefficient. Further,Ü n , U n and U n designate the nodal acceleration, nodal velocity and nodal displacement vectors at time t n , respectively. In correspondence with [28], the nodal accelerationsÜ n at time t n are approximated through a central difference scheme Inserting relations (21) into the equation of motion (20) results in the update for the nodal velocitiesU at time t n+ 1 2 : Subsequently, the displacement vector U n+1 is computed from the velocity vectorU n+ 1 2 as In the above relations a lumped (diagonal) mass matrix is used, for which the diagonal terms at nodes k = 1, . . . K of element e follow from [29] where K is the total number of nodes of element e, N k is the shape function referring to node k, andρ is the macroscopic density computed as the product between the particle density ρ and the packing volume fraction v m , i.e.,ρ = ρv m . The damping coefficient α appearing in Eq. (22) is adjusted in each iteration j of time step n as [30] α = α n, j = 2ξω n, j , with ξ the damping ratio and ω n, j a frequency parameter, computed as [31]: Finally, for warranting the stability of the solution, the time increment t is iteratively updated as [30] with γ a safety factor. This safety factor ensures that the time increment t can be prescribed sufficiently small to avoid numerical divergence characteristic of an explicit timemarching scheme. The solution of Eq. (20) is considered to be converged when the ratio between the systems' kinetic energyĒ n, j kin in the current iteration j of time step n and the maximal kinetic energyĒ n kin,max reached during the time step is less or equal than a prescribed tolerance tolĒ : with the kinetic energy computed as The macro-scale solution procedure has been implemented in the finite element code ESyS-Escript [32]; details on the solution algorithm are provided in Sect. 3.3. The package mpi4py 1 was used for parallelizing the DEM computations at the different integration points.

Dynamic relaxation
The micro-scale boundary value problem, consisting of the equilibrium Eqs. (7)- (9), the constitutive response of the particles (10), and the boundary conditions (12), is solved by applying a dynamic relaxation method, similar to the approach adopted for the macro-scale problem. For each particle i, with i = 1, . . . , P + Q, the generalized equation of motion can be expressed as where the mass matrix M i = diag [M i , I i ] includes the particle mass M i and the mass moment of inertia I i = M i R 2 i /2, with R i the particle radius. The termd i represents the generalized acceleration vector, calculated as the second time derivative of the generalized coordinate vector d i = [x i , θ i · e 3 ] T . The generalized coordinate vector contains the current locations of the particle centres x i and the particle rotations θ i . The vector p r = [f r , m r · e 3 ] T is the generalized force vector composed of the resultant force f r and moment m r acting on particle i. Analogously, 3 ] T is the generalized vector containing the resulting particle force and moment following from the artificial dissipation applied in the simulations. Based on [33], the artificial dissipative force f d and moment m d are defined as f d = −α|f r | sign(ẋ i ) and m d = −β|m r | sign(θ i ). The symbols α and β are damping values that are coupled to (signum) functions of the particle translational velocitẏ x i and rotational velocityθ i , respectively. The equations of motion (30) are integrated using an explicit, first-order finite difference scheme, the details of which can be found in [34]. The dynamic relaxation process is considered to be converged when the ratio between the kinetic energy E kin of the inner particles P p in the granular medium and their potential energy E pot is less or equal than a prescribed tolerance [35], i.e., with where f c n /k n and f c s /k s reflect the relative elastic displacements in the normal and tangential directions of particle contact c and N c is the total number of particle contacts.

Servo-control algorithm for micro-scale boundary conditions
The periodic displacement and periodic rotation boundary conditions discussed in Sect. 2.3.1 were implemented by means of a servo-control algorithm [12]. This algorithm is based on finding an iterative correction for the boundary particle displacements and rotations in order to reduce the difference between the measured and the required values of the boundary condition. Considering the anti-periodicity conditions (13) required for boundary forces and boundary moments, the corresponding residuals for the edge particles are a e = a + e + a − e , m e = (m + e + m − e ) · e 3 for e = 1 . . . , E/2.
Introducing the gain parameters g a and g m , the corrections for displacements and rotations of the edge particles can be obtained by multiplying the residuals by the corresponding gains: For the four corner particles no displacement correction is needed as their displacements are directly imposed as a function of the macroscopic deformation, i.e. x c =FX c . The rotations of the corner particles are updated in a similar way as done for the edge particles in relation (34), in accordance with the correction with This correction is added to the displacement and rotation calculated at the previously converged increment. The process is considered to be converged when the residual in terms of the particle forces and particle moments, defined in expressions (33) and (36), satisfies the criterion The residual r a is characterized by the normalized L 2 -norm of the incremental boundary forces a e and moments m e and m c of the edge and corner particles, see Table 2, and a is a pre-defined tolerance. The implementation of the micro-scale problem, involving the dynamic relaxation procedure and the servo-control algorithm for the boundary conditions, has been performed by using the open-source discrete element code ESyS-Particle [36,37]. The numerical algorithm for the multi-scale framework is discussed in detail in Sect. 3.3.

Multi-scale FEM-DEM coupling
In order to perform a multi-scale analysis of a granular system, the macroscopic continuum formulation treated in Sect. 3.1 is coupled to the discrete micro-scale model described in Sect. 3.2. The coupled macro-micro solution algorithm is summarized in Table 1.
The macroscopic domain is discretized in n e finite elements, with n i p integration points (IPs) per element. Before the loading is applied, an identical granular micro-structure is assigned to each integration point. The multi-scale coupling is realized by following a deformation driven procedure, consistent with the homogenization strategy presented in Sect. 2. As pointed out in Table 1, the external load is applied to the macroscopic domain in an incremental fashion, whereby for each macroscopic integration point the macroscopic deformation gradientF is incrementally updated from the nodal displacements, and subsequently imposed upon the frame of boundary particles of the micro-scale granular packing. At the onset of each increment the macroscopic deformation is applied homogeneously to all the boundary particles P q , and the particle assembly is dynamically relaxed to its equilibrium state. Subsequently, an iterative loop is entered, whereby the corrections for the displacement and rotation of the boundary particles are calculated following the servo-control methodology, see relations (33)- (36). After the application of these corrections, the particle system is again dynamically relaxed to the equilibrium state, and the residual r a is computed as a function of the current values of the boundary forces and moments. When the convergence criterion given by (37) is satisfied, the iterative loop is terminated. The macro-scale stressP is calculated from the micro-structural boundary forces in accordance with expression (14), and its value is transferred back to the corresponding integration point in the macroscopic domain. This procedure is performed for all integration points in the macroscopic domain. Table 2 summarizes the solution algorithm for the micro-scale DEM simulation.
Once the macroscopic stress is computed in the integration points, the macroscopic equations of motion (20) are solved using the iteratively adjusted values of the damping coefficient (25) and time step (27). The nodal velocities and nodal displacements are calculated from relations (22) and (23). When criterion (28) is satisfied, the current increment is considered to be converged. When this is not the case, the procedure above is repeated, until expression (28) holds.

Computational results
In this section the proposed FEM-DEM multi-scale framework is validated on a series of representative numerical simulations. A reference problem is defined first, for which a mesh convergence study is performed to establish the appropriate element size for the FEM model. Subsequently, the influence of various micro-structural properties on the macroscale response is investigated.

Definition of the reference problem
The macro-scale domain consists of a rectangular specimen of dimensions 10 mm × 20 mm, supported vertically at the complete bottom edge and horizontally in the lower left corner node. The domain is discretized into n e bilinear quadrilateral elements, with n i p = 4 integration points per element. The specimen is first subjected to isotropic compression with the stress magnitudeσ 0 = 0.15 MN/m applied in ten loading steps, see Fig. 3a. Next, a biaxial compression loading stage is initiated, whereby the vertical displacementd is increased incrementally up to a vertical strain ofε =d/h ic = 10% of the sample height h ic obtained after isotropic compression, see Fig. 3b. The contribution by the gravitational loading to the sample response is relatively small, and therefore may be ignored. Four different FEM discretizations of the macro-scale domain are considered, as  Table 2 • Prescribe periodic boundary conditions  (27) • Compute nodal velocities using (22) • Compute nodal displacements with (23)

5.
Check for convergence (28) in terms of kinetic energy • If converged ⇒ go to next increment 2 • If not converged ⇒ go to next iteration 3

Initialize boundary conditions by applying updated macro-scale deformation homogeneously
x q =FX q and θ q = 0 for q = 1 . . . , Q

Dynamic relaxation until convergence criterion (31) is satisfied.
Obtain boundary forces and moments.

Update particle configuration
Partition the boundary into corner c and edge e particles Calculate edge particles displacement u e and rotation θ e corrections via (33)- (34) and corner particles rotation corrections θ c via (35)-(36)

Dynamic relaxation until convergence criterion (31) is satisfied.
Obtain boundary forces and moments.

Calculate residual
with M k and R k the mass and radius of particle k; k ∈ {c, e}. 6. Check for convergence: r a ≤ a 6A. if converged ⇒ Save current configuration, compute macroscopic stressP with (14) and go to macro-scale simulation.  The constitutive behaviour in each macroscopic integration point follows from the effective response of an initially (almost) isotropic polydisperse packing composed of 228 particles, see Fig. 3c, d. The chosen number of particles is based on the convergence study performed in [12], which shows that under simple shear deformation and periodic boundary conditions the effective stress response for this number of particles does not significantly change under a further increase of the sample size. The initial packing structure is generated by a collision-driven molecular dynamics code and subsequently reconstructed into a geometrically periodic packing, see [12] for more details. The particle radii are taken from a uniform size distribution with polydispersity R max /R min = 1.5, with the minimal particle radius in accordance with R min /L = 0.03, where L is the side length of the initially square micro-structural domain. Note that from dimensional considerations it follows that the sample response is uniquely determined via the specific values chosen for the ratios R min /L and R max /R min . The volume fraction of the initial packing is v m = 0.846, and the initial coordination number (i.e, before the application of isotropic compression) isn 0 = 2C/(P + Q) = 3.42, with C the total number of particle contacts and P + Q the total num-2 Occasionally, the damping ratio ξ and safety factor γ were modified during the simulation to improve the convergence speed. ber of (inner + boundary) particles. The fabric anisotropy is measured as where 1 and 2 are the eigenvalues of the fabric tensor with N c p the total number of particles in contact with particle p and n pc the unit vector pointing from the centroid of particle p to the specific contact point c. The initial anisotropy of the reference packing is A 0 = 0.02. The normal and tangential stiffnesses in the stick-slip particle contact model (10) are k n = 10 4 N/m and k s = 2 · 10 3 N/m, representing relatively soft particles, and the friction coefficient equals μ = 0.4. The density of the particles is ρ = 2·10 3 kg/m 2 . The translational and rotational damping factors used in the dynamic relaxation procedure are α = β = 0.7, and the time increment is t = 10 −6 s. The gain parameters adopted in Eq. (34) for the application of the boundary conditions are (in dimensionless form) g a M/ t 2 = 1 · 10 2 and g m M i R 2 i / t 2 = 2 · 10 2 , with M i = ρπ R 2 i representing the mass of particle i and R i its radius. The tolerances in expressions (31) and (37) are taken as tol E = 10 −3 and a = 10 −4 , respectively. The parameters used at the macroscopic and microscopic levels of the multi-scale framework are summarized in Table 3.

Mesh convergence study
In order to explore the mesh sensitivity of the multi-scale approach, the macroscopic domain depicted in Fig. 3a has been discretized into four different meshes, which are characterized by the following number of elements: n e = [1 × 1, 1 × 2, 2 × 4, 4 × 8], with the corresponding number of integration points as n i p = [4,8,32,128], respectively. The specimen first undergoes isotropic compression, followed by biaxial compression, as described in Sect. 4.1.
The overall stress-deformation response for the considered meshes is shown in Fig. 4. The stress is measured in terms of the stress ratioσ , given bȳ which is taken as the average over all integration points, whereσ 11 andσ 22 are the normal components of the macroscopic Cauchy stress (16) in the axial and lateral directions of the specimen, respectively. Note that the biaxial stress ratiō σ given by Eq. (40) directly reflects the macroscopic friction angle φ m via the usual relation φ m = arcsin(σ ). The linear axial strainε is imposed on the specimen via the vertical displacementd: where h ic is the sample height obtained after the initial, isotropic compression stage.
The results of the multi-scale simulations are compared to those of two direct numerical simulations (DNS), in which the rectangular specimen is described by DEM models of 441 and 1681 particles. These models are referred to as DNS A and DNS B, respectively, and the corresponding samples are constructed by copying the square L × L micro-structural domain (used in the multi-scale simulations) multiple times along the sample height and width, see Fig. 7e, f. Accordingly, the length scale L defines the sample height and width as 2L × L for DNS A and 4L × 2L for DNS B. The particle radius is kept the same in the two DEM models, in accordance with the ratios R max /R min = 1.5 and R min /L = 0.03 adopted for the micro-structure used in the multi-scale simulations. The traction boundary conditions characterizing the initial isotropic compression stage and the displacement boundary condition defining the biaxial compression stage are applied in accordance with the servo-control procedure presented in [12]. Figure 4 illustrates that the average stress response of the coupled FEM-DEM models for all meshes considered is very close to the predictions of the direct numerical simulations, in particular in the pre-peak regime, 0 ≤ε ≤ 0.05. After this point, a moderate softening behaviour is observed, whereby the responses for the different meshes start to deviate from one other. The mesh size dependency of the FEM response in the post-peak regime is a well-known effect; to circumvent this problem in a multi-scale setting, kinematically enriched multi-scale frameworks have been proposed in the literature, see [38,39]. Since in the present work the focus is mainly on the pre-peak regime of the macroscopic response, the appli- cation of these frameworks for granular systems is considered as a topic for future research. The influence of the choice of the macroscopic mesh size on the effective response is further investigated by considering the evolution of the coordination numbern and the induced fabric anisotropy A (both averaged over all macroscopic integration points), see Fig. 5a, b. As a general trend, it can be observed that the results are only slightly sensitive to the mesh adopted. The coordination number somewhat decreases with increasing deformation, due to the horizontal expansion of the macroscopic domain. The anisotropy initially increases, since the packing deforms stronger in the vertical direction than in the horizontal direction, in correspondence with the macroscopic biaxial loading conditions applied. The induced anisotropy becomes maximal at about the same deformation stage as at which the peak strength is reached. In the softening regime, a small decrease in anisotropy is observed, which is caused by a dilating particle structure that develops under progressive shear failure.
In addition to examining the development of the average stress in the macroscopic domain, the variation in stress is considered by plotting the stress evolutions in the individual integration points of the FEM model. Figure 6a, b show the stress ratioσ in all the integration points of two specific meshes selected, which have n e = [2 × 4] and n e = [4 × 8] elements, respectively. The thick black line represents the corresponding average stress response taken from Fig. 4. In general, the spread in stress values in the relatively coarse mesh is smaller than in the fine mesh, which illustrates that the fine mesh more accurately describes the heterogeneous response of the macroscopic domain. The heterogeneity of the response clearly becomes rather strong in the softening regime, due to a localization of the macroscopic deformation pattern. This effect can be further explored by depicting the micro-structures at a vertical deformation ε = 0.09 in two different integration points, corresponding to the locations with the highest (bold dark grey line) and lowest (bold light grey line) stress levels at the macro-scale. The micro-structural responses clearly show differences in both the overall deformed shape and in the force chains developing within the particle structure (indicated by the red lines).
In conclusion, the results in Figs. 4 and 5 demonstrate that even the coarsest mesh with n e = [1 × 1] accurately captures the average multi-scale stress-deformation response. This supports the capability of the proposed multi-scale framework to adequately perform large-scale simulations with a relatively coarse mesh discretization of the macroscopic domain at limited computational cost. Obviously, a finer mesh allows for a richer description of the heterogeneous macroscopic response. This aspect, which has been noticed from the stress evolutions in Fig. 6, is further elaborated by comparing the deformed macroscopic configurations computed by the multi-scale method and the DNS at the end of the deformation processε = 0.1, see Fig. 7. Here, the variable used in the contour plots is the incremental equivalent strain ε eq , defined as the L 2 -norm of the increments of the components of the linear strain tensorε =F − I. It can be observed that the deformation of the macroscopic domain computed by the multi-scale simulations lies closer to those computed by DNS when the mesh becomes finer. In addition, the localization pattern becomes more prominent for a finer mesh. Based on these results, it is concluded that a good compromise between the accuracy in the description of the local macroscopic features and the computational time is found for the mesh of n e = [2 × 4] elements with n i p = 32 integration points. Hence, this mesh will be used for the forthcoming multi-scale simulations in which the influence of micro-structural parameters on the effective macroscopic response is analysed. As a final note, from the deformed configurations computed with DNS it can be observed that a few inner particles have been pushed through the frame of boundary particles during the dynamic relaxation procedure, thereby ending up outside the actual macroscopic domain. A way to prevent this from happening is by extending the thickness of the boundary frame with more particles. This solution, however, has not been implemented here, since the influence of this effect on the macroscopic response is negligible.

Influence of micro-structural parameters on the macroscopic response
The macroscopic geometry in Fig. 3 is discretized with the selected mesh of 2 × 4 finite elements, whereby the particle contact friction, the particle rotation, and the initial fabric anisotropy are varied. In the analysis of the results the attention will be focused on the pre-peak regime of the macroscopic response, during which the FEM results are independent of the mesh size, see Sect. 4.2

Particle contact friction
The effect of particle friction on the macro-scale response is investigated by considering three different values for the particle contact friction, namely μ = [0.2, 0.4, 0.6]. Note that μ = 0.4 is the friction coefficient of the reference particle packing studied in Sect. 4.2 The macroscopic stress evolution, expressed in terms of the stress ratioσ averaged over all the integration points, is presented in Fig. 8a. In addition, the evolution of the volumetric deformation Fig. 8b, whereF ic is the deformation gradient, evaluated at the end of the preliminary isotropic compression stage. As expected, from Fig. 8a it follows that in the pre-peak regime the stress ratio clearly is higher for a larger particle contact friction. In addition, the development of compaction (corresponding to a negative value ofJ rel ) during the initial stage of deformation, as typical for packings with relatively soft particle contacts, is larger for a higher value of the particle contact friction. Essentially, the normal contact forces become larger when the particle contact friction increases, which generates more particle overlap, and thus more compaction. This result is consistent with other DEM results presented in the literature [40]. Atε = 0.08 the structure starts to fail and develops into a dilated particle structure (corresponding to a positive value ofJ rel ). At this stage, the order of the deformation responses in Fig. 8b changes, and is expected to eventually show a trend whereby the dilation will be larger for a higher particle friction, see also [6,40]. However, this will happen at deformations falling outside the range considered here, whereby shear failure has substantially developed.
(a) (b) Fig. 8 a Average macroscopic stress ratioσ and b volumetric deformationJ rel as a function of the applied vertical strainε for different particle contact friction coefficients μ The evolutions of the mean coordination numbern and the mean anisotropy parameter A for simulations performed with different particle friction coefficients are illustrated in Fig. 9a, b, respectively. For all cases the coordination number decreases with increasing deformation; the loss of contacts can be ascribed to the horizontal expansion of the specimen, which occurs despite that the overall material structure compacts, see Fig. 8b. A similar observation was made in the DEM study presented in [41]. Moreover, the coordination number in general is lower for a higher particle contact friction, which is in correspondence with the DEM study presented in [40]. Since the coordination number is lower and the overall strength is larger at higher particle contact friction, the average particle contact force becomes larger when the particle contact friction increases. Observe further that the effect of the contact friction on the induced fabric anisotropy is minimal during the initial stage of the response. Close to progressive shear failure some differences emerge, whereby the induced anisotropy tends to grow when the particle contact friction, and thus the overall failure strength, becomes larger.

Particle rotation
The effect of particle rotation upon the macroscopic response of the granular structure is examined by comparing the responses of the reference particle structure defined in Table 3 for the cases with and without particle rotation. The latter case is obtained by prescribing the rotation of the inner particles P p of the granular micro-structure to be zero, θ p = 0, throughout the entire loading process. For the boundary particles P q the rotation follows from the periodicity requirement (12). Figure 10a shows the evolution of the stress ratioσ and Fig. 10b illustrates the development of the volumetric defor-mationJ rel for the two cases. The ultimate failure strength of the system with constrained particle rotation appears to be more than 30% higher than the ultimate failure strength of the system with unconstrained particle rotation. Essentially, limiting the particle rotation may be interpreted as a kinematic constraint that increases the shear strength [6,42]. Hence, granular materials composed of angular shaped particles, which are susceptible to interlocking, show limited particle rotation, and thus typically have a higher effective strength than granular materials composed of smooth, round particles [43]. The deformation behaviour plotted in Fig. 10b illustrates that the particle system with constrained rotation experiences more compaction than the particle system with unconstrained rotation. This result is consistent with the DEM study on relatively soft particle systems reported in [42], and is due to larger particle overlap from the higher normal contact forces generated under constrained particle rotation. Note that the above trends are analogue to what has been observed in Fig. 8 from constraining particle sliding by increasing the particle contact friction. Figure 11a, b show the evolution of the coordination numbern and the induced fabric anisotropy A, respectively, for the simulations with constrained and unconstrained particle rotations. The coordination number decreases with deformation, and is smaller for the case of constrained particle rotation. Furthermore, from Figs. 10a and 11b it is concluded that the shear strength is higher when the induced anisotropy is larger. More specifically, at an axial strain ε = 0.073 the granular assembly with unconstrained particle rotation has reached its maximal shear strength, which is 0.39/0.48 = 0.81 times lower than the corresponding strength of the granular assembly with constrained particle rotation. The anisotropy ratio for the two granular assemblies turns out to be similar, and atε = 0.073 equals 0.27/0.33 = 0.82, thus indicating here an almost linear relation between shear strength and induced anisotropy.

Initial anisotropy
The influence of the initial anisotropy of the granular micro-structure on the macro-scale response is assessed by considering, together with the reference particle packing defined in Table 3, two additional particle packings characterized by higher anisotropy values. Accordingly, the set of initial anisotropy parameters is A 0 = [0.02, 0.05, 0.08], and the corresponding microstructures and rose diagrams are sketched in Fig. 12.
The evolutions of the effective stress and the volumetric deformation for the different initial fabric anisotropies are shown in Fig. 13a, b, respectively. The value of the peak strength appears to be lower for a higher initial fabric anisotropy, and is reached at a smaller axial strainε. Further, in the pre-peak regime a higher initial fabric anisotropy creates less compaction of the granular packing, which agrees with the DEM study performed in [41]. Combining the results from both figures leads to the conclusion that a higher compaction level obtained under a lower initial anisotropy provides a higher effective shear strength of the sample.
The average coordination numbern and the induced anisotropy A are illustrated in Fig. 14a, b for the three different initial anisotropies. The structure with the lowest initial anisotropy A = 0.02 clearly has the highest coordination number close to the onset of localized failure at ε = 0.07, which corresponds to the highest compacting level and the largest overall strength, see Fig. 13. For the other two anisotropies considered this relation is less clear, but may become more apparent when the number of particles in the micro-scale domain is enlarged. Further, the link between the initial anisotropy and the induced anisotropy depicted in Fig. 14b can not be clearly established, although it can be observed that during progressive shear failure (i.e., in the deformation range 0.06 ≤ε ≤ 0.08) the order of the three curves is similar as for the shear strength depicted in Fig. 13a.

Conclusions
In the present contribution a multi-scale model for the analysis of granular systems has been proposed, which combines the principles of a coupled FEM-DEM approach with a novel servo-control methodology for the implementation of appropriate micro-scale boundary conditions. A mesh convergence study has been performed, whereby the results of a quasi-static biaxial compression test were compared to those obtained by direct numerical simulations. The comparison demonstrated the capability of the multi-scale method to realistically capture the macro-scale response, even for macroscopic domains characterized by a relatively coarse mesh; this makes it possible to accurately analyse largescale granular systems in a computationally efficient manner. The multi-scale framework has been applied to study in a systematic manner the role of individual micro-structural characteristics on the effective macro-scale response. The effect of particle contact friction, particle rotation, and initial fabric anisotropy on the overall response has been considered, as measured in terms of the evolution of the effective stress, the volumetric deformation, the average coordination number and the induced anisotropy. The trends observed are in accordance with notions from physics, and observations Since the proposed multi-scale framework is based on first-order homogenization principles, it can only be adequately applied for problems whereby microscopic length scale effects do not influence the macroscopic response. Examples whereby this separation of length scales holds are static (and dynamic) problems in which significant strain localization remains absent, and dynamic problems in which the time-dependent response is composed of non-dispersive, slowly varying low-frequency components. The extension of the proposed multi-scale FEM-DEM scheme for applications related to strain localization and high-frequency wave propagation is a topic for future research.