Formulation and numerical implementation of micro-scale boundary conditions for particle aggregates

The aim of this paper is to propose a novel methodology to deal with micro-structural boundary conditions for the analysis of granular materials. The response of the granular assembly is modelled through the discrete element method (DEM), consistent with a micro-macro homogenization scheme formulated in the transition from a continuous to a discrete setting. The three classical types of boundary conditions (displacement, periodic and displacement) are implemented with a servo-control method. Two types of algorithms are developed, which consider or not an initial prediction for the displacement field of the particles. Additionally, to simulate in a more realistic manner the effect of macroscopic deformation on the micro-structure, a mixed type of boundary conditions is introduced, i.e. a combination of classical boundary conditions. The effectiveness of the algorithms is tested on a series of DEM simulations. The influence of the micro-structural morphology and of the size of the micro-structure on the overall response is finally investigated.


Introduction
The accurate computation of the non-linear failure and deformation behavior of heterogeneous granular systems commonly requires a resolution of the complex mechanical interactions and deformation mechanisms at the particle scale, which can be adequately accounted for by using the discrete element method (DEM), see e.g., [1,2,3,4,5,6,7,8,9,10] and references therein. For practical granular systems composed of a vast number of particles, however, it is infeasible to simulate each particle as an individual discrete object, since this leads to DEM models with an enormously large number of degrees of freedom, and consequently, to impractical computation times. To circumvent this problem, advanced multi-scale frameworks have been developed, where the mechanical responses at the particle micro-scale and the structural macro-scale are hierarchically coupled in an computationally economical fashion. This is accomplished by simulating the macro-scale problem under consideration with the finite element method (FEM), whereby in every integration point the response to the corresponding deformation is calculated by means of a DEM model that accurately and efficiently represents the complex particle behavior at the micro-scale. Examples of coupled FEM-DEM approaches for granular materials can be found in [11,12,13,14,15], illustrating the use of various averaging theorems for relating force and displacement measures at the particle micro-scale to stress and strain measures at the structural macroscale. Specific aspects that should deserve more attention in FEM-DEM homogenization methods, but often are neglected for reasons of simplicity, refer to i) the Hill-Mandel micro-heterogeneity condition, which enforces consistency of energy at the micro-and macro-scales, ii) the effect of particle rotations in the formulation of micro-to-macro scale-transitions, and iii) a rigorous generalization of the multi-scale approach within the theory of large deformations.
The computational homogenization framework presented in [16] does include the three aspects mentioned above, and calculates the micro-scale response of a granular packing with a DEM model equipped with a frame of boundary particles at which the finite deformation following from the macro-scale is imposed. The formulation considers three types of micro-scale boundary conditions for the boundary particles, namely i) homogeneous deformation and zero particle rotation (D), ii) periodic particle displacements and rotations (P), and iii) uniform particle force and free particle rotation (T), where the abbreviations (D), (P) and (T) are adopted from analogous, classical boundary conditions used in continuum homogenization theories, referring to the displacement, periodic and traction boundary conditions, respectively. The numerical implementation of the boundary conditions in [16] is done via a penalty method, where the violation of the boundary conditions is punished by increasing the total virtual work of the particle packing, through the introduction of additional forces and moments on the frame of boundary particles. Although the algorithm presented in [16] has been nicely generalized for the three types of boundary conditions in a mathematically elegant and transparent fashion, due to the nature of the penalty method the expression for the homogenized stress of the particle packing becomes explicitly dependent on the value of the penalty parameter, and thereby looses its physical interpretation. In addition, in DEM models the penalty parameter may be difficult to control and must be chosen sufficiently large in order for the penalty function to be effective, which may induce numerical instabilities [17,18]. Another characteristic of the penalty method is that it requires the constraint equations to be satisfied "approximately" instead of "exactly", whereby the accuracy of the approximation is determined by the magnitude of the penalty parameter. As a consequence, the limit case at which the boundary conditions are met exactly is not rigorously retrieved from the formulation, since the homogenized stress of the particle packing then vanishes, see expression (44) in [16].
In order to improve on the algorithmic drawbacks mentioned above, in the present communication an alternative numerical algorithm is proposed for the implementation of the homogenization framework presented in [16]. This algorithm is based on a servo-control methodology, using a feedback principle comparable to that of algorithms commonly applied within control theory of dynamic systems [19]. Accordingly, the displacements and rotations of the particles of the boundary frame are iteratively adapted from a gradually diminishing discrepancy between the measured and desired values of the micro-scale boundary condition. A strong point of this approach is that it is relatively simple to implement, and only affects the interface communicating information between the macro-scale FEM and micro-scale DEM models. In other words, it does not require internal modifications of the FEM and DEM source codes, so that the approach can also be combined with commercially available software for which the user typically has no access to the source code. In addition to its simplicity, the servo-control algorithm preserves the physical interpretation of the homogenized stress measure derived for the particle packing, and furnishes a realistic value for the stress in the limit case at which the micro-scale boundary condition is met exactly. It is noted that the algorithm only considers the (P) and (T) boundary conditions, since for a macro-scale problem discretized with a displacement-based FEM code, the (D) boundary condition can be implemented in a straightforward fashion, without the use of iterations.
Apart from providing servo-control algorithms for the individual (P) and (T) boundary conditions, a novel formulation for mixed (D)-(P)-(T) boundary conditions is derived, and subsequently cast into a numerical formalism. The formulation is proven to satisfy the Hill-Mandel micro-heterogeneity condition, and therefore is very useful for i) a consistent derivation of macro-scale constitutive relations from standard material tests on particle aggregates subjected to any combination of (D)-, (P)-and/or (T)-type boundary conditions, and ii) the efficient computation of the homogenized response of large-scale particle aggregates characterized by a spatial periodicity in one or two directions, i.e., granular layers exposed to uniform (D) and/or (T) boundary conditions at their top and bottom surfaces. It will be demonstrated that the formulation allows to impose the (D) and (T) boundary conditions both at separate and identical parts of the layer boundary, where in the latter case the (D) and (T) contributions obviously need to be applied along different orthonormal directions.
The performance of the servo-control algorithms developed for the various micro-scale boundary conditions is tested by using monodisperse and polydisperse frictional and cohesive packings composed of two-dimensional, circular particles and subjected to various loading paths. These examples illustrate the basic features of each of the boundary conditions in full detail. Despite the focus on two-dimensional particle systems, it should be mentioned that the extension of the present framework towards three-dimensional granular systems is trivial, and can be made without the introduction of additional prerequisites.
The paper is organized as follows. Section 2 presents a review of the numerical homogenization framework for particle aggregates, and outlines the formulations of the micro-scale (D), (P) and (T) boundary conditions proposed in [16]. Section 3 discusses the numerical implementation of the micro-scale boundary conditions, where for the (P) and (T) boundary conditions two different servo-control algorithms are presented, which include or not an initial prediction of the displacements of the boundary particles based on their positions calculated at the previous loading step. In Section 4 the performance of the numerical algorithms is tested on monodisperse and polydisperse frictional packings subjected to various loading paths. The numerical results clearly illustrate the characteristic differences in response for the three types of boundary conditions, and show their response convergence behavior under increasing sample size. Section 5 presents the formulation for the mixed boundary conditions, and provides the details of the servo-control algorithm and its numerical performance for the cases of infinite frictional and cohesive granular layers loaded by a vertical compressive stress, and subsequently subjected to a horizontal shear deformation. Some concluding remarks are provided in Section 6.
In terms of notation, the cross product and dyadic product of two vectors are, respectively, designated as a × b = e ijk a i b j e k and a ⊗ b = a i b j e i ⊗ e j , where e ijk 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 product between two vectors is given by a · b = a i b i , and between two second-order tensors by A : The action of a second-order tensor on a vector is indicated as A · b = A ij b j e i . The superscript T is used to indicate the transpose of a vector or a tensor. Further, I = δ ij e i ⊗ e j denotes the second-order identity tensor, with δ ij the Kronecker delta symbol.
Since the present study focuses on two-dimensional particle aggregates, throughout the paper 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.
2 Micro-macro transitions for particle aggregates

