Geogrid pull-out modelling using DEM

With the advancement of the use of synthetic reinforcements in geotechnics, a greater understanding of the mechanisms involved in soil-reinforcement interaction is the focus of major research centres on the subject. The topic of this study is the shearing behaviour at interfaces between granular materials and geogrids. The main objective is to provide a more fundamental understanding of some micromechanisms present in this type of interface, which in turn are important to optimize the design of such reinforcement. The numerical modelling of these reinforced structures must deal with the complexity of the material-reinforcement interaction problem; therefore, it requires specific numerical models whose formulations admit localized behaviours in the contacts as well as the granular nature of the material (e.g., soil, gravel, ballast). A robust and flexible way of modelling this problem is through the Discrete Element Method (DEM). The DEM proposes to model this granular nature by representing the soil as interacting constituent particles, whose behaviour is ruled by physical laws defined at the contact points. The numerical approach is desirable since it allows, in an articulated and relatively fast way, studying closely different regions of the interface, in order to identify factors and variables that are important for the problem. The purpose involves the DEM for a 3D modelling of a geogrid pull-out test to calculate the magnitude of forces in different elements of the geogrid (i.e., nodes, longitudinal and transverse members). Preparation of numerical samples has a particular importance in the final results of simulations. Thus, the numerical techniques used to obtain better geometry for the geogrid and a granular assembly with a representative grain rolling effect are also presented in this paper.


Introduction
Discontinuous geomechanical problems with relevant deformation levels require specific numerical models whose formulations admit localized behaviours in the contacts. In fact, there is no numerical model capable of reproducing rigorously the soil behaviour, since they have complex arrangements and heterogeneous compositions. However, very realistic results can be obtained using models capable of incorporating the main soil characteristics required by the problem being studied.
The granular nature of the soil, as well as the interaction mechanisms, have a decisive impact on the mechanical behaviour of the soil layer. A powerful and flexible way of modelling this nature is by using the Discrete Element Method (DEM) [1], which makes it possible to characterize the soil as constituent particles with behaviours ruled by physical laws. This method is basically defined as a set of numerical processes for the calculation of motion and its respective effect on a given group of particles or elements.
This paper focuses on the material-reinforcement interaction, as well as the behaviour of synthetic reinforcements used in granular layers. Thus, the work uses the DEM for a three-dimensional modelling of a geogrid pullout test. The open-source framework Yade [2] was chosen * e-mail: marcus.moravia@3sr-grenoble.fr as the auxiliary tool for modelling. The discrete numerical approach allows for relatively easy access to information that is more difficult to obtain experimentally, such as the rotation mechanism of the granular particles and the deformation of the reinforcement members. As a contribution in this sense, the quantities measured and presented at this moment are the forces in the longitudinal and transverse members in different apertures of the geogrid. Although these are measures computed directly in synthetic reinforcement, they can be interpreted together with an expected soil-reinforcement interaction mechanism.
As in experimental work, the samples also need to be carefully assembled in numerical studies observing their state and incorporated properties. The preparation, therefore, involves assembly aspects of the samples as mechanical equilibrium and homogeneity. Because the reinforced layer is an interacting medium, its arrangement may interfere with the final results. Thus, preparation of the sample for the geogrid pull-out test needs special attention. As regards structural elements (i.e., synthetic reinforcements), their mechanical behaviour and the way they interact must be characteristic of a real element in order to obtain a more realistic model.
Finally, this paper is divided in three main parts. The first one summarizes the sample preparation and describes how the geogrid was modelled using DEM. The second part describes the results and the last one contains the final comments and conclusion.

Numerical sample
In DEM analyses, especially in quasi-static conditions that hold a strong influence of the initial state of the specimen, the sample preparation operation should be even more emphasized. As discussed by [3], there are numerical recipes to built a well-controlled, homogeneous and representative sample. These recipes involves observation of numerical and mechanical parameters considering the micro and macro domain of the sample, as well as the type of boundary condition (e.g., rigid walls and periodic cells). It is also important to note that in simulations involving structural elements, the mechanical behaviour of these elements and how they interact must be evaluated separately in order to ensure their correct use in the numerical model.
Aiming a geogrid pull-out test simulation, a packing of geogrid-reinforced granular assemblies characterizes the sample in this work. Therefore, it is important to obtain satisfactory conditions of homogeneity and equilibrium before starting the pull-out simulations. The following are the techniques used for the preparation of the numerical sample in question as well as the considerations regarding the mechanical behaviour of the geogrid. Note that the auxiliary tool for discrete modelling used here is the open-source framework Yade [2].

