FEM × DEM : a new efficient multi-scale approach for geotechnical problems with strain localization

The paper presents a multi-scale modeling of Boundary Value Problem (BVP) approach involving cohesive-frictional granular materials in the FEM × DEM multi-scale framework. On the DEM side, a 3D model is defined based on the interactions of spherical particles. This DEM model is built through a numerical homogenization process applied to a Volume Element (VE). It is then paired with a Finite Element code. Using this numerical tool that combines two scales within the same framework, we conducted simulations of biaxial and pressuremeter tests on a cohesive-frictional granular medium. In these cases, it is known that strain localization does occur at the macroscopic level, but since FEMs suffer from severe mesh dependency as soon as shear band starts to develop, the second gradient regularization technique has been used. As a consequence, the objectivity of the computation with respect to mesh dependency is restored.


Introduction
Since its appearance in the late 70s, the Discrete Element Method (DEM) has become well known as an effective method for modeling granular material at the grain scale.Despite the advantages of the DEM approach, it seems extremely challenging in most cases to model directly with DEM boundary value problems (BVP) involving the actual size of both the grains and the structure, because of the tremendous amount of grains and current abilities of computers.
On the contrary, the Finite Element Method (FEM) is appropriate for the modeling of BVP.However, seeking for a constitutive model that can account for all the relevant properties of granular materials, or a significant part of them, is difficult -and even not possible.Such mathematical (often very sophisticated) models are necessarily phenomenological, and involve a number of parameters for which the physical meaning can be tricky.
To face this issue in the field of granular matters, the FEM × DEM multi-scale approach has been recently proposed by various researchers [1][2][3].Its basic concept lies on replacement of the phenomenological model by the constitutive response of the material with the use of DEM computations.
Using this approach, the actual grain scale can be considered even with large size macroscopic problems; without facing the intractable issue of dealing with trillions of grains in a DEM-mapped full field problem.Micro-scale related features, such as inherent and induced anisotropy of the granular material, or softening/hardening, emerge naturally from the micro-scale DEM model.
Localized failure is a major issue to be considered for modeling BVP related to geotechnical engineering.In this study, we present a (plane strain-FEM) × (3D-DEM) model of BVP with cohesive-frictional granular materials and an emphasis on strain localization.To avoid FEM mesh dependencies for shear bands, an internal length is introduced by means of the 2 nd gradient regularization technique.Thanks to this technique, associated with a high performance computing ability of the code (MPI parallelization), the tool is now able to address some real-size geotechnical structures by accounting grain-scale features and strain localization.
The paper is laid out as follows: Section 2 presents the principle of FEM × DEM method.Section 3 gives some results related to biaxial tests and hollow cylinder cases modeled by this innovative approach.Some conclusions are drawn in Section 4.

Multi-scale modeling
The resolution of a BVP by the FEM requires a constitutive law used at each integration point of the mesh.This law is supposed to express the local stress as a function of the history of the local deformation gradient.In our case, this constitutive law is not formulated mathematically according to a plasticity theory, nor it comes from a mathematically formulated phenomenological law; instead, it results directly from the response of a DEM volume element (VE) as a function of the history of local deformation gradient provided by the FEM (see figure 1).The latter response can be considered to be similar to that of a constitutive law with a large (but finite) number of internal variables such as: grain geometries, contact network and contact forces.The procedure to build this numerical homogenized law is detailed in [3,4].

3D discrete element micro-model
We describe only some main features of the micro model.The numerical model of granular behavior is herein obtained by a DEM approach using periodic boundary conditions (PBC) (see chapter 7 in [5]).The VE associated with each Gauss point is a packing of spheres.All spheres interact via a linear elastic force-law in the normal direction and the Coulomb friction law in the sliding direction.Accordingly, the normal repulsive contact force f el is related to the overlap distance as f el = k n • δ , where k n is a normal stiffness.The intensity of the friction force f t evolves incrementally and proportionally to the tangential relative displacement -with a tangential stiffness k t -up to a limit constant value ±μ f n when sliding occurs (μ is the Coulomb friction coefficient).In order to model cohesivefrictional granular materials, a local cohesion is introduced between the grains by adding an attractive (negative) force f c to f el .As suggested by [6], this local cohesive force can be set in regards to the mean level of reduced pressure , where a is the mean diameter of grains and σ 0 is the 3D isotropic pressure of VE.
All these features put together provide a numerical homogenized law (NHL) with a Mohr-Coulomb-like response, but it is actually much richer since it is able to handle very subtle responses involving, e.g., volume changes, non-associated flow, hardening/softening, fabric dependency and cyclic loading.

Coupling with the FE macro-model
The FEM×DEM approach has been implemented in the FEM code Lagamine [7], which is able to perform finite strain analysis.The implementation consists in pairing the DEM code as another "constitutive law", and solving some specific difficulties linked to convergence issues.The numerical integration of the constitutive law requires solving a nonlinear system of equations when the mechanical behavior is non-linear, which is the case for granular materials.In order to solve the non-linear system of equations, an incremental-iterative Newton-Raphson strategy is adopted [8].This requires computing a consistent tangent operator C, where the word "consistent" relates to the consistency with respect to the stress update algorithm as defined in Eq. 1 [9,10].) The consistent tangent operator can be derived analytically for a simple law, but for more complex laws, numerical differentiation needs to be adopted.However, a convergence problem is encountered when the consistent tangent operators are calculated in a phase of pronounced softening.Indeed, the scattered response provided by DEM calculation in the VE (and inherent to the nature of discrete matters) may cause some convergence difficulties.A specific modified Newton-Raphson algorithm has been developed in the code Lagamine to overcome these issues.Instead of the tangent operator defined above, we use an operator extracted from the underlying granular structure of the VE as described in [11,12].Since the relation by Kruyt and Rothenburg is only valid for 2D problem, an extended solution for the 3D case is proposed (Eq.2) .It results in the following form of elastic stiffness operator: where V is the volume of granular packing.k n , k t are normal and tangential contact stiffness, respectively; l is the branch vector connecting two centers of particles in contact; n, t and w are three vectors that define the contact forces directions: n is the normal direction, t and w are arbitrarily positioned in the sliding plan so that w = n × t.