Micro-scale geometry
The initial micro-scale granular system is characterized by a two-dimensional square domain of P + Q rigid particles, which are partitioned into P inner particles P p , with p = 1, .., P , and Q boundary particles P q , with q = 1, .., Q, colored in yellow and red in Figure 1(a), respectively. 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 initial interior domain V comprises the inner particles with their center points as X p ∈ P p with p = 1, .., P . The boundary ∂V is defined by the boundary particles, whose center points in the initial configuration are X q ∈ P q with q = 1, .., Q. The macroscopic deformation of the granular microstructure is imposed via the frame of boundary particles P q , as a result of which the center points of the inner and boundary particles become located at x p and x q , respectively, see Figure 1(b). In the current configuration, the boundary particles P q are subjected to boundary forces a q , boundary moments m q , and particle contact forces f c q , see Figure 1(c), while the inner particles P p are subjected to particle contact forces f c p , see Figure 1(d), with the superscript c denoting a contact with a neighbour particle.
The macroscopic response of a granular assembly is derived by transforming relevant principles used in firstorder homogenization theories [20,21,22,23] from a continuous setting to a discrete setting. Accordingly, at the centroids of the boundary particles P q the finite area vectors A q and forces a q are derived from infinitesimal area vectors and forces, respectively, with N the vector pointing in the outward normal direction of the boundary ∂V of the initial particle volume V , t being the boundary traction, and ds indicating an infinitesimal part of the boundary surface. Various expressions for A q have been presented in the literature, see e.g., [24,25]. In the present communication the initial area vector is computed by accounting for the different radii of the boundary particles: where R q+1 , R q and R q−1 are the radii of adjacent boundary particles q + 1, q and q − 1, respectively. Further, e 3 represents the unit vector in the out-of-plane direction of the two-dimensional particle structure, see also Figure 1(a). Note that in (2) the boundary particles must be numbered in the anti-clockwise direction in order to obtain an area vector pointing in the outward normal direction of the boundary. Figure 1: (a) Two-dimensional particle aggregate of initial volume V and boundary ∂V . Yellow and red colors refer to inner P p and boundary P q particles, respectively; (b) Particle aggregate in the current configuration; (c) Boundary forces a q , boundary moments m q , and particle contact forces f c q acting on the boundary particles P q ; (d) Particle contact forces f c p acting on the inner particles P p .

Micro-scale governing equations 2.2.1 Equilibrium conditions
In the absence of body forces, the mechanical equilibrium of a granular microstructure can be formulated in terms of the boundary forces a q and moments m q acting on the boundary particles P q : with x q the current position vector of the boundary particles. The boundary forces a q and moments m q thus drive the overall, macroscopic deformation of the granular system via the frame of boundary particles. Note that, besides global equilibrium (3), local equilibrium conditions may be formulated for each of the inner particles P p , which interact through contact forces f c p at discrete contact points x c p on the particle surfaces: with the superscript c referring to a particle contact, N c p being the number of contact forces related to particle p and x p is the current position vector of the inner particle. Analogous conditions may be written for the boundary particles P q , for which discrete contact forces f c q act at contact points x c q on the particle surfaces. The frame of boundary particles is driven by boundary forces a q and boundary moments m q , where N c q is the number of contact forces for particle q. Note that the combination of expressions (4) and (5) is in correspondence with relation (3).

Particle contact laws
In order to solve the micro-scale problem, the constitutive response of the granular assembly needs to be defined through a relation between the contact forces f c i (or contact moments m c i ), with i = 1, ..., P + Q, and the corresponding contact displacements ∆u c i (or contact rotations ∆θ c i ). For the sake of clarity, in the following the superscript c and subscript i will be dropped. Two types of particle contact interactions will be considered, which are referred to as frictional contact and cohesive contact.
in accordance with the following definitions where ∆u c n and ∆u c s are the relative displacements in the normal and tangential direction of particle contact c and N c is the total number of particle contacts. Note that for the cohesive contact law given by equations (7) and (8) the potential energy in (12) needs to be extended with the rotational term k θ (∆θ c ) 2 /2.
Obviously, for deriving the solution of a boundary value problem, the equation of motion (9) and the constitutive response of the particles (6) and (7) should be complemented by the appropriate boundary conditions. As mentioned in the introduction, the numerical implementation of the micro-scale boundary conditions is based on the formulation presented in [16], and the main equations are summarized in Section 2.3 for the sake of clarity.