Granular assembly
The arrangement and interlocking of the granular assembly can significantly influence its mechanical behaviour. For example, [4] found greater shear strength in arrangements of elements with more angular or non-convex shapes. This fact is correlated with the rolling constraint provided by the geometric characteristics of the elements. The bending moment generated in the interaction between elements contributes to increase the overall resistance of the granular assembly. On the other hand, it is known that complex geometries can penalize the task of contact detection, which can lead to extremely time-consuming simulations.
Spheres compose the granular medium of the numerical model presented in this work. This choice was made to initiate interface measures from known geometries that allow a more explicit interpretation of the interaction mechanism. In addition, the results obtained with spheres can be used in comparisons with future analyses involving more complex samples with different geometric characteristics.
One of the characteristics of the DEM is the determination of motion from the calculation of the forces in the contacts. For this purpose it can be taken into account in each of the existing contacts friction force, elastic force, as well as a viscous damping force. In the case of the constitutive law used for the granular assembly of this work, only friction and elastic forces were considered. Figure 1 illustrates the scheme for the contact between two spheres.  1. Scheme of the contact between two spheres, adapted from [5] The law used for sphere-sphere interaction of the granular assembly (i.e., Law2_ScGeom_FrictPhys_CundallStrack [6]) computes contact forces and moments based respectively on relative displacement or overlapping with respect to the contact point and relative rotations. The contact normal force F n and contact incremental shear force dF s (i.e., in order to keep the shear component coherent during the movement of the interacting elements) can be expressed by the following equations: where k c n and k c s are the normal and tangential contact stiffness, respectively, u n is the relative displacement (i.e., overlap between spheres),u s is the relative shear velocity, and dt is the time step.
The shear component of the contact force F s incorporates in its calculation a Coulomb-type friction law without cohesion and an inter-particle friction angle φ. Thus, the shear force is limited by: Assuming both the normal and tangential contact as two springs arranged in series and with lengths equal to the respective radii of the spheres, the stiffness values in these points can be calculated by the harmonic mean of the individual rigidities of the springs. These contact stiffness values are defined as: where E 1 and E 2 are Young's moduli of each sphere in contact, R 1 and R 2 are the spheres radii, ν 1 and ν 2 are the Poisson's ratios. From this, the total force F in the contact can be applied to the interacting spheres in such a way that each particle acquires forces and moments for each step. Eq. (6) gives the expression for the moment of force in each sphere. These values are then used to integrate motion equations for each element separately, obtaining position and orientation data.
where τ is the moment of force, R i is the is the radius relative to the sphere, δ is the penetration depth of the contact, n is the normal unit vector (i.e., in the direction of the interaction axis), and F is the force applied at the contact.