Examples of BVPs
Although the NHL DEM model implemented in the study is a 3D model, it can be used in plane strain, plane stress and 3D FEM models.In this paper, we discuss plane strain BVPs, using 2D finite elements.

Biaxial Tests
Figure 2a presents the geometry and boundary conditions of biaxial test simulations.The constant confined pressures σ 2 is applied on the lateral sides and a vertical displacement at constant velocity is imposed on the top surface of the sample.
As shown in previous work of the same authors [3] or in ref. [13], it is well known that implementing strain softening constitutive laws in FEM produces mesh dependency: the deformation concentrates in zones as narrow as the mesh permits, independent of any material parameter.In order to restore a mesh independent behavior in such computations, an enrichment of the constitutive model following [13,14] is used.This enrichment involves the addition of a second gradient term to the constitutive model, the classical (first gradient) part being still given by the NHL DEM described above.The domain is discretized by using 128 (1 st case) or 512 (2 nd case) 2D elements second gradient (see figure 2b).These elements use eight nodes for the macro kinematic field (displacement field u i ), four nodes for the micro-kinematic field (micro deformation gradient v i j ) and one node for a Lagrange multiplier λ i j which is used to impose the equality of the microscopic and the macroscopic deformation gradient.Details about the second gradient method can be found in [13] or [14].
At the micro scale level, identical VE of 1000 spheres are used at each Gauss point.The volume is isotropically loaded and all contacts are compressive but with adhesion.This represents a cohesive-frictional granular material.The microscopic parameters have been chosen according to classical dimensional analysis: (i) the stiffness ratio is k n /k t = 1, (ii) the cohesion force f c is defined so that p * = 1, (iii) the normal stiffness k n is defined to obtain a stiffness parameter κ = k n /σ 0 = 500 and (iv) the intergranular coefficient of friction is μ = 0.5 The macroscopic responses of the biaxial tests are given in figure 3, with comparison to DEM simulation on the single VE of 1000 spheres.The responses show that the 2 nd gradient model regularize the post-peak responses, especially after ε 11 = 5% , as the softening response of  4 that uses the 2 nd invariant of strain tensor to exhibit the shear band at ε 11 = 10% .In first gradient model, shear bands tend to concentrate into a band of one element width [3].As result of 2 nd gradient model effect, shear band consists now of few elements.From this figure, it can be checked that the observed width of the shear band is intrinsic and independent of the mesh size.

Hollow Cylinder
These simulations predict the response of the material in the mid-place section of a pressuremeter, assuming plane strain condition in the axial direction.Because of the symmetry of the problem, only one fourth of a plane section is modeled.The domain is discretized using either first gradient element (Q8: 8 nodes, 4 Gauss points) or 2 nd gradient element (see figure 5d).The computation was performed as follows: starting from a homogeneous state of isotropic compression, the internal pressure (σ int ) was increased up to 4 times of the initial isotropic stress, while the external pressure (σ ext ) was kept constant.Numerical model setup and other related informations are given in figure 5.
Figure 6 shows the deformation mode in the model at the end of the loading in terms of the second invariant of strain tensor: (a) mesh with 952 standard Q8 elements (1 st gradient, i.e. no regularization); (b,c) enriched elements (2 nd gradient), two different meshes using 870 (b) or 1280 (c) elements.In the three cases presented in this figure, strain localization has taken place, organized in shear bands originated at the internal wall and progressing significantly inside the cylinder.This is the result of the inherent strain softening exhibited by the materials.In (b) and (c), the thickness of the shear bands does not depend on the mesh size, while in (a) it concentrates into a single layer of elements.These results demonstrate the capacity of this novel multi-scale approach to account for complex failure patterns.It is shown again that the second gradient method mitigates efficiently the mesh dependency problem observed with classical (first gradient) constitutive models.

Conclusion
A numerical law based on Discrete Element Modeling has been presented.This is accomplished in the framework of numerical homogenization applied to granular media.The model has been implemented in a finite element code.Due to scattering in the response of a model by DEM, some difficulties in the determination of consistent tangent operators occur, inducing convergence problems in the Newton-Raphson interactive processes in the FEM code.An alternative solution using an elastic stiffness operator instead of the tangent operator is shown to solve the difficulty.Thanks to the multi-scale FEM × DEM numerical approach, examples of BVPs involving strain localization is proposed, simulating biaxial tests and pressurized hollow cylinder cases.Strain localization is observed at macro level.The second gradient method is shown to solve efficiently the mesh dependency problem.

Figure 1 .
Figure 1.A sketch of the concept of DEM constitutive law

Figure 2 .
Figure 2. (a) Spatial discretization of the problem; (b) element enriched for the 2 nd gradient technique[15]