Micro-scale kinematics and boundary conditions
Consider a rigid particle i within a granular assembly. The current location x of an arbitrary material point, located within the initial particle volume at X, is defined through the non-linear deformation map where x i and X i are the current and original positions of the center of particle i, and Q i is the second-order particle transformation tensor. For plane problems defined with respect to the orthonormal tensor basis {e k ⊗ e l } 2 k,l=1,2 , the transformation tensor of particle i can be expressed as Q i = cos θ i e 1 ⊗e 1 −sin θ i e 1 ⊗e 2 +sin θ i e 2 ⊗ e 1 + cos θ i e 2 ⊗ e 2 , with θ i the magnitude of the particle center rotation θ i = θ i e 3 , where e 3 is the unit vector normal to the plane. In addition, the current position of the particle center x i can be expressed as the sum of a contribution affine to the macroscopic deformation gradientF and a local, micro-scale fluctuation w i : In homogenization schemes for continuous media, the macro-to-micro scale transition is enforced by requiring the macro-scale deformation gradient to be equal to the volume average of the micro-scale deformation gradient. In a discrete setting, this is equivalent to the condition Relation (15) can be derived by transforming the volume average of the macro-scale deformation into a surface integralF with N the vector normal to the outer boundary of the original particle volume, and subsequently performing the transition from a continuous to a discrete setting with the aid of (1) 1 . Equation (15) needs to be satisfied by applying specific boundary conditions to the boundary particles of the granular microstructure. For continuous media, this goal is typically accomplished by applying one of the three classical types of boundary conditions, namely i) a homogeneous deformation, also known as the displacement boundary condition and thus abbreviated as (D), ii) periodic displacements (P), and iii) a uniform traction (T), see, e.g., [21,22,23]. For discrete particle structures, however, additional conditions need to be imposed on the rotations or moments of the boundary particles. Correspondingly, along the lines of [16], the three boundary conditions mentioned above are extended as i) homogeneous deformation and zero rotation (D), ii) periodic displacement and periodic rotation (P), and iii) uniform force and free rotation (T), of which the formulations are presented below. The abbreviations (D), (P) and (T), although typically used in continuum homogenization theories, are maintained here for reasons of consistency. In addition to the three classical boundary conditions, a novel combination of these boundary conditions has been derived, which will be referred to as "mixed boundary conditions". The corresponding formulation is proven to satisfy the consistency of energy between the microscopic and macroscopic scales of observation, known as the Hill-Mandel micro-heterogeneity condition, and the details are provided in Section 5.

Homogeneous deformation and zero rotation (D)
In accordance with this boundary condition, all the boundary particles P q are prescribed to have zero microscale displacement fluctuations and zero rotations: where the first expression follows from (14) with the displacement fluctuations as w q = 0. Due to the second condition in (17) the boundary moments do not vanish, i.e., The homogeneous deformation and zero rotation boundary condition is expected to result in a relatively stiff macroscopic response of the particle aggregate.

Periodic displacement and periodic rotation (P)
For this boundary condition, both the displacements and rotations of the boundary particles P q are related by periodicity requirements: where the superscripts + and − refer to corresponding particles on opposite boundaries of the granular assembly. From the viewpoint of equilibrium, the forces and moments on opposite boundaries need to be anti-periodic, thus satisfying the relations

Uniform force and free rotation (T)
The boundary forces a q of the boundary particles P q are here determined from the product of the macroscopic first Piola-Kirchhoff stressP and the discrete area vectors A q introduced in equation (1) 1 : In addition, no constraint is applied to the boundary rotations, so that the boundary moments vanish: The uniform force and free rotation boundary condition is expected to provide a relatively soft macroscopic response of the particle aggregate.

Macro-scale stress and Hill-Mandel condition
The first Piola-Kirchhoff stress at the macro scale is defined in terms of the boundary forces a q acting on the particle aggregate:P The Hill-Mandel micro-heterogeneity condition expresses the equality between the volume average of the virtual work applied at the boundaries of the micro-structure and the virtual work of a macroscopic material point [31]. For a discrete particle system this condition specifies intō The macroscopic stressP given by (23) must satisfy the energy consistency between the two scales. Accordingly, considering definition (15), the following identity holds Alternatively, by making use of the definition of the macro-scale stress (23), the inner productP : δF can be expanded asP Subsequently, reformulating equation (24) as and substituting equations (25) and (26), leads to Invoking the micro-scale displacement fluctuations in accordance with relation (14) turns expression (28) finally into Note that the recast form (29) of the Hill-Mandel condition is satisfied for all three types of boundary conditions introduced above: For the (D) boundary condition, the combination of equations (14) and (17) results in δw q = 0. For the (P) boundary condition, the periodicity of the micro-fluctuations of the boundary displacements w + q = w − q and the anti-periodicity of the boundary forces a + q = −a − q , following from equations (19) and (20), respectively, make that their products in expression (29) vanish for opposite boundaries. For the (T) boundary condition, relation (21) leads to a q −P · A q = 0.
It should be mentioned that the Hill-Mandel condition elaborated above only accounts for the influence of contact forces acting on inner particles, and does not include the effect of contact moments. Although in the frictional contact law the contact moments are absent, see expression (6), for the cohesive contact law they contribute both to the elastic behavior and the strength criterion, see expressions (7) and (8), respectively. The extension of a contact law with a contact moment contribution formally introduces a couple stress in the macroscopic response of the particle aggregate, which is energetically conjugated to the gradient of the overall rotation, see e.g., [32,33,34,35,36,37]. These higher-order stress and deformation measures correspond to higher-order natural and essential boundary data [37], which are known to be difficult to measure in experiments, and commonly are (substantially) lower in magnitude than the classical boundary data. For these reasons, and from the fact that the cohesive contact law defined by equations (7) and (8) is used only for one example discussed at the end of this communication, see Section 5.3, an extension of the Hill-Mandel condition with the effect of a contact moments is omitted here, but may be considered as a topic for future research. Correspondingly, for the cohesive contact law a consistency in energy between the micro-and macro-scales of a particle aggregate can only be warranted in an approximate fashion.
In Sections 4 and 5, the results of the DEM analyses will be presented in terms of components of the macroscale Cauchy stress tensorσ. This stress measure can be derived from the first Piola-Kirchhoff stressP computed through (23) by using the common transformation rule:

Numerical implementation of micro-scale boundary conditions
The micro-scale boundary conditions outlined above were implemented by using the open-source discrete element code ESyS-Particle [38,39]. The numerical algorithms developed for this purpose are described below.

Homogeneous deformation and zero rotation (D)
The homogeneous deformation and zero rotation boundary condition (D) given by equation (17) can be implemented straightforwardly by imposing this condition in an incremental fashion on the boundary particles P q . After moving the boundary particles in accordance with the incremental update of the deformationF , dynamic relaxation is applied to reach the equilibrium state of the particle aggregate, during which the displacements imposed on the boundary particles remain fixed. The particle configuration corresponding to the equilibrium state is stored, and the next deformation increment is applied. This process is repeated until the total number of deformation increments i tot is reached. The details of the algorithm are summarized in Table 1.
1. DEM simulation. Increments 0 ≤ i ≤ i tot 1.1 Apply updated boundary conditions x q =F X q and Q q = I for q = 1..., Q 1.2 Dynamic relaxation 1.3 Save current configuration and go to 1 (next increment i + 1) Table 1: Algorithm for the (D) boundary condition.