Geogrid
Purely discrete analyses were chosen with the objective of reproducing the rolling mechanisms between the particles and the geosynthetic reinforcement during the pullout test. In this context, some authors (e.g., [7], [8], [9]) have modelled the reinforcement by means of clumps formed by rigid spheres connected together. In this case, the deformation behaviour of the geogrid is given by introducing internal degrees of freedom into the predefined sets of rigid spheres. As a drawback of using sphere clumps to represent the geosynthetic reinforcement is the significant increase in the number of elements, reducing the computational efficiency of the numerical model. In addition, an artificial numerical roughness is produced at the interface between the reinforcement and the granular material in contact.
Aiming at an approach that captures the continuous nature of the elements that compose the reinforcement but without the problems described above, the numerical representation of the geogrid was performed based on the method proposed by [10] for the discrete modelling of deformable objects with arbitrary geometries. The method is an extension of the technique introduced by [11] and later used by [12] in three-dimensional modelling of plant roots. It consists in applying the Minkowski sum [13], also known as morphological dilation, in polytopes and round bodies to describe the topology of objects. The idea is the representation of objects by a set of other geometric forms considered primitive as spheres, segments, walls, among others. From this approach, two types of elements were used for the geogrid discrete modelling, cylinder and particle facet elements.
Cylinder element is the result of the morphological dilation of a sphere and a segment. This object introduced in [12] behaves as a regular discrete element capable of interacting with other bodies. In the Yade code, cylinders are built by means of two nodes and a connection between them. The nodes are defined as centres of spheres and the connection generates the cylindrical shape. They are designed as non-flexible elements that can deform only in their longitudinal direction. Despite this fact, interconnected cylinders can be used to simulate materials with an elastic perfectly plastic behaviour as beams under tensile, compressive, shearing, bending or twisting loadings [12]. The constitutive relations are defined in terms of the nodes located at the element extremities called GridNodes.
In most cases, cylinders use a sphere-sphere law to compute interactions. The feasibility of using such a contact law lies in the trick of defining a virtual sphere within the cylinder with its centre in the projection of the contact point (i.e., contact between two interacting elements) on the cylinder axis. From these conditions, all virtual sphere status update (i.e., translational and rotational velocities) is interpolated and the efforts are then distributed on the nodes of the cylinder. It should be noted that displacements on the cylinder surface are calculated assuming that they vary linearly between the nodes at the element extremities. Figure 2 illustrates interactions involving virtual sphere in cylindrical elements. Therefore, the first class provides the law governing cylinder interaction with other elements and the second class the constitutive law for the cylinder mechanical behaviour. Law2_ScGridCoGeom_FrictPhys_CundallStrack is practically the same than Law2_ScGeom_FrictPhys_CundallStrack with the physical quantities divided and applied on the cylinder nodes (i.e., GridNodes). Considering the element mechanical behaviour, the following equations give the different values of stiffness used to calculate the internal forces and moments: where k n , k s , k r , k tw are the normal, tangential, bending, and torsional stiffness, respectively, E is Young's modulus of the cylinder, G is the cylinder shear modulus, R c is the radius of cylinder cross-section, and L c is the cylinder length. Law2_ScGeom6D_CohFrictPhys_CohesionMoment has limit conditions for the internal forces and moments so that it is possible to consider a plasticity behaviour of the cylindrical elements. The maximum tensile and shear forces are setted up to respective values of adhesion forces. If the normal or shear adhesion is achieved, the cylinder links break, which means that the adhesion values go to zero. In this situation, if a tensile force is still present, contact between nodes will be lost. The normal and shear adhesion forces are functions of the normal cohesion and shear cohesion parameters, respectively, which in turn can be used to change the element strength. Maximum values of bending and twisting moments are also defined. The threshold values for the moments are assumed by Eqs. (11) and (12): M max tw = 2 · σ max n · I tw R c (12) where M max r is the elastic limits for the bending moment, M max tw is the elastic limits for the twisting moment, σ max n is the maximum normal stress, I r is the bending moment of inertia (i.e., π · R c 4 /4), and I tw is the polar moment of inertia (i.e., π · R c 4 /2). In addition to the cylindrical elements, the numerical representation of the geogrid in this work also involved particle facet elements, so-called PFacet in Yade [10]. The use of this type of elements allows modelling deformable smooth objects in the sense that no interface roughness is introduced numerically. Aiming a representative geogrid shape, PFacet was used together with cylinders to model longitudinal and transverse members of the geogrid with flatter geometry.
The particle facet element is geometrically created by the Minkowski sum of a triangular facet and a sphere, resulting in three nodes (i.e., centres spheres) and three connections (i.e., cylinders). The constitutive behaviour of PFacet uses the informations stored in the nodes, applying the same idea used in cylindrical elements. [10] presented the methodology for the calculation of contact forces for the different types of interaction involving particle facet elements. The most relevant for the geogrid pull-out test simulation is the sphere-PFacet interaction. A virtual sphere is also introduced to calculate the PFacet interaction. It has the same radius as the cylindrical connections and its centre P is placed at the contact point projection on the plane containing the PFacet nodes. From the concept of barycentric coordinates, the following conditions are defined for the interaction occurring in the element facet area: where p 1 , p 2 , p 3 are the barycentric coordinates of the contact point projection on the plane of the PFacet nodes From these conditions, it is possible to note that p 3 is a function of p 1 and p 2 . The values of p 1 and p 2 are calculated as weighting functions for the corresponding PFacet nodes by the equations: where N 1 , N 2 , N 3 are the PFacet nodes, Ψ N 1 N 2 N 3 is the area of the triangle formed by the PFacet nodes, and Ψ PN 2 N 3 , Ψ N 1 PN 3 are the areas of two sub-triangles formed by P and PFacet nodes.
Yade uses the method detailed by [14] for a more efficient implementation of Eqs. (13) to (16). The scheme is similar to that of the virtual sphere in the cylindrical element. The virtual sphere displacement and rotation are linearly interpolated between the three PFacet nodes. Thus, the velocities of this sphere (i.e., translational and rotational) are calculated from the respective values obtained at the nodes N 1 , N 2 , and N 3 [i.e., Eqs. (17) and (18)]. The contact force F computed for the virtual sphere is then distributed on the PFacet nodes by Eq. (19).
where υ vs is the translational velocity of the PFacet virtual sphere, υ N i is the translational velocity of the node N i , ω vs is the rotational velocity of the PFacet virtual sphere, ω N i is the rotational velocity of the node N i , F N i is the distributed contact force at N i with i ∈ {1, 2, 3}, and F is the virtual sphere contact force. Figure 3 shows the final aspect of the modelled geogrid using cylinders and PFacets elements. Nodes, connections and facets can be seen in the enlarged detail. The flat shape of the numerical representation captures geometric attributes of real geogrids that are important for the rolling mechanisms present in the interaction of the reinforcement with granular media.