Periodic displacement/periodic rotation (P) and uniform force/free rotation (T)
The periodic displacement and periodic rotation boundary condition (P) and the uniform force and free rotation boundary condition (T) were numerically implemented by means of a servo-control algorithm, which uses a feedback principle similar to that of algorithms commonly applied within control theory of dynamic systems [19]. More specifically, the algorithms iteratively correct the boundary particle displacements and rotations from a gradually diminishing discrepancy between the measured and the required values of the boundary condition.
For the periodic displacement and periodic rotation boundary condition (P), the boundary forces and boundary moments should satisfy the anti-periodicity conditions presented in equation (20). Accordingly, the corresponding residuals for the edge particles are Multiplying the residuals by corresponding gain parameters g p a and g p m results into the following displacement and rotation corrections for the edge particles: which are added to the particle locations and rotations from the previous iteration. Note that the four corner particles straightforwardly follow the macroscopic deformationF , by prescribing their displacements in accordance with equation (17). Hence, for these particles no displacement correction is needed. The rotations of the four corner particles will be updated similarly to (32), using the corrections For the uniform force and free rotation boundary condition (T), the boundary forces ensue from the applied macroscopic stress through relation (21), whereby the imposed macroscopic deformation is given by expression (15). The corresponding residuals can be formulated as The correction for the displacement of the boundary particles is derived by multiplying the force and deformation residuals in (34) by the gain parameters g t a and g t F , respectively, leading to When performing numerical simulations, the specific values of the gain parameters g p a , g p m , g t a , g t F need to be fine-tuned from accuracy and stability considerations of preliminary numerical benchmark tests.
The corrections for the displacement and rotation of the boundary particles were implemented by means of two different algorithms, which consider or not an initial prediction of the position of the boundary particles based on their positions calculated at the previous loading step. These algorithms are therefore given the labels "with initial displacement prediction" and "without initial displacement prediction". The algorithms are discussed below, and their effect on the computational results will be investigated in Section 4. The specific parts of the algorithms that refer to the periodic displacement and periodic rotation boundary condition will be denoted by the symbol (P), while the symbol (T) indicates the uniform force and free rotation boundary condition. Finally, the residuals defined in expressions (31) and (34), which relate to the particle force, particle moment and macroscopic deformation gradient, are evaluated at each iteration by subjecting their dimensionless form to a convergence check. The dimensionless forms are obtained through, respectively, a normalization by the following parameters:ã with k = c, e, q referring to corner, edge, and boundary particles, respectively. In (36), M k is the mass of particle k, R k is its radius and ∆t is the time increment used in the dynamic relaxation procedure.

Algorithm with initial displacement prediction
The macroscopic deformation is imposed in i tot steps on the boundary particles P q via the incrementally updated deformation gradientF . In correspondence with the algorithm presented in Table 2, in the initialization step, i = 0, the boundary particles are moved in accordance with a homogeneous deformation, and for the periodic boundary also a zero rotation, similar to equation (17). Subsequently, the granular assembly is dynamically relaxed to the equilibrium state, keeping the translational and, for the periodic boundary, rotational degrees of freedom of the boundary particles fixed. The iterative loop is entered, and the actual values of the forces and moments of the boundary particles are recorded. For the (P) boundary condition, the corrections for obtaining periodic particle translations and rotations at the boundary are calculated for the corner and edge particles separately, in accordance with relations (31)- (33). For the (T) boundary condition, the boundary moments vanish and the displacement corrections are computed via (34)- (35). The residuals are computed and compared with prescribed tolerances. For the (P) boundary condition, the residual is based on boundary forces and moments. For the (T) boundary condition, two residuals are calculated, which are based on the boundary forces and on the imposed macroscopic deformation. If the norm of the residual(s) is(are) smaller than the tolerance(s) (referred to as a for the force criterion and F for the deformation criterion), the iterative loop is terminated and the next loading step is applied. If the convergence criterion is not satisfied, the corrections are computed again and the residual is iteratively re-examined, until convergence is reached.
After the initialization step is concluded, the responses for subsequent increments, 1 ≤ i ≤ i tot , are calculated, see Table 2. For the (P) boundary condition, the corner nodes are moved by straightforwardly imposing the updated macro-scale deformation in accordance with relation (17). For the edge particles, their current position is determined from a prediction based on the particle position in the previous loading step i − 1. More specifically, this prediction is a function of the position a particle would have in case of a homogeneous deformation (using the displacement boundary condition (17)), plus the difference, multiplied by an inheritance factor n f , between the final particle position at the previous increment and the position the particle would have at the previous increment under a homogeneous deformation. The inheritance factor lies between 0 and 1, and its optimal value depends on the loading conditions applied and the characteristics of the particle assembly. For the (T) boundary condition, the prediction occurs in an analogous fashion and is applied to all the boundary particles. After the boundary particles are translated in accordance with the predicted values of their positions, the granular assembly is dynamically relaxed to its equilibrium state. Subsequently, the iterative loop is entered, which invokes the previously described correction procedure of the displacements and rotations, in correspondence with the servo-control methodology.

Algorithm without initial displacement prediction
Similar to the algorithm with initial displacement prediction, for the algorithm without initial displacement prediction the macroscopic deformation is imposed in i tot steps to the boundary particles P q . However, as pointed out in Table 3, all increments are now treated in the same fashion. The boundary particles are initially moved in accordance with the updated homogeneous macroscopic deformationF , similar to expression (17), after which the particle assembly is dynamically relaxed to its equilibrium state. The iterative loop is started, in which the corrections for the displacement and rotation of the boundary particles are calculated based on the servo-control methodology. For the (P) boundary condition, the boundary is partitioned into corner and edge particles, whereby relations (31)-(33) are applied. For the (T) boundary condition, equations (34)-(35) are employed. Subsequently, the particle system is relaxed to the equilibrium state, and the current values of the boundary forces and moments are recorded and used to compute the residuals. If the norms of the residuals are smaller than the corresponding tolerances adopted, the iterative loop is terminated and the next loading step is applied. It can be confirmed that the algorithm without displacement prediction can be obtained as a limit case of the algorithm with initial displacement prediction by setting the inheritance factor equal to zero, n f = 0, whereby the algorithmic structure provided in Table 2 reduces to the more compact and simpler algorithmic structure presented in Table 3.

Computational results for regular and irregular packings
The algorithms proposed above for the implementation of the micro-scale boundary conditions are tested on a series of DEM simulations on regular, monodisperse and irregular, polydisperse particle packings. =F X e and the inheritance factor 0 < n f ≤ 1

2.1.B if (T) =⇒
Prediction of the positions of boundary particles: =F X q and the inheritance factor 0 < n f ≤ 1 2.  Table 2: Algorithm for the (P) and (T) boundary conditions with initial displacement prediction.