Sample preparation
It is important to mention that the geogrid pull-out test analyses in this work are not intended to reproduce a real or experimental case. At this time the central point is to increase understanding of the interaction behaviour between a granular medium and a generic geogrid. This study concerns the way the grains roll in the reinforcement and how the forces are mobilized locally in the longitudinal and transverse members of the geogrid.
The sample preparation consisted of obtaining different reinforced granular volumes for a given confinement prior to the geogrid pullout. The methodology for obtaining controlled compacted samples consisted basically of three main steps defined as generation of elements, application of the REFD (Radius Expansion -Friction Decrease) method, and preparation for the simulation.
As the name suggests, the first step is to generate the geogrid and the granular sample with predefined micromechanical properties in the sample volume. The geogrid is initially arranged fully stretched in the centre of the sample between layers of the granular material. A static condition is imposed on the geogrid before the next step, which is the use of the REFD method in the granular assembly. REFD [15,16] is a method for obtaining a dense assembly with specific porosity. It basically consists of increasing the elements radius by controlling the pressure in the contours of the sample volume up to a stipulated value, then the coefficient of friction is reduced and the elements radius increased in order to maintain the confinement condition until a desired porosity value. Note that there is no friction between the contours and the granular assembly. The third step involves, in this order, the recovery of the friction angle to the initial value, the removal of the static condition of the geogrid and the imposition of the gravitational force on the model, followed by the stabilization of the sample. Table 1 shows the parameters used for the samples.

Geogrid pull-out simulations
Three samples containing 3000, 6000, and 9000 spheres respectively were created under the same confining pressure and minimum porosity. Thus, the partial pull-out of the geogrid was performed by imposing constant displacement at the left end of the reinforcement (Figure 4). It is important to note that the model does not comprise interactions between the box and the geogrid, so that the only confinement exerted on the reinforcement is due to the granular layers. The first information obtained through the simulations was the force-displacement curve ( Figure 5). The mobilization of the initial resistance observed in the first peak of the graph is mainly caused by the friction between the granular material and the reinforcement and the secondary peaks are essentially due to the difficulty of movement imposed by the transversal members of the geogrid (i.e., the transverse members act as a buttress). The secondary peaks are more pronounced in the simulation with 3000 spheres. Smaller number of spheres results in larger diameters, consequently, the rolling constraint is greater in these larger mass spheres, resulting in greater pullout resistance. Another important point to note in these results is that the increase in the number of elements does not have a significant impact on the value of the initial peak resistance. Efforts in the regions of two geogrid openings were also tracked during pull-out. Figure 6 shows the position of the geogrid apertures analysed in particular. Thus, the axial stresses of the longitudinal and transverse elements identified respectively as A and C were plotted for cells 1 and 2 and shown in Figures 7 and 8.

Pull-out force (kN/m)
The measured values in the longitudinal elements of cell 1 are always larger than in cell 2. This behaviour is explained by the distribution of the anchoring forces, the greater efforts are concentrated in the region of application of the pull-out force and tend to zero at the contrary end of the geogrid. Although the trend of the values observed in cell 1 for the simulation with 3000 spheres is slightly lower than in the simulations with 6000 and 9000 spheres, the same behaviour can not be characterized in cell 2. Inversions of the values as well as exceptional peaks obtained between the simulations for the longitudinal elements can be explained by some specific arrangement of the spheres during geogrid pull-out, resulting in a greater confinement of the reinforcement element.
Differently from the longitudinal elements, in the transverse elements the values of the axial stresses in cells 1 and 2 follow the same range of values. This behaviour suggests that the main interaction mechanism of these ele-ments is similar to that of a buttress, restricting the spheres horizontal movement. The mean value of the stresses in the transverse elements is clearly lower compared to the longitudinal, where the resistive component due to the friction is larger. It is also important to note that the initial peak resistance in the transverse elements although lower is not negligible, with their mean value being approximately 14% of the value obtained for the longitudinal elements. In addition, secondary peaks can reach four times the calculated mean value.

Conclusion
The obtained results are coherent with the behaviour expected and presumed for particular regions of the geosynthetic reinforcement. In large deformations, the axial stress peaks in longitudinal and transverse members of the geogrid can reach four times the mean value. The stresses mobilized on the geogrid transverse members during the simulation are significantly smaller than those of the longitudinal members, but they are not negligible. Despite the determination of contact parameters still a difficulty for real cases, the discrete numerical model presented collaborates and shows great potential for future studies aiming at micromechanisms of interaction between granular materials and geogrid.