Regular monodisperse packing
In this section the responses of three different regular, monodisperse particle packings are considered, which consist of circular particles of radius R = 1. ∆a q · ∆a q /ã 2 q and r F = Q q=1 ∆F q · ∆F q /F 2 q 1.6 Check for convergence: r a ≤ a for (P); r a ≤ a and r F ≤ F for (T) 1.6.A if converged =⇒ Save current configuration and go to 1 (next increment i + 1) 1.6.B if not converged =⇒ Return to 1.3 Table 3: Algorithm for the (P) and (T) boundary conditions without initial displacement prediction.
The particles obey a frictional contact law, in correspondence with relation (6). Assuming relatively soft particles, the normal and tangential stiffnesses are chosen as k n = 10 4 N/m and k s = 2 · 10 3 N/m, 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. The packings are subjected to a combined biaxial compression-true shear deformation F = I +F 11 e 1 ⊗ e 1 +F 12 e 1 ⊗ e 2 +F 21 e 2 ⊗ e 1 +F 22 e 2 ⊗ e 2 , withF 11 =F 22 = −0.03 andF 12 =F 21 = −0.3, which is applied in i tot = 300 loading steps. For reaching the equilibrium state at each loading step, the particle system is subjected to dynamic relaxation steps of constant time increments ∆t = 10 −6 s. The gain parameters (in dimensionless form) used for the correct application of the boundary conditions are: for (P) g p a M/∆t 2 = 1 · 10 2 and g p m M R 2 /∆t 2 = 2 · 10 2 ; for (T) g t a M/∆t 2 = 1 · 10 2 and g t F R 2 = 2 · 10 −5 , with M = ρπR 2 representing the mass of the particles. The force and deformation tolerances are taken as a = 10 −6 and F = 10 −2 , respectively. For the dynamic relaxation process, a value of 10 −3 is adopted for tol E , whereby equation (11) must be minimally satisfied for a pre-defined, continuous period of 20∆t, in order to ensure a rigorous dynamic relaxation to the equilibrium state. An overview of the model parameters is given in Table 4.

Responses for algorithms with and without initial displacement prediction
In order to investigate the performance of the two algorithms presented in Tables 2 and 3, the packing of 25 particles is considered first. The stress responses under the combined biaxial compression-true shear loading were computed with Eq. (30), and plotted as a function of the applied macroscopic shear deformationF 12 . Figures 2(a) and (b) show the results for the algorithms with (solid line) and without (dot-dashed line) an initial displacement prediction for the (P) and (T) boundary conditions, respectively. The normal and shear components of the Cauchy stress are normalized asσ 11 =σ 11 R/k n andσ 12 =σ 12 R/k n , respectively, wherē σ 11 andσ 12 are the macroscopic normal and shear Cauchy stresses of the particle aggregate. For the periodic boundary conditions (P), Figure 3 illustrates the packing structures at specific macroscopic shear deformations F 12 = −0.01 (a),F 12 = −0.113 (b), andF 12 = −0.28 (c). The red lines plotted in the deformed particle aggregates indicate the network of normal contact forces between the particles.
For the algorithm with initial displacement prediction, the local minimum of the normal stressσ 11 near F 12 = −0.113, as shown in Figure 2(a) for the periodic boundary conditions (P), can be ascribed to a joint localized sliding of all the boundary particles, see Figure 3(b), top. This localization mechanism does not arise for the algorithm without initial displacement prediction, which furnishes a shear response that is much more homogeneous, see Figure 3(b), bottom. It may be therefore concluded that the response of the packings is
rather sensitive to bifurcations in the equilibrium path followed, which here become evident due to the relatively low number of particles present in the packing. Under continuing deformation towardsF 12 = −0.28, the inner particles of the aggregate also develop substantial sliding, such that for both algorithms the particle structure gradually reaches its densest packing structure, at which the deformation as well as the normal contact force network become strongly homogeneous. Both for the normal and shear stress components the responses computed by the two algorithms at this stage have coalesced, and steadily grow under further increasing deformation.
A similar trend can be observed for the normal and shear stress responses of the particle aggregates with the (T) boundary condition, see Figure 2(b). The discrepancies in the responses computed by the two algorithms appears to be less than for the (P) boundary condition. Since the computational robustness and efficiency of the two algorithms is comparable in the present simulations, no clear distinction can be made in terms of their overall performance. Hence, as an arbitrary choice, the forthcoming DEM results are computed with the algorithm without initial displacement prediction.

Responses for the (D), (P) and (T) boundary conditions
The influence of the choice of the boundary condition on the overall packing response is illustrated in Figure  4 rotation boundary condition (P) with a dot-dashed line, and for the uniform force/free rotation boundary condition (T) with a dashed line. The stress response computed for the (P) boundary condition is bounded by the stiffer and softer responses measured for the (D) and (T) boundary conditions, respectively, a result that is in agreement with the numerical studies performed in [16]. It may be observed that the initial stress value corresponding to the (T) boundary condition is somewhat smaller than the value computed for the other two boundary conditions. This is, because after the sample preparation procedure was finished, the boundary forces generated by a 4% particle overlap do not exactly satisfy equation (21), and therefore in the first loading increment are slightly relaxed by the algorithm in order to meet this condition. Consider the average normalized particle overlap ∆ū n , defined as where N c is total number of particle contacts,R c is the average radius at contact c and ∆u c n is the particle overlap at contact c. Figure 4(b) depicts ∆ū n as a function of the applied shear deformation −F 12 . Observe that the trend for the average particle overlap is similar to that for the macroscopic normal stress in Figure 4(a). This can be explained as follows: The macroscopic stress is represented by the volume average of all contact forces generated within the granular microstructure. Since the assumed normal contact stiffness is considerably larger than the shear contact stiffness, k n >> k s , see Table 4, the normal contact forces f c n , which are proportional to contact overlaps ∆u c n , see expression (6), provide the main contribution to the macroscopic stress response.

Responses for different sample sizes
The effect of the sample size on the macroscopic stress response is considered by plotting the computational results for packings composed of 25, 100 and 225 particles for the (D), (P) and (T) boundary conditions in Figures 5(a), (b) and (c), respectively. Generally, for a larger sample the normal stressσ 11 decreases. The (D) and (P) boundary conditions show a close resemblance in the responses for 100 and 225 particles, from which it may be concluded that for a sample of about 225 particles the stress response has more or less converged. Conversely, for the (T) boundary condition the stress for a sample of 225 particles shows a substantial relative drop in value up to a deformation ofF 12 ≈ 0.20. This softening behavior appears to be governed by strongly localized deformations emerging at the boundaries of the particle system, a phenomenon that also has been reported for continuum homogenization methods equipped with this relatively soft boundary condition, see [40].

Irregular polydisperse packing
The irregular polydisperse packings analyzed in this section are composed of circular particles, with the particle radii arbitrarily taken from a uniform size distribution with polydispersity R max /R min = 2, where R min = 0.67 mm. A collision-driven molecular dynamics code described in [41] is used to randomly generate irregular packings with the initial number of particles equal to n 0 p = [25,100,200,400,600], as shown in Figure 6(a). Subsequently, these packings are reconstructed into geometrically periodic packings by copying each of the boundary particles intersecting with the edges of the square particle volume to corresponding positions at the opposite boundaries. This results into the packing geometries shown in Figure 6(b), with the final particle numbers as n p = [37,120,228,444,650]. The initial volumes of the particle aggregates are equal to V = [100, 400, 818, 1462, 2156] mm 2 , respectively. The rose diagrams of the particle assemblies are sketched in Figure 6(c), clearly indicating that the packings become more isotropic when the particle number increases. The particle volume fraction of the packings varies in the range v ∈ [0.833, 0.850], where the smallest and highest values correspond to the packings with the smallest and highest number of particles, respectively. The corresponding coordination numbers lie in between 2.97 and 3.47.
(a) (b) (c) Figure 6: Characteristics of the five different irregular polydisperse packings studied: (a) Initial packings generated by a collision-driven molecular dynamics code [41], (b) geometrically periodic packings with the number of particles equal to n p = [37,120,228,444,650], and (c) the rose diagrams.
The particle packings are subjected to a simple shear macroscopic deformation withF 12 = 0.5, which is applied in i tot = 100 loading steps. For the irregular packings the same physical and algorithmic parameters are used as for the regular packings, see Table 4, except for the tolerance F = 10 −1 and the two gain values for the (T) boundary condition, which here relate to g t a M i /∆t 2 = 5 and g t F R 2 i = 2 · 10 −6 , with M i = ρπR 2 i being the mass of particle i and R i its radius. Note that for an irregular polydisperse packing the specific gain values depend on the characteristics of the actual particle i.

Responses for the (D), (P) and (T) boundary conditions
The response of a packing with 228 particles is considered first. The normalized macroscopic stressesσ 11 = σ 11R /k n andσ 22 =σ 22R /k n , with the average radius asR = Figure 7 as a function of the applied macroscopic shear deformationF 12 . The solid, dot-dashed and dashed lines refer to the (D), (P) and (T) boundary conditions, respectively. Since for the packing of 228 particles the particle structure is rather isotropic, see Figure 6(c), it can be confirmed that the responses for the two normal stresses σ 11 andσ 22 indeed are similar. As for the regular monodisperse packing, the (D) and (T) boundary conditions provide the upper (stiffest) and lower (softest) bounds for the particle system response, and thereby encapsulate the response calculated for the (P) boundary condition. Although not depicted here, the responses for the normalized shear stressσ 12 =σ 12R /k n under the (D), (P) and (T) boundary conditions follow similar trends as observed for the normal stressesσ 11 andσ 22 , with the magnitude of the shear stress being about one third of that of the normal stresses. red lines. Obviously, this is in correspondence with the largest and smallest stress levels for the (D) and (T) boundary conditions, as depicted in Figure 7. Figure 9 illustrates the average normalized particle overlap ∆ū n , defined by relation (38), and the average particle rotationθ both as a function of the applied macroscopic shear deformationF 12 . As for the regular monodisperse packings, the average normalized particle overlap is the largest for the (D) boundary condition and the smallest for the (T) boundary condition, and shows a similar evolution as observed for the normal stresses, see Figure 7. As expected, the average rotation shows the opposite trend, being the largest for the soft (T) boundary condition and the smallest for the stiff (D) boundary condition.

Convergence behavior of macroscopic response under increasing sample size
The convergence behavior of the apparent macroscopic response of the particle aggregate towards its effective response under increasing sample size is studied by subjecting the five microstructures depicted in Figure 6(b) to a simple shear deformation given by (39). In convergence studies, this type of loading condition occasionally is characterized as "critical", because of a relatively slow convergence behavior towards a representative volume element (RVE). The convergence behavior is evaluated here by means of the L 2 -norm of the normalized, homogenized Cauchy stress tensorσ, integrated along the entire deformation path Figure 10 illustrates the stress norm σ L2 as a function of the sample size, expressed in terms of the number of particles. It can be observed that for the stiff (D) and soft (T) boundary conditions the stress norm, respectively, decreases and increases with increasing sample size, while for the periodic (P) boundary condition it remains approximately constant. These trends are typical for a change in apparent properties under increasing sample size, see e.g., [21]. However, for arriving at an RVE the curves for the (D) (P) and (T) boundary conditions must coincide [31], which indeed is not the case for the largest sample of 650 particles. As already indicated above, the minimum size of the RVE depends on the type of loading condition applied, which is known to be relatively large under a macroscopic shear deformation. From the approximately constant stress value observed in Figure 10 for the (P) boundary condition, it may be expected that the stress response of the RVE will be close to σ L2 ≈ 0.011. Hence, in multi-scale simulations on granular materials, the computational costs can be kept manageable by adopting the (P) boundary condition for a micro-structural sample of relatively small size, for which the homogenized response is similar to that of the minimal RVE with a much larger size.

Mixed boundary conditions
In this section the formulation and numerical implementation of mixed (D)-(P)-(T) boundary conditions is presented. The homogenization framework proposed satisfies the Hill-Mandel micro-heterogeneity condition, and thus can be used for i) a consistent derivation of macro-scale constitutive relations from standard material tests on particle aggregates subjected to any combination of (D)-, (P)-and/or (T)-type boundary conditions, and ii) the efficient computation of the homogenized response of large-scale particle aggregates characterized by a spatial periodicity in one or two directions, i.e., granular layers exposed to uniform (D) and/or (T) boundary conditions at their top and bottom surfaces. To the best of the authors' knowledge, the formulation presented is novel in the field of granular materials.

Formulation
For the formulation of the mixed boundary conditions, the basic particle configuration sketched in Figure 1 is considered, with the boundary being split up into the top part ∂V t , the bottom part ∂V b , the left part ∂V l  Figure 10: Stress norm σ L2 versus particle number for irregular polydisperse particle packings subjected to (D), (P) and (T) boundary conditions, in accordance with a macroscopic simple shear deformationF 12 . and the right part ∂V r . It is emphasized that the main concepts of the mixed formulation are general, and can be applied to arbitrary boundary value problems. The concepts are elaborated here for the specific case of an infinite horizontal layer of particles loaded by a constant vertical pressure,P 22 =P * 22 , and subsequently subjected to a shear deformationF 12 in horizontal direction. The reason for choosing this boundary value problem is that it includes all the three (D), (P) and (T) boundary conditions discussed previously, with their combinations entering the formulation both at separate and identical parts of the layer boundary. This allows for highlighting the characteristics of the mixed formulation in full detail. Accordingly, the macroscopic deformation of the particle aggregate is imposed via a combined (D)-(T) condition in which the macroscopic shear stressP 21 is measured from the response of the particle assembly. Note that (42) 1 implicitly accounts for the conditionF Furthermore, the first contribution in the right-hand side of (42) 2 typically is relatively small, since most particles q at the top boundary ∂V t are characterized by A q,1 << A q,2 , with A q,1 vanishing for the specific case of an ideally horizontal boundary composed of identical particles. Since the shear deformationF 12 is imposed after the application of the vertical stressP * 22 , in (42) 1 the reference positions X q of the boundary particles relate to the particle configuration obtained after the vertical stress has been applied. In summary, the boundary conditions for the particle aggregate are specified as follows: • For the particles that are part of the bottom boundary, q ∈ ∂V b , the homogeneous deformation and zero rotation boundary condition (D) is applied in accordance with expression (17). The vertical boundary displacements are constrained to construct a rigid support for the layer, and the horizontal boundary displacements follow the shear deformation given by expression (42) 1 .
• For the particles that are part of the left and right boundaries, q ∈ ∂V l ∪∂V r , the periodic displacement and periodic rotation boundary condition (P) is applied, as given by expression (19). This boundary condition reflects the horizontal confinement of the particles within the infinite layer.
• For the particles that are part of the top boundary, q ∈ ∂V t , free rotations are assumed, in correspondence with the (T) boundary condition 1 . For the description of the particle displacements, the boundary is split up along the two orthonormal directions e 1 and e 2 indicated in Figure 1. Along the e 1 -direction, the (D) boundary condition (42) 1 is applied to simulate the horizontal macroscopic shear deformation.
Along the e 2 -direction, a constant macroscopic pressureP * 22 is imposed via the (T) boundary condition (42) 2 , for which the corresponding components of the macroscopic deformation gradient,F 21 andF 22 , in accordance with the general form (15), turn intō Note that the two deformation components above should be considered as a computational result obtained by prescribing the stress componentP * 22 .
The macroscopic deformation gradientF , which is followed by the four corner nodes of the sample, now is fully specified through its "(D)-type components" provided by (42) 1 and (43), and its "(T)-type components" given by (44) 1,2 . The corresponding macroscopic Piola-Kirchhoff stress tensor is defined by equation (23). For the adopted mixed-boundary conditions it will now be demonstrated that this stress definition satisfies the recast Hill-Mandel condition given by expression (29), i.e., the energy consistency between the macro-and micro scales. Accordingly, relation (29) is first split up with respect to the different boundary parts considered above, i.e., For the bottom boundary ∂V b , boundary condition (D) holds, which, by comparing equations (17) and (14), lets the micro-fluctuations of the displacements vanish, w q = 0. Hence, the first term in (45) is equal to zero. At the left and right boundaries q ∈ ∂V l ∪ ∂V r , the (P) boundary condition is imposed, for which the microfluctuations of the displacements are periodic, w l q = w r q , see (14) and (19). Together with the anti-periodicity of the boundary forces a l q + a r q = 0, see (20), the second term in (45) vanishes. Finally, for the top boundary ∂V t , the last term in (45) may be further developed as q∈∂Vt a q,1 −P 1j A q,j δw q,1 + a q,2 −P 2j A q,j δw q,2 = 0 .
Along the e 1 -direction, the micro-fluctuations of the boundary particle displacements vanish, w q,1 = 0, in correspondence with equation (42) 1 , by which the first term in (46) becomes zero. Along the e 2 -direction, the boundary forces are uniform, a q,2 −P 2j A q,j = 0, see equation (42) 2 , so that the second term in (46) becomes zero. With this result, the Hill-Mandel condition (45) is proven to be satisfied for the mixed boundary conditions.

Numerical implementation
The numerical algorithm for the implementation of the mixed boundary conditions is outlined in Table 5, and is based on a combination of the algorithms presented in Section 3 for the (D), (P) and (T) boundary conditions, without an initial displacement prediction. During stage 1 of the loading process, the vertical compressive stress P 22 = P * 22 is applied to the particle aggregate in a stepwise fashion 2 , using a total of i vs loading increments, with the subscript vs designating "vertical stress". After initiating the displacement and rotation boundary conditions at the top ∂V t and bottom ∂V b boundaries, the vertical stress is incrementally updated and subsequently used to compute the displacement and rotation corrections at the left and right boundaries with expressions (31)- (32), and the displacement correction at the top boundary boundary with ∆u q,2 = g t a ∆a q,2 with ∆a q,2 = a q,2 −P 21 A q,1 −P * 22 A q,2 .
The expression above is derived from (42) 2 , whereby during the incremental application of the vertical stress P * 22 the value ofP 21 is prescribed as zero, in order to avoid the initial development of a shear stress. After the particle aggregate has reached its equilibrium state under dynamic relaxation, the boundary forces and moments of the particles at the top, left and right boundaries are recorded and employed to compute the corresponding residuals. When all residuals are lower than the prescribed values of the corresponding tolerances, the iterative loop is stopped and the next vertical stress increment is applied. Otherwise, the iterative loop is entered again, until a converged solution is found. After the application of i vs increments the vertical stress has reached the desired value, and stage 1 of the loading process has completed.
During stage 2 of the loading process, the horizontal shear deformationF 12 is imposed on the particles at the top ∂V t and bottom ∂V b boundaries of the granular assembly, by displacing these in a stepwise manner using i tot − i vs increments. The rotations of the particles at the top boundary are free, and the vertical displacement and rotation of the particles at the bottom boundary are fully constrained. In a similar way as explained above for stage 1, the boundary forces and moments in the relaxed equilibrium state are used to compute the displacement and rotation corrections at the periodic left and right boundaries ∂V l and ∂V r , and at the top boundary boundary ∂V t . However, the only difference is that in (47) the shear stressP 21 here is not prescribed as zero, but is calculated from the homogenized response of the particle assembly using equation (23). After the dynamic relaxation procedure has completed, the residuals are computed in the same way as during stage 1, and compared against the corresponding tolerances. The iterative process is terminated when the convergence criterion is satisfied, after which the shear deformation is incremented and the response to the next loading step is computed. This procedure is continued until all loading increments i tot are applied.  Table 5: Algorithm for the application of the mixed boundary conditions. The loading process consists of stage 1, during which the vertical stress is incrementally applied, and stage 2, during which the horizontal shear deformation is incrementally imposed.

Computational results
The performance of the algorithm used for the implementation of the mixed boundary conditions is demonstrated by means of two DEM simulations of an irregular polydisperse packing of 449 particles, with the particle radii taken randomly from a uniform size distribution with polydispersity R max /R min = 1.5, where the minimum radius equals R min = 0.8 mm. The initial particle volume is V = 1517 mm 2 , with the particle volume fraction of the packing being equal to v = 0.849, and the average particle coordination number as 3.55. The two simulations consider different particle contact laws, namely the frictional contact law and the cohesive contact law described in Section 2.2.2. Assuming relatively hard particles, the normal and tangential contact stiffnesses for the frictional contact law are set as k n = 10 5 N/mm and k s = 4 · 10 4 N/mm, respectively, and the friction coefficient equals µ = 0.6. The normal contact stiffness k b n and the tangential contact stiffness k b s for the cohesive contact interaction are assumed to be equal to those of the frictional contact law, and the bending contact stiffness is taken as k b θ = 2 · 10 4 Nmm. The normal, shear and bending strengths have the values f b,u n = 300 N, f b,u s = 60 N and m b,u θ = 200 Nmm, respectively. The density of the particles is ρ = 10 · 10 3 kg/m 2 . The macroscopic vertical (compressive) stress isP * 22 = −1.05 · 10 6 N/m, which is applied in i vs = 6 increments. The total macroscopic shear deformation equalsF 12 = 0.2, which is imposed on the particle aggregate in i tot − i vs = 100 increments. The translational and rotational damping factors used in the dynamic relaxation procedure are α = β = 0.7, and the time increment equals ∆t = 10 −5 s. The dimensionless values of the gain parameters are g t a M i /∆t 2 = g p a M i /∆t 2 = 3 · 10 4 and g p m M i R 2 i /∆t 2 = 6 · 10 4 , and the corresponding tolerances are equal to t a = p a = p m = 2 · 10 −10 . The above model parameters are summarized in Table 6. Nmm Density ρ 10 · 10 3 kg/m 2 Translational damping α 0.7 -Rotational damping β 0.7 -Time increment ∆t 10 −5 s Tolerance force (P) p a 2 · 10 −10 -Tolerance moment (P) p m 2 · 10 −10 -Tolerance force (T) t a 2 · 10 −10 -Gain force (P) g p a M i /∆t 2 3 · 10 4 -Gain moment (P) g p m M i R 2 i /∆t 2 6 · 10 4 -Gain force (T) g t a M i /∆t 2 3 · 10 4 -Tolerance dynamic relaxation tol E 10 −3 - Table 6: Physical and algorithmic model parameters for the simulations with mixed boundary conditions. Figure 11 shows the macroscopic response of the particle aggregates as a function of the applied shear de-formationF 12 , with the dot-dashed and solid lines referring to packings with cohesive and frictional particle contact interactions, respectively. In Figure 11(a) the stress ratioσ 12 /σ 22 is depicted, while Figure 11(b) illustrates the relative volumetric change det(F ) (using the packing obtained after the application of the vertical stress as the reference state), and Figure 11(c) sketches the average particle rotationθ, in accordance with expression (40). Furthermore, in Figure 12 the particle configurations of the cohesive and frictional packings are plotted at four different deformation levels, namely (a)F 12 = 0.002, (b)F 12 = 0.05, (c)F 12 = 0.1 and (d) F 12 = 0.015. Here, the red lines between the particles represent cohesive contact forces, while the blue lines indicate frictional contact forces. It can be observed from Figure 11(a) that up to a deformationF 12 = 0.02 the cohesive and frictional packings show a similar response, whereby the stress increases approximately proportionally with deformation. Upon continuing deformation, in the frictional packing a large number of contacting particles meets the failure criterion (6) and starts to slide, such that the stress ratioσ 12 /σ 22 reaches a maximum atF 12 ≈ 0.06. The maximal stress ratio is accompanied by a volumetric increase of the particle structure, commonly referred to as "dilation", see Figure 11(b). After passing the peak value, the stress ratio for the frictional packing slightly drops in magnitude, which is caused by a substantial rolling of particles. The effect of particle rolling can be clearly observed from Figure 11(c), where at the onset of shear deformation the increase in average particle rotation is only mild, but steadily grows towards a more or less constant value atF 12 = 0.07, both for the frictional and cohesive packings. Note from Figure 11(c) that the initial value of the average particle rotation is due to the application of the vertical stressP * 22 , and for the frictional packing appears to be somewhat larger than for the cohesive packing. For the cohesive packing the maximal value of the stress ratioσ 12 /σ 22 is about 1.5 times larger than for the frictional packing, and is reached atF 12 ≈ 0.07, see Figure 11(a). At this stage a significant number of particle bonds are broken, in correspondence with the failure criterion (8). With continuing deformation, the broken particle bonds of the cohesive packing become frictional, as indicated by the thick blue lines in Figure 12(c), whereby the stress ratioσ 12 /σ 22 of the packing drops to a similar level as for the frictional packing, see Figure 11(a). Observe from Figures 12(c) and (d) that for the cohesive packing the frictional contacts indicated by the blue lines are established only along local force chains in the particle structure, whereby the rest of the contacts remain cohesive, as indicated by the red lines. This implies that the overall deformation of the packing towards the end of the loading process becomes governed by a localized failure zone, which indeed is associated to a strong softening behavior in the stress response, see Figure 11(a). Close to a shear deformation ofF 12 = 0.20, both the stress ratioσ 12 /σ 22 and the relative volumetric change det(F ) of the frictional and cohesive packings become approximately constant, characterizing the occurrence of a so-called "critical state". The overall residual strength at the critical state isσ 12 /σ 22 ≈ 0.28. This value is only about half of the value of 0.6 adopted for the local particle contact friction µ, which can be explained from the fact that at the end of the deformation process the dilating particle structure, instead of sliding, is predominated by rolling of particles, see [7] for a more detailed discussion on this aspect.

Parameter
As a final note, it is mentioned that the contact moments in the cohesive packing determine about 10% of the total potential energy. The reason that this contribution does not become manifest in the stress ratioσ 12 /σ 22 depicted in Figure 11(a), is because the Cauchy stress is determined by contact forces only, see expressions (23) and (30). As already mentioned in Section 2.4, the overall effect by the contact moments on the particle aggregate can be accounted for through the introduction of a couple stress, but this stress measure is left out of consideration in the present study.

Conclusions
Novel numerical algorithms have been presented for the implementation of three types of classical boundary conditions for a particle aggregate. The micro-scale boundary conditions are formulated within the discrete element method using large deformation theory, and, along the lines of [16], are imposed on a frame of boundary particles of the particle packing, in accordance with i) a homogeneous deformation and zero particle rotation (D), ii) a periodic particle displacement and rotation (P), and iii) a uniform particle force and free particle rotation (T). The algorithms can be straightforwardly combined with commercial discrete element codes, thereby enabling the determination of the solution of boundary-value problems at the micro-scale only, or at multiple scales via a micro-to-macro coupling with a finite element model. The performance of the algorithms has been tested by means of discrete element method simulations on regular monodisperse packings and irregular polydisperse packings composed of frictional particles, which were subjected to various loading paths. The simulations provide responses with the typical stiff and soft bounds for the (D) and (T) boundary conditions, respectively, and illustrate for the (P) boundary condition a relatively fast convergence of the apparent macroscopic properties under an increasing packing size. Finally, a homogenization framework has been presented for the formulation of mixed (D)-(P)-(T) boundary conditions that satisfy the Hill-Mandel micro-heterogeneity condition on energy consistency at the micro-and macro-scales of the granular system. The numerical algorithm for mixed boundary conditions has been developed and tested for the case of an infinite layer subjected to a vertical compressive stress and a horizontal shear deformation, whereby the response computed for a layer of cohesive particles is compared against that for a layer of frictional particles. The results illustrate that the failure response for both contact laws is characterized by the development of a dilated particle structure, which at large deformation gradually turns into a critical state with an approximately constant residual strength and specific volume. The application of the present algorithms for multi-scale FEM-DEM analyses on granular systems with a large number of particles, and their extension towards a dynamics homogenization framework, are topics for future studies.