Daniel Dias

Hussein Mroueh

Márcio Almeida

Claire Silvani

Paulo Ivo

José Schiavon

Keywords: Geogrids, Discrete Element Method, Reinforcement Layers, Geogrid Pullout Test

for agreeing to be part of the examining committee, as well as for the fruitful discussion and their helpful comments.

Bruno Chareyre for clarifying critical issues in my work involving the YADE frame-

Résumé

Le travail de thèse porte sur la modélisation numérique par la méthode des éléments discrets (MED) du renforcement des sols granulaires par des géogrilles. L'objectif principal de ces travaux est d'améliorer la compréhension des mécanismes d'interaction entre la géogrille et les matériaux granulaires, en particulier ceux liés aux micro-mécanismes. En perspective, ces travaux peuvent contribuer au développement et à l'optimisation de nouveaux produits pour des applications spécifiques. Les géogrilles, en tant qu'éléments de renforcement des couches granulaires ont, par le passé, démontré leur efficacité à travers le monde. Dans de nombreux cas, ces structures apportent une solution le plus souvent économique, sûre et respectueuse de l'environnement. Cependant, le développement et l'optimisation de cette solution doivent tenir compte du comportement d'interaction complexe entre le renforcement et le matériau granulaire. Ainsi, l'utilisation de méthodes numériques capables de reproduire le comportement mécanique de la géogrille (étirement, flexion, cisaillement et torsion) et le comportement du sol à l'échelle du grain au voisinage de la géogrille (effet de roulement, de frottement et de butée) peut contribuer grandement à une meilleure compréhension de l'interaction sol-géogrille. Les modèles numériques basés sur la MED offre cette possibilité tout en ne nécessitant qu'un nombre relativement restreint de paramètres. Toutefois, la mise au point d'un modèle numérique de géogrilles représentatif des produits actuels reste une demande non satisfaite. Cette thèse comprend une nouvelle proposition de modélisation des géogrilles par la méthode des éléments discrets, basée sur l'utilisation d'éléments déformables. Les développements numériques ont été réalisés à partir du code open source YADE. Ils avaient pour objectif la définition d'une géométrie de géogrille générique constituée de brins de forme aplatie, comme c'est généralement le cas dans la pratique. Des simulations numériques d'essais d'extraction de géogrilles ont permis d'analyser le comportement d'interaction entre la grille et le sol granulaire et d'étudier la mobilisation progressive des forces de tension dans le renforcement. Des échantillons numériques comportant des particules sphériques ou non sphériques (clumps constitués de plusieurs sphères imbriquées) de différentes tailles ont été utilisés. Les résultats obtenus montrent que la forme des particules et leur taille relative (comparativement à la taille des ouvertures de la géogrille) ont une influence significative sur la mobilisation et l'intensité des forces d'extraction. Il a également été possible de vérifier que les efforts mobilisés dans les brins transversaux de la géogrille au cours des tests d'extraction sont nettement inférieurs à ceux des brins longitudinaux, mais qu'ils ne sont pas négligeables. Bien que la détermination des paramètres de contact reste une difficulté pour des applications à des cas réels, le modèle numérique discret semble être un outil pertinent et très prometteur pour l'étude de l'interaction entre tout type de géogrille et matériaux granulaires. 5.7 -Ratio of clump-particle sample to sphere-particle sample for mean values of stress-displacement curves in corresponding elements. . . . The ratio between the division of pullout force by the surface area of the geogrid and the vertical stress at the geogrid interface for the pullout test samples. . . . . . . . . . . . . . . . . . . . . . . . . . .

Mots clés :

List of Tables

Introduction 1.Topic and background

Deformations induced by loads can be considered one of the most important aspects of geotechnics. The good performance of geothecnical structures is directly related to their ability to withstand deformations produced by loads. In this context, reinforcements are commonly used, for example, in embankment foundations to increase their strength in order to avoid ruptures due to excessive deformation or shear mechanisms. The Figures 1.1 and 1.2 exemplify ultimate and serviceability limit states, respectively, which are normally considered in design of basal reinforced embankments. Among the synthetic reinforcements, geogrids have a successful application history throughout the world and are considered as a very interesting solution from the economic and environmental point of view. The mechanical improvement of a geogrid-reinforced granular layer is related not only to the stress-strain behavior of the reinforcement, but also to its grid-like geometry. Geogrid openings promote mechanisms of interaction with the soil or other granular material that help to restrict the deformation of the reinforced layers. However, there is still no complete knowledge of these mechanisms, especially micromechanisms, which are crucial for an optimized design of reinforced structures.

Structures reinforced with geogrid have complex behavior, and it is often difficult to experimentally study the interaction mechanisms between the soil particles and reinforcement members such as friction, grain rolling, or soil abutment at the vicinity of the geogrid due to stress concentration during pullout test, for example. On the other hand, the numerical approach represents a powerful alternative for the evaluation and development of this type of problem. In this case, numerical modeling must deal with the complexity of the material-reinforcement interaction problem, which requires specific numerical models whose formulations admit localized behaviors in the contacts as well as the granular nature of the material (e.g., soil, gravel, ballast). Thus, this thesis mainly concerns the numerical approach that can take into account the interactions of granular particles with geosynthetic reinforcement to model geogrid reinforced layers.

The granular nature of the material has a decisive impact on the mechanisms of interaction with the reinforcement and, consequently, on the mechanical behavior of layers reinforced with geogrid. A powerful and flexible way of modeling this nature is by using the Discrete Element Method -DEM (CUNDALL; STRACK, 1979), which makes it possible to characterize materials as constituent particles with behaviors 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. The use of this method by computer codes is a technique with great potential for problems involving par-ticle interaction. The continuous advancement of the processing technology has improved the use of this technique, allowing numerical models increasingly complex.

In many cases, reinforcement modeling is still challenging. The behavior of geosynthetics such as geotextiles or geogrids, which is related to their three-dimensional configuration, is not simple to be reproduced numerically, and simplifications may imply a less realistic mechanical response of the reinforced structure. Two modeling manner of this type of structure are predominantly found in the literature. One is the multi-domain coupling of different numerical methods, most commonly the DEM is coupled to the Finite Element Method -FEM [START_REF] Villard | Coupling between finite and discrete element methods for the modelling of earth structures reinforced by geosynthetic[END_REF][START_REF] Tran | A finite-discrete element framework for the 3D modeling of geogrid-soil interaction under pullout loading conditions[END_REF][START_REF] Tran | Three-dimensional analysis of geogrid-reinforced soil using a finite-discrete element framework[END_REF]. Thereby, the geotextile or geogrid are modeled as continuous solids discretized by finite elements, whereas the granular material is modeled by discrete elements. The other manner is a purely discrete modeling with the reinforcement defined by agglomerates composed of connected spherical elements [START_REF] Zhang | A few considerations of pullout test characteristics of geogrid reinforced sand using DEM analysis[END_REF][START_REF] Chen | Discrete element modelling of cyclic loads of geogrid-reinforced ballast under confined and unconfined conditions[END_REF][START_REF] Lackner | Microscopic interaction effects of prestressed geogrids in a reinforced soil element[END_REF].

Although the multi-domain coupling between DEM and FEM has advantages, the difficulty in precisely defining a contact behavior between the two methods is an inconvenience to the problem of geogrid interaction with granular materials. Another drawback is the additional parameters imposed on the model through interface elements, necessary to transmit the efforts between discrete and finite elements (DANG; MEGUID, 2013). It is also important to note that, depending on the model, the number of degrees of freedom for the geosynthetic or geogrid elements may be very large, resulting in high computational cost.

In the case of exclusively discrete models with geosynthetics or geogrids formed from agglomerates of spheres, the computational cost is also a difficulty, since the reinforcement modeling by this method introduces a considerable increase in the amount of elements of the problem. [START_REF] Chareyre | Dynamic spar elements and discrete element methods in two dimensions for the modeling of soil-inclusion problems[END_REF] also indicate two important limitations of this alternative. First, the artificial numerical roughness created on the surface of geosynthetic sheet or geogrid members by the spheres that compose it. The second is related to the deformation behavior of the reinforcement, if the axial strain in the members is large, significant voids will be created between the members spheres contacts. [START_REF] Effeindzourou | Modelling of deformable structures in the general framework of the discrete element method[END_REF] described a method for the discrete modeling of deformable objects with arbitrary geometries. The main idea is to apply the Minkowski sum [START_REF] Skiena | The algorithm design manual[END_REF], also known as morphological dilation, in polytopes and round bodies to describe the topology of objects. The method is an extension of the technique introduced by [START_REF] Chareyre | Dynamic spar elements and discrete element methods in two dimensions for the modeling of soil-inclusion problems[END_REF] and later used by [START_REF] Bourrier | Discrete modeling of granular soils reinforcement by plant roots[END_REF] in three-dimensional modeling of plant roots. With this idea it is possible to capture the continuous nature of the elements that compose the reinforcement but without the problems described above. Considering this, the present thesis presents a new proposal of geogrid modeling in discrete models, using two types of elements introduced in the works of [START_REF] Bourrier | Discrete modeling of granular soils reinforcement by plant roots[END_REF] and [START_REF] Effeindzourou | Modelling of deformable structures in the general framework of the discrete element method[END_REF]. This new proposal is used to compose the scope of the thesis, which is the study of aspects of the interaction between granular materials and reinforcement by three-dimensional DEM modeling of the geogrid-reinforced layer. The auxiliary tool for DEM used in this scope was the open-source framework Yet Another Dynamic Engine -YADE ( ŠMILAUER et al., 2015c).

Aim and objectives

The main motivation of this study is precisely the interface between reinforcement and granular materials. The lack of a better understanding of the stress distribution in the reinforced layer makes it difficult to accurately determine the material-reinforcement interface parameters in macroscopic models. Therefore, this negatively affects the analysis as well as the optimized design of geosynthetic reinforced structures. Thus, there is much interest in making progress in understanding this problem.

The interaction between reinforcement and granular material depends on many factors such as granular particle heterogeneity, reinforcement characteristics, confining pressure and porosity of the reinforced layer, reinforcement spacing in multilayer reinforcement solutions, among others. The numerical approach with a well articulated model allows agile studies on the specific importance of each of the factors involved and the identification of those with the greatest capacity to govern the behavior of the reinforced structure.

The aim of this work is to improve the understanding of the interaction between granular materials and geosynthetic reinforcement in geogrid-reinforced layers, through an approach by the Discrete Element Method. The perspective is that this thesis can contribute with new studies related to the use of geogrids or even help to skip steps in new works that choose to develop the approach presented here.

A new proposal for discrete geogrid modeling using deformable elements is presented as well as details for the preparation of numerical samples of geogrid-reinforced layer. In addition, geogrid pullout simulations are developed to study mechanisms of interaction between granular material and reinforcement. In the scope of this work, the following objectives can also be highlighted:

• To study possible numerical approaches to the geogrid-reinforced layer and to examine their respective strengths and weaknesses;

• To provide a numerical representation that captures the attributes of actual geogrids that are important to adequately describe their mechanical behavior (i.e., stretching, bending, shearing, and twisting) and their interaction with granular media particles;

• To develop a well-articulated numerical model that can take into account interaction between particles and geogrid elements;

• To analyze the behavior of the modeled geogrid for classical load cases and to compare with analytical solutions to validate the numerical model response;

• To discuss the possibilities that the proposed model offers to study the problem;

• To compare geogrid pullout simulations involving samples with different amounts of both spherical and agglomerated particles;

• To evaluate forces in longitudinal and transverse geogrid yarns as well as in different segments or members of the reinforcement structure by geogrid pullout simulation;

• To synthesize the individual importance of the parameters that govern the behavior of the geogrid-reinforced layer;

• To recommend improvements to the numerical model from its verified potential;

• To provide a study that can help improve numerical modeling of geogrid, which in turn contributes to solutions involving this type of reinforcement material.

Thesis outline

This thesis is divided into five main parts. The parts are separated into chapters as summarized below:

• Chapter 2 This section describes applications of geosynthetic reinforcements and some aspects considered important for the thesis topic. Notions of different approaches to modeling geogrid reinforced layers found in the literature are also described in this part. The chapter closes with a brief explanation of the Discrete Element Method.

• Chapter 3

The auxiliary tool for the Discrete Element Method used in this thesis was the open-source framework Yet Another Dynamic Engine -YADE. This chapter introduces this auxiliary tool, highlighting aspects of DEM in YADE that define the simulation loop. The first simulations performed in YADE in order to know important points for the development of the thesis numerical model are also included here.

• Chapter 4

This chapter is entirely concerned with the numerical model of the thesis. The modeling of both the granular assembly and geogrid used in sample preparation is detailed in this section. All processes and settings for the geogrid pullout simulation are also presented.

• Chapter 5

Simulations used to demonstrate and validate the numerical geogrid's response to classical load cases are presented. An overview of the capabilities of the developed numerical model is given prior to the presentation of geogrid-reinforced layer analyzes. Geogrid pullout simulations were performed to support these analyzes. The section concludes with the presentation and discussion of the results of three-dimensional geogrid pullout simulations for samples varying quantities and geometry of granular particles.

• Chapter 6

This last part of the thesis deals with the overall conclusion of the presented studies and final comments. A discussion about possible improvements of the developed approach and proposals for future works ends the present work.

2 Approach to the Problem

Geosynthetic reinforcement considerations

In the case of embankments on soft compressible or low strength soils, stability is mainly related to the shear strength of the foundation material, and the construction of these structures is generally conditioned to a load capacity problem, as noted in BSI (2010). Thus, solutions involving geosynthetic reinforcements (e.g., geotextiles and geogrids) for the purpose of improving load distribution in the compressible layer or uniforming the embankment foundation settlement are important alternatives to be rigorously developed.

Geosynthetic reinforcements with increasingly stiffness have been required to meet high embankments load levels, as analyzed by [START_REF] Blanco | Geosynthetic-reinforced embankments on soft soils: numerical analysis of the strain mobilization in the reinforcement at failure, service condition and staged construction[END_REF]. Geosynthetic reinforcements with high stiffness have restricted or even specific production and can result in economically uninteresting solutions. Therefore, the importance of reinforcement optimization in geosynthetic reinforced layer designs should be emphasized, which in turn depends mainly on understanding the interaction between materials.

It is worth mentioning that the geosynthetic reinforced layer technique does not only cover the use of coarse granular materials. For example, [START_REF] Riccio | Field monitoring and analyses of the response of a block-faced geogrid wall using fine-grained tropical soils[END_REF] confronts the measurement and prediction (i.e., by analytical and numerical methods) of reinforcement tension in an instrumented section of block-faced geogrid wall built using fine-grained tropical soils as backfill. There are several studies in the literature on the effect of geogrids as reinforcement in fine soils. [START_REF] Almeida | Embankment supported on piles with biaxial geogrids[END_REF] presented data from field monitoring of a pile-supported embankment reinforced with a bi-directional geogrid. In this case, the geogrid is between a pre-existing embankment and a backfill that comprises wellgraded sandy-clayey compacted soil. The investigations presented by [START_REF] Palmeira | Behaviour of geogrid reinforced abutments on soft soil[END_REF] showed the beneficial effect of the geogrid on the stability of two bridge abutments, including the reduction of lateral displacements. The filling material used to build the abutments was a random Brazilian soil with grain sizes varying from fine sand to blocks of rock. [START_REF] Palmeira | Uplift resistance of buried pipes anchored with geosynthetics[END_REF] and [START_REF] Tupa | Geosynthetic reinforcement for the reduction of the effects of explosions of internally pressurised buried pipes[END_REF] evaluated geogrids embed-ded in fine sandy soils. [START_REF] Palmeira | Uplift resistance of buried pipes anchored with geosynthetics[END_REF] examined the use of geogrid in buried pipes to increase maximum uplift resistance due to pipe bucking or buoyancy. [START_REF] Tupa | Geosynthetic reinforcement for the reduction of the effects of explosions of internally pressurised buried pipes[END_REF] investigated the soil-reinforcement technique to protect and minimize the consequences of eventual collapses in pressurized tubes.

In geosynthetic reinforced soil -GRS walls with well-graded sand backfill material, composed of crushed quartz powder with a significant amount of fines (i.e., 19 % passing the No. 200 sieve), [START_REF] Mirmoradi | Experimental evaluation of the effect of compaction near facing on the behavior of GRS walls[END_REF] conducted experimental studies on the effect and importance of compaction conditions close to the wall facing. Another example is the comparison presented by [START_REF] Santos | Performance of two geosynthetic reinforced walls with recycled construction waste backfill and constructed on collapsible ground[END_REF] between the performance of two GRS built with recycled construction waste backfill (i.e., approximately 12 % passing the No. 200 sieve).

Even in silty and clayey soils, the use of geogrids as reinforcement of the layer shows relevant improvements in mechanical behavior. For example, data from the instrumented pullout test of two types of geogrids shown by [START_REF] Sieira | Displacement and load transfer mechanisms of geogrids under pullout condition[END_REF] indicate this behavior for silty sand and clayey silt soils. [START_REF] Palmeira | Interaction between soils and geosynthetic layers in large-scale ramp tests[END_REF] presented theoretical and experimental studies of the interaction between soils and geosynthetic layers in large scale ramp tests. This study included a comparison between different types of soil and geosynthetic reinforcements. One of the soils used was a residual clayey soil with 70 % by weight passing through the No. 200 sieve. Among the reinforcements, there were two polyethylene terephthalate -PET geogrids with a square aperture width of 20 mm and 30 mm, respectively. França and Bueno (2011) performed creep tests using, among other geosynthetics, a geogrid with a square aperture width of 28 mm. These tests included the geogrid embedded in clayey sand soil. The results indicated that soil confinement (i.e., in-soil tests) considerably reduced the creep strain of the geosynthetic reinforcement.

In addition to the uses mentioned above, geosynthetic reinforcements also have applications in road and rail works. [START_REF] Garcez | Aplicação de geossintéticos como reforço de base em obras viárias[END_REF] describes the methods proposed by [START_REF] Giroud | Geotextile-reinforced unpaved road design[END_REF], Milligan et al. (1989a) e Milligan et al. (1989b), [START_REF] Sellmeijer | Design of geotextile reinforced paved roads and parking areas[END_REF], Houlsby and[START_REF] Houlsby | Design of reinforced unpaved roads for small rut depths[END_REF][START_REF] Burd | Analysis of membrane action in reinforced unpaved roads[END_REF] for the design of geosynthetic-reinforced bases. [START_REF] Ferreira | Avaliação estrutural de pavimento experimental reforçado com geogrelha[END_REF] discusses the benefits of this type of design solution and introduces newer design methods.

The geosynthetic reinforcement solution is also applicable in situations of shallow foundations supported by low bearing capacity, highly deformable, and low shear strength soils. [START_REF] Fabrin | uso de geossintéticos no reforço de fundações -aplicação ao caso de tubulações enterradas[END_REF] analyzed the applicability of geosynthetics to reinforce pipeline foundations in regions with low strength soils. The author indicated that the inclusion of the reinforcement reduces both the initial and consolidation strain as well as increasing the bearing capacity of the foundation layer.

It is important to note that simply the presence of basal reinforcement in embank-ments or foundation layers does not normally change the consolidation of the compressible layer. Nevertheless, the deformation behavior of the geotechnical structure is modified, particularly when the technique employed results in greater stiffness in the embankment base layer or foundation layer. For example, building a reinforced granular mattress prior to the embankment. Studies by [START_REF] Vidal | Simulation numérique non-linéaire de la construction de remblais en matériaux isotrope et anisotrope[END_REF], [START_REF] Williams | Design of reinforced embankments for great yarmouth by-pass[END_REF], and [START_REF] Dash | Model studies on circular footing supported on geocell reinforced sand underlain by soft clay[END_REF] have indicated that the stiffness properties provided by this layer at the base of the embankment influence the deformation of the low strength foundation soil, ensuring a more uniform load distribution and stress field for the foundation.

A more homogeneous stress field on the foundation layer is a desirable condition, especially when low strength soils are present, such as embankments on soft compressible soils. Therefore, the behavior for the reinforced granular mattress mentioned above suggests that including multiple geosynthetic reinforcement basal layers interspersed with granular material may be an advantageous solution. Multiple reinforcing basal layers, along with granular material, may result in a mattress layer with improved stiffness in addition to including tensile strength elements. In this sense, the use of more than one basal reinforcement layer with properties (i.e., stiffness and tensile strength) lower than those calculated considering a single layer may result in a more optimized and economically interesting solution.

Theoretically, it is possible to distribute the tensile force of a geosynthetic reinforcement between two or more layers of reinforcements overlapping each other. It should be noted, however, that in cases involving large strains and multiple reinforcements with different tensile strengths, field testing indicates that the higher strength reinforcement, generally the one placed at the base, mobilizes higher forces (ROWE; LI, 2001). Moreover, in multiple layers of overlapping reinforcements with identical properties, the installed lower layer tends to receive the highest tensile forces [START_REF] Koerner | Designing with geosynthetics[END_REF]. However, the precise distribution of forces in this type of structure is not fully understood. It depends on the interaction between geosynthetic layers and granular material. Many factors affect this behavior, for example, the attributes of the reinforcements, the distance between them, and the characteristics of the filler material employed.

One of the most recent publications on recommendations and guidance for multi-layer reinforcement problems is the standard BS 8006-1 (BSI, 2010). This standard prescribes, in multi-layer solutions, reinforcement layers of equal strength and stiffness and assigns individual reduction coefficients in the design strength calculation. Despite these considerations, BSI (2010) recommends, wherever possible, to use only one sufficiently strong reinforcement layer, justifying that the precise distribution of forces in such structures is not fully understood.

Indeed, the study of the behavior of structures involving a single reinforcement layer by FEM already presents difficulties in the correct determination of the interface parameters between the reinforcement and the granular material. To understand this subject matter, it is necessary first to explain the mechanisms that govern the material-reinforcement interaction. Therefore, the study of the behavior of geogrid-reinforced layers, showing the material-reinforcement interaction, is a relevant topic to be developed. Any contribution that enhances understanding of this interaction in both single-layer and multi-layer solutions is considered very useful not only in cases involving embankments on poorly supported soils, but also in cases of foundation reinforcement, reinforced soil wall, and road and rail infrastructure work. Knowledge gained from understanding the mechanisms involved in this problem contributes to the ability to evaluate stress concentration in both soil and reinforcement, improve numerical modeling, and develop better design solutions from a technical, economic, and environmental point of view.

Possible approaches to the problem

Although continuous numerical methods, for example, the Mohr-Coulomb linear elastic and perfectly plastic model of the classical soil mechanics, are widely used in engineering, there are problems where they are not suitable, such as those where discontinuities influence overall behavior. Discontinuous nature problems with significant levels of deformation require specific numerical models whose formulations admit localized behaviors in the contacts. Soil-reinforcement interaction has a complex behavior, especially under loading and unloading conditions, whose study requires the use of specific numerical models that can take into account the granular nature of the soil and its interaction with the reinforcement.

In reality, there is no numerical model capable of perfectly reproducing the behavior of soils with all their complex arrangements and heterogeneous compositions. However, the use of models that can take into account the main characteristics of soils or granular materials required by a given problem can lead to very realistic behaviors, providing proper results and assisting the development of geotechnical engineering.

The granular nature of the material and the mechanisms of interaction between the granular material and the reinforcement have a decisive impact on the mechanical behavior of the geogrid-reinforced layer. From this perspective, aiming at a numerical study that involves important aspects to reproduce a realistic behavior of the modeled structure, three approaches can be highlighted:

• Critical state-based advanced models;

• Cosserat continuum mechanics;

• Discrete Element Method.

Critical state-based advanced models calculate deformations of a given element from its stress state. These models may be valid for two-dimensional and three-dimensional stress states and are considered generalized (BRITTO; [START_REF] Britto | Critical state soil mechanics via finite elements[END_REF]. The Massachusetts Institute of Technology -MIT models MIT-E3 and MIT-S1 are relatively recent examples.

MIT-E3 was developed to describe the behavior of overconsolidated clays as well as clays under cyclic loading. This model aims to group important characteristics of the nonlinear, inelastic, and anisotropic behavior of the previously mentioned soil types (WHIT-TLE; KAVVADAS, 1994). [START_REF] Aubeny | Rotational interpretation of in-situ tests in cohesive soils[END_REF] obtained realistic estimates of shear stresses with the MIT-E3 model for an ideal cylindrical cavity expansion condition.

MIT-S1 is a generalized model capable of simulating elastoplastic properties of materials and was developed to analyze the anisotropic behavior of a wide variety of soils. The main idea that distinguishes it from other models (i.e., constitutive laws) is the explicit inclusion of effective stress and void ratio as independent variables controlling the mechanical response of the model, including anisotropy directions based on the orientation of the bounding surface and its evolution with rotational hardening. In this type of approach, new expressions can be introduced to describe nonlinearities in small deformations, which allows more realistic evaluations in both clay and sand, as described in [START_REF] Pestana | Formulation of a unified constitutive model for clays and sands[END_REF].

In the context of generalized continuous media, the Cosserat theory or micropolar elasticity is appropriate for modeling rupture mechanisms in materials with strong microstructure influence on their overall behavior. Cosserat continuum adds rotational degrees of freedom to the conventional continuum, which makes it possible to take into account moments at any point in a given material, in addition to the usual stress field (COSSERAT; COSSERAT, 1909). Each particle is comparable to a small-size continuum around a point that characterizes it. Because of this, in the kinematics of the material point, additional rotational degrees of freedom occurs, and the gradient of these rotations is associated with a stress moment tensor by the principle of virtual works [START_REF] Figueiredo | Analysis of bookshelf mechanisms using the mechanics of Cosserat generalized continua[END_REF]. Another aspect is the constitutive description of the material that considers intrinsic lengths. The introduction of these lengths in the constitutive relations allows indirect consideration of the particle size and geometry in the macroscopic behavior of the medium (VARDOULAKIS; SULEM, 1995; VARDOULAKIS, 2019). [START_REF] Ebrahimian | Modeling shear localization along granular soil-structure interfaces using elasto-plastic Cosserat continuum[END_REF] simulated shearing at the interface between a layer of granular material and a rough structure using an improved model with Cosserat theory. Although the lack of both numerical and experimental studies that investigate deformations in the soil-structure interface regions makes it difficult to validate this model, the advances obtained by the author indicated the relevance of this approach. An interesting study was the comparison between parameters that control the thickness of the formed shear band. Between the initial void ratio, average grain size, and vertical pressure, the first two indicated greater influence on the shear band. It is worth mentioning that the authors assumed simplifications that affect the model's behavior. For example, balancing factors, which reflect aspects of the microstructure in the constitutive model, such as slipping and shearing between particles, were defined as equal to one by simplification. The work also assumed full bonding in the interface region, so there are no relative horizontal displacements between the bottom surface of the granular layer and the rough surface. Regarding the kinematic boundary conditions in the interface, the study considered two antagonistic cases. The first one considers the model with zero couple stresses, that is, the free Cosserat rotation is assumed, which is a condition equivalent to that of classical continuum mechanics. The second case assumes fully constrained Cosserat rotation (i.e., zero rotations). The two cases are antagonistic, so a study that includes an intermediate boundary condition to the described cases could represent a behavior closer to that expected in geotechnical soil-structure interaction problems. Simplifications are usually necessary to introduce very complex problems, and considering the lack of research on the subject, the authors' numerical results contribute to further advances in this topic.

The Discrete Element Method, also called the Distinct Element Method, models the material as constituent particles in which contacts can change during deformation. This method is essentially a set of numerical processes for calculating motion and its effect on a particular group of particles or elements whose behavior is governed by physical laws. [START_REF] Cundall | A computer model for simulating progressive large scale movements in blocky rock systems[END_REF], as cited in [START_REF] Cundall | A discrete numerical model for granular assemblies[END_REF], originally proposed this method to study problems related to Rock Mechanics. However, its application has extended to the study of micromechanisms in granular media, as in [START_REF] Cundall | A discrete numerical model for granular assemblies[END_REF], [START_REF] Hori | Micromechanical analyses on granular column formation and macroscopic deformation[END_REF], and [START_REF] Mirghasemi | Numerical simulation of assemblies of two-dimensional polygon-shaped particles and effects of confining pressure on shear strength[END_REF], and in the behavior of clays, as in [START_REF] Anandarajah | Discrete element modeling of leaching-induced apparent overconsolidation in kaolinite[END_REF]. The simulation of a medium by a discrete body system (i.e., elements) is a highly dynamic process with periodic changes of the forces acting on the contacts. Thus, the computational implementation of DEM is fundamental and also makes possible the use of this numerical method in different types of problems. Consequently, the computational cost of a given problem becomes a crucial point for the discrete numerical approach.

A moderate increase in the number of elements for a more realistic numerical model can produce a significant increase in computational tasks. There is a high computational cost inserted in the task of updating the contact status between elements since the model performs this task whenever there are body displacements that, in turn, change the contact status. Many DEM codes take advantage of parallel processing capabilities (i.e., a particular coupled form of distributed computing) to extend the number of elements in a simulation. With the enhancement of computer processing capabilities and the use of more agile numerical algorithms, the simulation of increasingly complex problems by DEM has become feasible and accessible.

Unlike critical state-based advanced models and Cosserat continuum mechanics ap-proaches, DEM allows modeling the complex behavior of granular materials and their interaction with structural elements by using a relatively small number of parameters. For this reason and due to the discrete feature of the method, which provides effective modeling of the granular nature present in geogrid-reinforced layers and also the interaction mechanisms between material and reinforcement, this thesis focuses on the use of DEM approach. The tool chosen to support the discrete modeling was YADE. It is an extensible open-source framework for discrete numerical models focusing on DEM. The use of this tool is mainly due to the robustness of the codes and the possibility of implementing algorithms and functions with full source code control.

Discrete Element Method -DEM

A dry granular medium, for example, is composed of a large number of particles that can move separately and interact with each other at the contact points. This discrete aspect produces complex media behavior under loading and unloading conditions. Discrete modeling makes it feasible to study this behavior in an articulated manner. As mentioned in the previous section, the method is composed of numerical processes that calculate motion and its effect on a group of particles or elements.

In a way, DEM can be considered similar to Molecular Dynamics, which studies the physical motion of atoms and molecules from the interaction potential between particles and the equations governing their movement [START_REF] Haile | Molecular dynamics simulation: elementary methods[END_REF], but including rotational degrees of freedom, contact state, and more complex geometries (e.g., polyhedra). The first studies involved applying the method to the simulation of progressive movements in rock masses (CUNDALL, 1971 apud CUNDALL;[START_REF] Cundall | A discrete numerical model for granular assemblies[END_REF]. Many studies today use DEM to model a wide variety of materials and applications have increasingly shown that the method is a helpful, powerful, and necessary tool for geotechnical problems involving discontinuous media.

Despite the application of the method to problems involving mechanical interaction between particles, DEM can also incorporate physical-chemical interaction laws for analysis in clays, for example. [START_REF] Anandarajah | Discrete-element method for simulating behavior of cohesive soil[END_REF] developed an approach using DEM to simulate cohesive soil behavior. This study ignored the attractive force but modeled the repulsive force based on a study on the double-layer repulsive force between inclined particles. Other examples are the three-dimensional DEM model developed by [START_REF] Yao | Three-dimensional discrete element method of analysis of clays[END_REF], which comprises cuboid particles for numerical analysis of assemblies of clay particles, and the DEM framework developed by [START_REF] Lu | A discrete element model for kaolinite aggregate formation during sedimentation[END_REF] to quantitatively model the formation of kaolinite aggregates during a sedimentation process. The possibilities offered by DEM make it capable and with the potential to also deal with problems in tropical soils with fines.

According to [START_REF] Cundall | A discrete numerical model for granular assemblies[END_REF], the basic principle for the formulation of the method is the idea of a dynamic process based on the laws of motion of the English physicist Sir Isaac Newton, in which, from the propagation of perturbations applied to the limits or the elements themselves, individual movements in the particles are generated, which in turn result in contact forces and subsequent displacements until an equilibrium condition. In the numerical description of this dynamic process, accelerations and velocities are considered constant at the specified time intervals. The assumed time interval is so short that the above consideration is acceptable. The method also assumes that over a time interval, applied or generated perturbations can only propagate to immediate neighboring elements. Therefore, at each time step, the resultant force on any element are calculated exclusively by the sum of the contact forces and the field forces imposed on it, such as gravitational, magnetic, and electrostatic forces [START_REF] Bharadwaj | Using DEM to solve bulk material handling problems[END_REF]. The flowchart of Figure 2.1 presents the basic scheme of a DEM simulation, which is characterized by the balance of forces on the elements at each time increment to reproduce the complex behavior of the material. The DEM calculation steps alternate between applying a force-displacement law and Newton's second law of motion. The first law provides the contact forces resulting from displacements caused and the second law the movement of each element from the forces acting on it. To illustrate the calculation cycle and to show how the method calculates forces and displacements, Figure 2.2 replicates the example presented by [START_REF] Cundall | A discrete numerical model for granular assemblies[END_REF]. The simplicity of this case makes the explanation of the method didactic. It is important to note that since the strain of a single particle is considerably lower than the whole sample strain, it is not necessary to accurately determine the deformations of each element to obtain adequate mechanical behavior of the numerical sample. Thus, in the formulation presented by [START_REF] Cundall | A discrete numerical model for granular assemblies[END_REF], the overlapping elements indicate particle strain during contact. Because of that, at time t1, the elements overlap at the contact points A and C, as shown in Figure 2.2.

In Figure 2.2, after the time increments, the contact point A is defined as half the distance between the projected points on the disk and the wall and defined by A D and A W , respectively. Thus, the normal relative displacement at contact and time t1 is the displacement of point A W relative to that of point A D and can be calculated by:

∆n A t1 = v • ∆t (2.1)
where ∆n A t1 = relative displacement at contact point A occurring at time t1; v = wall velocity; ∆t = time increment.

Likewise, it is possible to calculate the relative displacement at point C and the re-spective contact forces. Then, a force-displacement law applies, which can be as follows:

∆F n = k n • ∆n t1 = k n • v • ∆t (2.2)
where ∆F n = increment in normal force;

k n = normal stiffness coefficient; ∆n t1 = normal relative displacement at contact point at time t1 (i.e., v • ∆t).

Taking positive forces in the directions of the axes shown in Figure 2.2 and knowing that the contact forces are zero at time t0, the resultant forces for disks x and y at time t1 are calculated by:

F x t1 = k n • ∆n t1 (2.3) F y t1 = -k n • ∆n t1 (2.4)
where F x t1 = resultant force on disk x at time t1;

F y t1 = resultant force on disk y at time t1.

From the contact forces on each element, it is possible to calculate the acceleration generated on the disks at time t1 by applying Newton's second law (i.e., Eqs. 2.5 and 2.6). ẍ1 = F x t1 /m x (2.5)

ÿ1 = F y t1 /m y (2.6)
where ẍ1 = acceleration of disk x in the 1-axis direction; ÿ1 = acceleration of disk y in the 1-axis direction;

m x = disk x mass; m y = disk y mass.
The calculated accelerations in Eqs. 2.5 and 2.6 are constant over time between t1 and t2 for very short time intervals and can be integrated to obtain velocities (i.e., Eqs. 2.7 and 2.8). ẋt2 = (F x t1 /m x ) • ∆t (2.7) ẏt2 = (F y t1 /m y ) • ∆t (2.8) where ẋt2 = disk x velocity in the time interval between t1 and t2; ẏt2 = disk y velocity in the time interval between t1 and t2.

Thus, the relative displacement increments at the contacts (i.e., contacts points A, B, and C) at time t2 are calculated Eqs. 2.9 to 2.11, considering positive displacements for compression:

∆n

A t2 = [v -(F x t1 /m x ) • ∆t] • ∆t (2.9) ∆n B t2 = [(F x t1 /m x ) • ∆t -(F y t1 /m y ) • ∆t] • ∆t (2.10) ∆n C t2 = [(F y t1 /m y ) • ∆t -(-v)] • ∆t (2.11)
where ∆n A t2 = relative displacement increment at contact A at time t1;

∆n B t2 = relative displacement increment at contact B at time t1;

∆n C t2 = relative displacement increment at contact C at time t1.

In summary, the basic DEM calculation cycle involves firstly calculating the contact forces caused by some disturbance in the elements, secondly the accelerations associated with these forces, thirdly the velocities as a function of accelerations and finally the resulting displacements. The integration scheme divides the time step into two parts so that it is possible to obtain forces and accelerations for the first half and the velocities and displacements for the end of the time step. This cycle composed of a force-displacement law and Newton's second law repeats until an imposed condition or time limit, usually representative of a state of equilibrium.

In the general case, involving many elements and three-dimensional configuration, the force-displacement law applies to each of the contacts, and both normal and shear components can result. It is possible to incorporate a Coulomb-type friction law into the formulation so that the model provides shear strength in the contacts. In this case, the vector sum of contact forces of an element calculates the resultant force acting on it. It is also possible to introduce contact damping forces into the formulations. This damping enters as a value in the sum of the contact forces. The method is articulated and also allows to include other types of force components in the model, for example, global damping in motion calculations, acting directly on the absolute velocities of the elements.

Discrete modeling of geosynthetic reinforcements

The soil-structure interaction in problems involving large deformations is one of the appealing and challenging topics of geotechnical engineering. This topic is broad and can include tunneling, foundations, retaining structures, soil reinforcement, etc. In cases involving geogrid-reinforced granular materials, FEM-only numerical study cannot adequately capture the interlocking of granular materials in the openings of the geosynthetic reinforcement and therefore does not consider proper interface features of such structures. In contrast, DEM has proven its efficiency in modeling micromechanical problems. In this context, some authors present alternatives involving discrete models for the study of this interaction of granular materials with the reinforcement layer.

One approach used is the FEM and DEM coupling (O ÑATE; ROJEK, 2004) to take advantage of both methods, modeling in one single model the continuous behavior of structural elements by FEM and the discrete nature of granular materials by DEM. Although this multi-domain approach is already established, including open-source availability for it (STR ÁNSK Ý; JIR ÁSEK, 2012; STR ÁNSK Ý, 2013), there are still few studies focusing on modeling geosynthetic reinforcements (e.g., geotextiles and geogrids). Works such as those presented by [START_REF] Villard | Coupling between finite and discrete element methods for the modelling of earth structures reinforced by geosynthetic[END_REF], [START_REF] Tran | A finite-discrete element framework for the 3D modeling of geogrid-soil interaction under pullout loading conditions[END_REF][START_REF] Tran | Three-dimensional analysis of geogrid-reinforced soil using a finite-discrete element framework[END_REF] still seek to improve understanding of interactions between granular materials and reinforcement elements. [START_REF] Villard | Coupling between finite and discrete element methods for the modelling of earth structures reinforced by geosynthetic[END_REF] presented a dedicated numerical model to study soil interaction with a continuous and homogeneous geosynthetic reinforcement layer. The authors proposed a specific contact law for the soil-geosynthetic interface considering the case of a reinforced embankment over an area subject to potential cavities (i.e., subsidence). [START_REF] Tran | A finite-discrete element framework for the 3D modeling of geogrid-soil interaction under pullout loading conditions[END_REF] investigated the pullout behavior of a biaxial geogrid in a granular medium. [START_REF] Tran | Three-dimensional analysis of geogrid-reinforced soil using a finite-discrete element framework[END_REF] proposed multi-domain coupling to numerically investigate two cases, one related to the efficiency of geosynthetic reinforcement in increasing soil loading capacity for shallow foundations and the other concerning the use of geogrid in a hypothetical subsidence case. [START_REF] Tran | A coupled finite-discrete element framework for soil-structure interaction analysis[END_REF] presents the mathematical and numerical development of a FEM and DEM coupling algorithm for the study of soil-structure interactions. The framework used by [START_REF] Tran | A coupled finite-discrete element framework for soil-structure interaction analysis[END_REF] and [START_REF] Tran | Three-dimensional analysis of geogrid-reinforced soil using a finite-discrete element framework[END_REF] to couple the two numerical methods is an extension of the work originally presented by [START_REF] Dang | An efficient finite-discrete element method for quasi-static nonlinear soil-structure interaction problems[END_REF]. These works used the dynamic relaxation method presented by [START_REF] Dang | Evaluating the performance of an explicit dynamic relaxation technique in analyzing non-linear geotechnical engineering problems[END_REF] in the part that applies FEM.

It is important to highlight that although multi-domain coupling is an interesting alternative to multi-scale problems, as presented by [START_REF] Cheng | A simple multiscale model for granular soils with geosynthetic inclusion[END_REF], in this approach it is difficult to precisely define the contact behavior between discrete elements and finite ele-ments, especially in cases involving complex three-dimensional geometries. It is necessary to use interface elements to make this contact between numerical methods, which in turn impose new parameters on the model. Moreover, when the FEM domain is composed of volumetric elements, the number of degrees of freedom of this domain can be very large, resulting in high computational cost simulations.

Another approach to the problem found in the literature is purely discrete models. Authors have modeled geogrids from rigid aggregates or agglomerates formed by less complex geometry elements linked together, such as spheres. Figure 2.3 illustrates a numerical representation of a geogrid by this approach. The introduction of internal degrees of freedom in the rigid aggregate of spheres allows the modeled structures to be deformable. Purely discrete modeling applications involve diverse studies. [START_REF] Konietzky | Use of DEM to model the interlocking effect of geogrids under static and cyclic loading[END_REF] studied the effect of geogrid interlocking for static and cyclic loading cases. [START_REF] Mcdowell | Discrete element modelling of geogrid-reinforced aggregates[END_REF] applied discrete modeling for a geogrid-reinforced ballast layer. [START_REF] Zhang | A few considerations of pullout test characteristics of geogrid reinforced sand using DEM analysis[END_REF] studied the influence of soil compaction, considering different porosity values, in geogrid pullout simulations. [START_REF] Lobo-Guerrero | Fibre-reinforcement of granular materials: DEM visualisation and analysis[END_REF] analyzed the mechanical improvement of granular materials with the inclusion of synthetic fibers. [START_REF] Bhandari | Investigation of geotextile-soil interaction under a cyclic vertical load using the discrete element method[END_REF] investigated soil-geotextile interaction in cases of cyclic vertical loading. [START_REF] Chen | Discrete element modelling of cyclic loads of geogrid-reinforced ballast under confined and unconfined conditions[END_REF] compared geogrid-reinforced ballast layer behavior under confined and unconfined conditions. [START_REF] Ferellec | Modelling of ballast-geogrid interaction using the discrete-element method[END_REF] examined the influence of ballast particle geometry on ballast-geogrid interaction. [START_REF] Han | DEM analysis of stresses and deformations of geogrid-reinforced embankments over piles[END_REF] studied the behavior of the geogrid in solutions involving piled embankment with basal reinforcement. [START_REF] Chen | A study of geogrid-reinforced ballast using laboratory pull-out tests and discrete element modelling[END_REF] simulated the pullout test of geogrid reinforcement in layers of railway ballast. [START_REF] Lackner | Microscopic interaction effects of prestressed geogrids in a reinforced soil element[END_REF] studied mechanisms of interaction between soil particles and geosynthetic reinforcements.

As an example of more recent studies, [START_REF] Chen | Discrete element modelling of geogrids with square and triangular apertures[END_REF] presented two three-dimen-sional models of geogrids with square and triangular openings aiming at a realistic reinforcement shape. Another example is the work in which Chen et al. (2019) performed two-dimensional numerical pullout tests with discretely modeled geogrids and evaluated the effect of geogrid tensile stiffness on the micro-mechanical behavior of the reinforced layer. Both cases model the geogrids connecting spherical elements. They use the parallel bond contact model to create different sphere arrangements. This contact model is like two parallel surfaces positioned in the contact plane that provides mechanical properties of an elastic bond between the two contacting spheres. Thus, it is possible to take into account the transmission of forces and moments between the particles that make up the geogrid.

It is important to highlight two downsides of this type of numerical approach. The first is that modeling the geogrid or geosynthetic reinforcement from sphere aggregates, for example, results in a significant increase in the total number of elements considered, which reduces the computational efficiency of the numerical model. The other one is an artificial numerical roughness on the reinforcement surface caused by the agglomerated elements, which may result in unrealistic behavior of the numerical model. The use of deformable elements indicates a refined way of numerically simulating the geogrid, capturing the continuous nature of the elements that make up this type of reinforcement. Considering the potential that this approach has in modeling complex deformable structures interacting with other elements and its relevance to problems involving geogrid-reinforced layer behavior, further studies based on this technique are devoted to this thesis.

Pullout testing of geogrid reinforcements

This thesis includes geogrid pullout simulations to study aspects of the interaction between granular material and reinforcement. Despite this, the proposal is a numerical model for the geogrid-reinforced layer that is versatile enough to be suitable for other types of tests or different loading conditions. This section highlights the general characteristics of force-displacement curves obtained by pullout tests (i.e., laboratory data found in the literature). The intention is to present the expected shape and general tendencies for this type of curve. [START_REF] Palmeria | Scale and other factors affecting the results of pull-out tests of grids buried in sand[END_REF] grouped several results of the grid pullout tests to investigate factors that affect the results. Figure 2.5 compares the effect of different top boundary conditions. In this case, one test comprised a rigid rough plate to apply vertical pressure to the sample, while the other used a flexible bag filled with water. The discussion presented by [START_REF] Palmeira | The study of soil-reinforcement interaction by means of large scale laboratory tests[END_REF] shows the influence of the initial distance from the geogrid to the box wall on the pullout test results (i.e., Figure 2.6). In the literature, there are several studies on other effects and details to be observed to minimize them. However, the force-displacement curves of pullout tests usually show common tendencies caused by the load transfer mechanisms between the soil and reinforcement. The pattern follows an increase in pullout force to a peak value as the reinforcement starts to move inside the sample. After the peak, pullout forces tend to decrease because the length of the reinforcement inside the sample also decreases. [START_REF] Palmeira | Bearing force mobilisation in pull-out tests on geogrids[END_REF] evaluated the mobilization of forces on the grids using pullout tests. The test results comprising different grids indicated that those with smaller aperture size and greater tensile stiffness result in curves with more noticeable peaks (i.e., Figure 2.7). Sample confinement can also have a significant effect on the pullout resistance of the reinforcement. Although the increase in this condition reduces soil tendency to dilate, higher confining pressures also lead to an increase in the passive soil resistance on the transverse ribs; consequently, the geogrid pullout forces increase [START_REF] Farrag | Pull-out resistance of geogrid reinforcements[END_REF]. Figure 2.8 shows the results of a pullout experiment comprising two confining conditions of the sample. Although many factors influence the results of pullout tests, the force-displacement curves of geogrids show common behaviors. The curves presented in this section can be a reference for force-displacement curves obtained by the geogrid pullout simulations of this thesis. They are useful to show that the curves obtained by the present work have similar tendencies (i.e., concerning general behavior).

Yet Another Dynamic Engine -YADE

YADE is an extensible open-source framework for discrete numerical models. It has an environment designed for scientists or other users to implement computational algorithms and allows for easy reuse, exchange, and enhancement of code while preserving, developing, and disseminating acquired knowledge (KOZICKI; DONZ É, 2009). The project began in 2004 at the Laboratoire 3SR of the Université Grenoble Alpes, France, and derived from the SDEC code initially developed by Frédéric Victor Donzé and introduced by Donzé and Magnier (1995). Many research institutes currently promote the continued development of this tool that has an active and productive community.

The conception of YADE involves a layered environment. These layers consist of structured and ordered parts with some of them having interdependence (e.g., an upper layer depending on a lower layer). The innermost or lower layer is the library layer. Its function is to make functions and methods available to other parts of the code. Above this layer, there is a generic layer representing the YADE core that is responsible for providing the interfaces for simulations (e.g., functional or engine components, discrete bodies, and modes of interaction between elements). One level up is the common layer which, as its name implies, contains components commonly used by different types of simulation, including for simulations coupled with other methods such as the Finite Element Method -FEM [START_REF] Zienkiewicz | The finite element method: its basis and fundamentals[END_REF] and Smoothed Particle Hydrodynamics -SPH (LIU; LIU, 2003). Examples of these components are Newton's and Hooke's laws, time-domain numerical integration algorithms [e.g., Runge-Kutta 4th order method (BUTCHER, 2008), Newmark's method (NEWMARK, 1959), and leapfrog method (FINCHAM, 1992)], damping methods [e.g. Cundall's non-viscous damping (CUNDALL; STRACK, 1979)], collision detection algorithms, boundary conditions (e.g., displacement imposition and gravitational force), data classes (i.e., information about elements and interactions), methods for drawing common geometries (e.g., Open Graphics Library -OpenGL). At the next level is the specialized layer that contains code packages specific to numerical methods. There are many distinct packages, such as those defined for DEM, FEM, SPH, and lattice geomet-rical model [START_REF] Kozicki | Application of discrete models to describe the fracture process in brittle materials[END_REF]. Finally, the top layer is a graphical user interface. Kozicki and Donzé (2009) gives more details about YADE's layer structuring.

YADE's open-source code1 is well structured in C ++ programming language, which makes it flexible to adaptations or even additions of specific functionality by implementing independently new algorithms. Python programming language is used directly for simulations and allows the development of scripts for models and mechanisms, debugging process, manipulation and control of procedures and data, and post-processing of simulation results. Summarily, discrete bodies and their interactions constitute running simulations. Regarding interactions, algorithms and functions based on physical laws, so-called engines, govern them. Simulations can be defined as iteration processes that typically run more than once or in time steps. Iteration processes are fundamentally cycle of actions that are defined by engines and functions. The sequence of a typical DEM simulation involves the reinitialization of the forces acting on the discrete bodies (i.e., only for the second cycle onwards), identification of potential interactions, collision type detection (i.e., classes of interacting elements), accurate collision detection, imposition of conditions on the model (e.g., gravity), computation of forces on interacting elements, and the modification of the positions of the elements by integrating the equations of motion, depending on the forces applied on them. Before each time increment that characterizes the last action of the simulation cycle, it is possible to include additional algorithms for the model to perform desired actions(e.g., data mining functions and tasks for coupling between numerical methods).

As described earlier, YADE uses C ++ that is an object-oriented programming language. Thus, the code structure benefits from the conception of classes and objects in programming. In the simulation cycle (i.e., Figure 3.1), YADE clearly distinguishes between two class families, providing code flexibility. Each simulation time step involves essentially two classes, the so-called data components and function components. The first class, as its name implies, contains only data, while the second class comprises functions that operate on the data.

The data component class includes the discrete bodies, their interactions, and the generalized forces. Discrete bodies are elements that interact with each other and constitute the object of simulations. Shape, material, state, and bound properties characterize each element. The shape defines the geometry of elements such as spheres, facets, walls, etc. The material sets up characteristics related to mechanical behavior such as Young's modulus and density. This property is independent of the geometry and size of the elements and can be shared by different discrete bodies. The state defines the spatial configuration and internal variables of elements, for example, linear and angular velocities and accelerations. State attributes are recalculated with each simulation cycle. Collision detection tasks between elements use the bound component. It describes the part of the space occupied by the element and may be larger than the element itself (e.g., approximate collision detection).

Interactions always involve a pair of discrete bodies. Depending on the proximity between elements, the collision detector creates the interactions. It is also possible to create interactions manually (i.e., explicitly) independently of the distance between elements. Each YADE interaction has two other components, the so-called IGeom and IPhys. The first component is related to the geometric setting of the interaction, in other words, it interprets the geometry of interacting elements (e.g., the contact between two spheres, or between a sphere and a cylinder). The IPhys component specifies the mechanical properties of interactions, such as contact stiffness, usually determined as a function of the mean between Young's modulus values of the interacting elements.

Generalized forces container includes forces and torques acting on each element as well as the storage of both forced displacement and rotation. This data is temporarily stored in the container during the corresponding time step and reset each new simulation cycle. Whenever an element has more than one interaction for a time step in progress, balancing the forces of those interactions provides the resultant force in the element.

Concerning the class of functional components, processes in charge of actions in the simulation characterize this class. This class includes essential functions (i.e., engines), process dispatchers or simply dispatchers, and functors (i.e., mappable data or functions). There are two fundamental types of engines characterized as global and partial. Global engines act throughout the simulation, for example, the mechanism that controls the gravitational force in the simulation cycles from the mass values of each element. Partial engines operate only on preselected discrete bodies, for example, a mechanism that applies a constant force to a specific group of elements.

Finally, dispatchers do not perform computations by themselves. As its name suggests, dispatchers trigger other functions represented by function objects or functors. Each functor is capable of handling a particular type of object and can be dispatched only if the dispatcher also considers the corresponding object. For example, creating the IGeom component of the interaction between a sphere and a facet uses dispatchers and functors specific to this collision configuration.

As an illustration, in collision detection, it is necessary to calculate the bound of all elements of the simulation. For an approximate collision detection task where the bounds are of the axis-aligned bounding box type (i.e., AABB), the algorithms for calculating these bounds differ depending on the shape of the elements. Thus, simulation requires specific functors for each type of object in the model and dispatchers for all possible combinations for interaction between these objects. The following sections describe how simulations occur in YADE and present preliminary models that provided a consistent understanding of code operation.

DEM simulations in YADE

Interaction procedure in simulations

The interaction between elements is what governs any numerical model involving DEM. In YADE, for each simulation cycle, the interaction between two discrete elements consists of the steps:

For new interactions

• to detect collision between elements;

• to create new interactions and to determine their properties (i.e., values derived from the parameters of interacting elements);

For existing interactions

• to evaluate displacements and strains;

• to calculate forces as a function of displacements;

• to apply forces on the elements and to determine the motion.

The following sections contain the description of each of these steps in YADE. This explanation applies the same symbols used by Šmilauer et al. (2015a). Also, in the last part of this topic, two additional sections deal respectively with numerical damping and calculating critical time steps.

Collision detection between elements

Considering a generic pair of discrete bodies symbolized by i and j, with their geometric shapes represented by the sets of points P i and P j , respectively, it is possible to detect the collision between these bodies by analyzing overlapping points. Two stages involving potential and exact collisions describe this analysis:

1. First, fast collision detection is performed by using less refined shape point sets, P i and P j . A rough characterization of the space occupied by the elements and the statements of Eqs. 3.1 and 3.2 define these two sets of points.

∀x ∈ R 3 : x ∈ P i ⇒ x ∈ P i (3.1) ∀x ∈ R 3 : x ∈ P j ⇒ x ∈ P j (3.2)
P i and P j point sets represent bounds in YADE (i.e., element bounding volume). They follow the implications of Eqs. 3.3 and 3.4.

(P i ∩ P j ) = 0 ⇒ P i ∩ P j = 0 (3.3)

P i ∩ P j = 0 ⇒ (P i ∩ P j ) = 0 (3.4)
2. The previous process involving fast detection refuses collisions that cannot occur (i.e., Eq. 3.4) and an exact (i.e., computationally more expensive) collision detection algorithm analyzes the remaining potential collisions by filtering out the fictitious ones by the condition of Eq. 3.5.

P i ∩ P j = 0 ∧ (P i ∩ P j ) = 0 (3.5)
The algorithms that operate the sets of points are capable of handling combinations between different element geometries. The collision detection task is divided into two steps to computationally speed up this task. YADE applies the sweep-and-prune -SAP (ERICSON, 2004) algorithm, which uses the axis-aligned bounding box, to detect collision between elements. Šmilauer et al. (2015a) also describes this algorithm. Figure 3.2 illustrates possible scenarios of the described procedure. Considering a two-dimensional example with two disks represented by the P i and P j point sets and their less refined geometric spaces represented by the circumscribed squares, P i and P j , the collision will only happen if there are overlapping areas between sets of points that define the geometric shapes of the elements along the axes (i.e., potential collisions for the intersection of P i and P j and exact collisions for the intersection of P i and P j ). In the scenario of Figure 3.2b, there was only overlap on the x-axis, so no collision occurs. In Figure 3.2c, the fast detection stage indicates a virtual collision because of the overlapping of P i and P j ; however, the points that accurately represent the geometry of the elements do not overlap (i.e., P i ∩ P j = 0). It is important to emphasize that collision detection is approximate in the first stage, resulting from this process only potential interactions. The exact calculation of the collision depends on the individual geometry of the elements and is performed separately using specific algorithms for each type of combination. Šmilauer et al. (2015c) details the scheme of these algorithms.

Interaction establishment

A feature of DEM is to relate the mechanical behavior of a medium to the properties of its component particles. The establishment of interactions also depends on nongeometrical properties such as stiffnesses and friction acting on the interacting elements. A usual YADE's interaction has normal stiffness, shear (i.e., tangential) stiffness, and inter-particle friction angle. Figure 3.3 provides a representation of the contact between two interacting spheres without a damping effect (i.e., the model can also incorporate this property). In this case, the parameters K N , K T , and µ are the normal stiffness, the shear stiffness, and the coefficient of friction, respectively. FIGURE 3.3 -Simplified representation of the contact between two spheres. Adapted from [START_REF] Coetzee | Calibration of granular material parameters for DEM modelling and numerical verification by blade-granular material interaction[END_REF].

By the spring contact model (i.e., Figure 3.3), it is possible to relate the Young's modulus of the elements in interaction with the normal stiffness, K N . Thus, the shear stiffness of the contact, K T , is a given fraction of the computed K N and the Poisson's ratio of the arrangement is the K N /K T ratio (i.e., by dimensional analysis).

In YADE, a series arrangement of two springs equal in length to the interacting sphere radii characterizes the normal stiffness that appears in the interaction (i.e., Figure 3.4). It is important to note that this representation depends on the algorithm used, thus depending on the type of interaction, YADE can use other algorithms to calculate stiffness values. For Figure 3.4, the calculation hypotheses define the spheres as rigid bodies. Thus, the elements overlapping verified by reducing the length l (i.e., Figure 3.4), for example during a collision, indicates strain values proportionally distributed in each of these elements. For the case involving particles of a medium, it is common on this scale to refer to displacement (i.e., ∆l 1 e ∆l 2 for the given example) rather than strain. This choice depends on the conceptual model.

In Figure 3.4, the contact force, F N , is equal to the resultant forces on each element (i.e. F N 1 e F N 2 ) and is calculated by the product between the contact stiffness and the sum of particle strains. Normal force at contact is a function of stiffness in the same direction as defined in Eq. 3.6.

F N = F N 1 = F N 2 = K N (∆L 1 + ∆L 2 ) (3.6)
Each sphere has stiffness associated with its mechanical properties. Thus, it is possible to calculate the normal stiffness of Eq. 3.6 as a function of Young's moduli and lengths proportional to the radii of the spheres (i.e., Eq. 3.7).

K N = E 1 l1 • E 2 l2 E 1 l1 + E 2 l2 (3.7)
where E 1 , E 2 = Young's modulus of spheres; l1 , l2 = length proportional to radius and used in calculating spheres stiffness, that is,

K i = E i li with i ∈ {1, 2}.
In YADE there are different classes for calculating interaction properties. One of the most used is the class that characterizes frictional materials (i.e., YADE C ++ class Ip2 -FrictMat FrictMat FrictPhys). This class considers the lengths l1 and l2 equal to the respective sphere diameters, that is, li = 2 • r i with i ∈ {1, 2}. Some classes define an equivalent cross-section that appears in li , for example, YADE C ++ class Ip2 CpmMat -CpmMat CpmPhys for models involving concrete ( ŠMILAUER, 2010). In this case, li is equal to the quotient between the equivalent cross-section and the length given by the contact point to the center of the element. These examples illustrate the possibility of applying different considerations depending on the type of formulation used in the model.

Concerning non-elastic interaction parameters, the respective average values of the material properties of the elements generally define them (e.g., Ip2 CpmMat CpmMat -CpmPhys class averages most quantities). An exception to this procedure is in class Ip2 FrictMat FrictMat FrictPhys itself, which defines the internal friction angle as the minimum value found between interacting elements. In this case, this is preferable as it makes it possible to model elements without friction.

Finally, it is important to highlight that interactions between elements also depend on non-geometrical properties defined by the characteristics of the materials involved. Therefore, specific functions are responsible for assigning these properties for each new interaction.

Displacement and strain evaluation

The evaluation of displacements and strains considers the degrees of freedom existing in interactions between elements. Considering a three-dimensional space, each body or element can have up to six degrees of freedom (i.e., 6 DOF). That is, in a 3D space represented by the x, y, and z axes, displacement and rotation can occur in all three directions. In the case of interactions between elements, they can also move in space with six degrees of freedom, as illustrated in Figure 3.5 (i.e., basic interaction configuration defined by two spheres in contact). Normal strains occur when there is, along the axis of interaction (i.e., axis passing through the center of the elements), the relative displacement between contacting elements. For this condition, there is 1 DOF (i.e., Figure 3.5a). In the case of shearing caused by orthogonal relative displacements to the interaction axis, including those generated by angular velocity, there can be 2 DOF, as shown in Figure 3.5b. Lastly, twisting has 1 DOF and bending 2 DOF. Figures 3.5c and 3.5d show displacements caused by the angular velocity difference parallel to the interaction axis (i.e., twisting) and perpendicular to the interaction axis (i.e, bending), respectively. In YADE, most of the constitutive laws used in element-to-element contacts do not include twisting and bending components. Therefore, this section focuses on the general description of normal and shear component evaluation.

The typical case in which two spheres are initially in contact with some overlap between them (i.e., Figure 3.6) helps to exemplify the calculation of the normal component. From the initial position of the respective sphere centers and the interaction reference point, defined in the middle of the elements or the center of the overlapping area, it is possible to describe the following state variables that update as the spheres move:

n o = C o 2 -C o 1 = C o 2 -C o 1 |C o 2 -C o 1 | (3.8) C o = C o 1 + d 1 - d 0 -|C o 2 -C o 1 | 2 n o (3.9)
where n o = the unit vector in the direction of the interaction axis at the current calculation time step; It is important to note that the geometric radii r 1 and r 2 differ from interaction radii d 1 and d 2 when elements overlap. From Figure 3.6 and considering C o 1 , C o 2 , and d 0 , the normal displacement u N and the normal strain ε N can be defined, respectively, by the following equations:

C o 1 , C o 2 =
u o N = |C o 2 -C o 1 | -d 0 (3.10) ε o N = u o N d 0 = |C o 2 -C o 1 | d 0 -1 (3.11)
where u o N = normal displacement at the current calculation time step;

ε o N = normal strain at the current calculation time step.

In order that the repulsive force does not remain finite as the centers of both spheres approach, preventing them from penetrating each other, it is common to use logarithmic strain. Eq. 3.12 presents this adjustment for normal strain so that in the case of massive compression, the strain tends to -∞ instead of -1.

ε o N =        log |C o 2 -C o 1 | d 0 , if |C o 2 -C o 1 | < d 0 |C o 2 -C o 1 | d 0 -1 , otherwise (3.12)
Eq. 3.12 has the drawback of effectively increasing the stiffness of the contact, causing this value to tend to infinity as the elements overlap. Thus, it is necessary to control the time step of the calculation cycle so that the simulation does not become unstable. It is possible to control it dynamically through the YADE's engine GlobalStiffnessTimeStepper ( ŠMILAUER et al., 2015b), which acts as a time step stiffness-based manager.

In the case of the shear component, characterized by the displacement perpendicular to the interaction axis, its calculation must consider the mutual spatial motion of the elements, so that the value is consistent in cases where the elements move with constant configuration between them. YADE C ++ class yade.wrapper.ScGeom uses the classical incremental algorithm to calculate the shear component. Several works describe this algorithm, which is commonly used in DEM codes, for example, [START_REF] Luding | Introduction to discrete element methods: basic of contact force models and how to perform the micro-macro transition to continuum theory[END_REF] and Alonso-Marroquín et al. (2004).

The incremental algorithm updates the shear displacement u T by applying an increment whose determination can be divided into two parts. The first part depends on the interaction movement, characterized by the state changes of the interaction reference point and the unit vector in the direction of the interaction axis (i.e., rotation). The second part takes into account the relative displacement perpendicular to the interaction axis.

The change in the state of interaction is due to the movement of the elements (i.e., displacement of element centers). Eqs. 3.8 and 3.9 provide the unit vector in the direction of the interaction axis and the current state of the interaction reference point. Figure 3.7 illustrates the update of the shear component from the previous state u - T to the current state (i.e., u o T ). The u - T component is perpendicular to vector n -and, after updating the interaction state, u o T must also be perpendicular to n o . Besides, the increment caused by the relative rotation must be added to obtain the current component. Eq. 3.13 provides the perpendicularity condition and Eq. 3.14 the increase due to the difference in rotation between the elements.

∆u Concerning the relative displacement of interacting elements in the direction perpendicular to the interaction axis (i.e., shearing), the calculation depends on the linear velocity component in the same direction. Thus, the relative linear velocity and its perpendicular component to n o can be calculated by Eq. 3.15 and Eq. 3.16, respectively.

T 1 = -u - T × n -× n o (3.13) ∆u T 2 = -u - T × 1 2 ω 1 + ω 2 • ∆t • n o • n o (3.
v 12 = v 2 + ω 2 × (-d 2 • n o ) -v 1 + ω 1 × (d 1 • n o ) (3.15) v ⊥ 12 = v 12 -(n o • v 12 ) • n o (3.16)
where v 12 = relative linear velocity between the two elements;

v 1 , v 2 = respective linear velocities of elements with index 1 and 2;

v ⊥ 12 = relative linear velocity component perpendicular to the interaction axis.

Therefore, the increase in shear displacement corresponding to the relative displacement between the two elements in the direction of the interaction axis is calculated by Eq. 3.17. Eq. 3.18 gives the shear displacement at the current calculation time step. Finally, it is possible to calculate the shear strain by the ratio between the result of Eq. 3.18 and d 0 (i.e., Eq. 3.19).

∆u

T 3 = v ⊥ 12 • ∆t (3.17) u o T = u - T + ∆u T 1 + ∆u T 2 + ∆u T 3 (3.18) ε o T = u o T d 0 (3.19)
where u o T = shear displacement at the current calculation time step;

ε o T = shear strain at the current calculation time step.

Force calculation

The calculation of forces depends on the nature of the modeled material. Thus the attainment of these values results from a constitutive law. Constitutive laws can be expressed using formulations involving stresses and strains as well as forces and displacements. In YADE, all constitutive laws derive from the base class LawFunctor ( ŠMILAUER et al., 2015b). One of the most commonly used in DEM formulations is the law originally proposed by [START_REF] Cundall | A discrete numerical model for granular assemblies[END_REF] for non-cohesive elastic behavior with Coulombtype friction (i.e., no lubricating fluid), which is implemented in YADE C ++ class Law2 -ScGeom FrictPhys CundallStrack.

In the calculation of forces, when there is a new interaction between elements, the materials fixed for these elements define the physical properties of the interaction. Considering, for example, the simple case of frictional materials in contact, the YADE C ++ class Ip2 FrictMat FrictMat FrictPhys creates a new FrictPhys instance (i.e., interaction between frictional materials) that defines normal stiffness K N , shear stiffness K T , and friction angle φ. From these interaction parameters, the previously defined normal and shear displacements, and the law written in the class Law2 ScGeom FrictPhys CundallStrack, it is possible to calculate normal and shear forces by Eqs. 3.20 and 3.21, respectively. In Eq. 3.21, a threshold value is considered using a Coulomb-type friction law.

F N = K N • u o N • n o (3.20) F T =      F t T • |F N |•tan(φ) F t T , if |F T | > |F N | • tan (φ) F t T , otherwise (3.21)
where

F N = normal contact force;
F T = contact shear force;

F t T = trial contact shear force (i.e., F t T = K T • u o T ).
When u o N is greater than zero, the term |C o 2 -C o 1 | is greater than the sums of interaction radii (i.e., Eq. 3.10). In this case, the algorithm removes the contact, and no force is generated. Assuming there is contact, the sum of the forces calculated in Eqs. 3.20 and 3.21 gives the total force F applied to each of the two elements (i.e., Eqs. 3.22 and 3.23). This total force acts at the interaction reference point, and therefore can generate torque on both elements as defined by Eqs. 3.24 and 3.25. Thus, each element of the model can concentrate forces and torques (i.e., Eqs. 3.22 to 3.25) acting on them in each simulation cycle.

F + 1 = F (3.22) F + 2 = -F (3.23) T + 1 = d 1 • (-n o ) • F (3.24) T + 2 = d 2 • n o • F (3.25)
where F + 1 , F + 2 = respective total forces on elements with index 1 and 2 to be applied in the next calculation time step; F = total force (i.e., F = F N + F T );

T + 1 , T + 2 = respective torques on elements with index 1 and 2 to be applied in the next calculation time step.

Motion calculation

As described above, during interactions, generalized forces (i.e., forces and torques) are generated, which are applied individually to the contacting elements. From these forces, it is possible to integrate equations of motion for each element or particle separately. For time-domain numerical integration, YADE applies the leapfrog scheme [START_REF] Fincham | Leapfrog rotational algorithms[END_REF], which is widely used in numerical methods. Leapfrog integration uses the idea that derivation at midpoint provides a better approximation than that obtained at extreme points. This feature makes it possible to know position and orientation derivatives at full-step points and velocity derivatives at mid-step points in a sequence of calculation steps (i.e., the values leap over each other). This integration algorithm is implemented in the YADE C ++ class yade.wrapper. NewtonIntegrator ( ŠMILAUER et al., 2015b).

The new position of the element relative to the current state (i.e., next time step) depends on its calculated acceleration to the current time step, as described in Section 2.3 of Chapter 2. From this and Newton's laws of motion, it is possible to write: By applying the Finite Difference Method -FDM, it is possible to obtain an approximation for acceleration (i.e., the second derivative of position with respect to time), allowing truncation in the second-order term. Eq. 3.27 provides this approximation, in which the assumed time interval is small enough to make the truncation error negligible.

üo = u + + u --2u o ∆t 2 (3.27)
where u + = element position at the next calculation time step (i.e., t + ∆t);

u -= element position at the previous calculation time step (i.e., t -∆t);

t = current time; ∆t = time step.
Isolating the term relative to the position of the element for the next step gives:

u + = 2 • u o -u -+ üo • ∆t 2 = u o + ∆t • u o -u - ∆t + üo • ∆t (3.28)
where u o = element position at the current calculation time step.

The average velocity of the element calculated at time t -∆t 2 (i.e., previous step) can be approximated by:

u ∼ = u o -u - ∆t (3.29)
where u = average element velocity defined at the mid-step point of the previous step (i.e., t -∆t 2 ).

The term in parentheses in Eq. 3.28 defines the average velocity of the element for the current step (i.e., the calculation involves the mid-step point of the next step). Using the relationship given in Eq. 3.29 in the term in parentheses (i.e., Eq. 3.28) allows writing Eq. 3.30. u⊕ = u + üo • ∆t (3.30) where u⊕ = average element velocity defined at the mid-step point of the next step (i.e., t + ∆t 2 ).

Therefore, Eq. 3.28 results:

u + = u o + u⊕ • ∆t (3.31)
To summarize, the motion algorithm first calculates the average element velocity required to calculate the next step and then calculates the new element position for the next step by Eq. 3.31. Thereby, this algorithm allows knowing positions at (i • ∆t) and velocities at (i • ∆t + ∆t/2) with i ∈ {1, 2, ..., n cs } and n cs equal to the number of calculation cycles.

It is possible to update the orientation of the elements analogously to the position by calculating the angular acceleration which, in turn, depends on the calculated torque (i.e., Section 3.2.5). Considering spherical elements, which have a diagonal inertia tensor in any orientation, it is possible to calculate the angular acceleration of the element by Eq. 3.32 (i.e., Newton's law for rotation).

ωo i = T i I 11 , with i ∈ {1, 2} (3.32)
where ωo i = element angular acceleration at the current calculation time step;

T i = current torque on element;

I 11 = moment of inertia (i.e., diagonal elements of the particle inertia tensor which for spherical particles satisfies I 11 = I 22 = I 33 ).

The same scheme described above calculates the angular velocity of the element by:

ω ⊕ i = ω i + ωo i • ∆t (3.33)
where ω ⊕ i = element angular velocity defined at the mid-step point of the next step (i.e., t + ∆t 2 ); ω i = element angular velocity defined at the mid-step point of the previous step (i.e., t -∆t 2 ).

The number system first presented by [START_REF] Hamilton | On quaternions; or on a new system of imaginaries in algebra[END_REF] and described by several authors [e.g., [START_REF] Goldman | Rethinking quaternions: theory and computation[END_REF]], so-called Quaternions, makes it possible to represent the rotation vector given by ω ⊕ i • ∆t through the quaternion ∆q. Thus, the terms of ∆q regarding rotation and direction of rotation are given respectively by:

∆q θ = ω ⊕ i (3.34) ∆q u = ω ⊕ i (3.35)
where ∆q θ = quaternion rotational term; ∆q u = quaternion directional term;

ω ⊕ i = normalized vector ω ⊕ i .
Finally, the rotation composition of Eq. 3.36 calculates the updated orientation of the element.

q + = q o • ∆q (3.36) where q + = element orientation at the next calculation time step (i.e., t + ∆t);

q o = element orientation at the current calculation time step; ∆q = quaternion of the rotational vector.

In YADE, motion calculation is also possible for elements modeled by rigid aggregates of particles, so-called clumps, and non-spherical elements. These cases include the same idea presented above with some adaptations and an extension of the leapfrog integration algorithm for calculating element orientation ( ŠMILAUER, 2010; ŠMILAUER et al., 2015a).

Numerical damping

Assuming elastic behavior, the energy transformed during the interaction between elements may not dissipate. Therefore, numerical damping may be required to achieve a static or steady-state solution in some cases. Aiming at dissipating the kinetic energy of discrete elements, particularly in quasi-static case simulations, it is possible to adjust the equations of motion to allow damping in the model.

The main idea of artificial numerical damping is to decrease the forces that increase particle velocities by considering their acceleration and velocity directions ( ŠMILAUER et al., 2015a). This process makes the damping scheme non-physical, as it is not invariant to rotation in a coordinate system; on the other hand, it is not difficult to calculate. Thus, it is possible to introduce a force decrement ∆F in Eq. 3.26, defined from a dimensionless damping coefficient as follows:

∆F ξ F ξ = -λ d • sgn F ξ • u ξ , with ξ ∈ {x, y, z} (3.37)
where ∆F ξ = force decrement components;

F ξ = force components acting on the element;

λ d = dimensionless damping coefficient;
u ξ = element velocity components defined at the mid-step point of the previous step (i.e., t -∆t 2 ).

The coefficient λ d must be between zero (i.e., no damping) and one (i.e., no acceleration). YADE uses a version adapted for Eq. 3.37, which replaces the previous mid-step velocity u ξ with its current time step estimate uo ξ , as shown in Eq. 3.38. This avoids lockin effect that can occur if the velocity changes its sign due to force application at each step (i.e., when the element oscillates around the position of equilibrium over a period 2∆t). et al., 2015b). This formulation has some conveniences, such as acting on forces (i.e., accelerations) without constraining uniform motion, being independent of particle eigenfrequencies (i.e., all particles are equally damped), and depending only on the dimensionless parameter λ d that does not require calculation.

∆F ξ F ξ = -λ d • sgn F ξ • uo ξ = -λ d • sgn F ξ • u ξ + üo ξ • ∆t 2 , with ξ ∈ {x, y, z} (
YADE also features viscous damping [START_REF] Ng | Input parameters of discrete element methods[END_REF] as used in the studies by [START_REF] Albaba | Modeling the impact of granular flow against an obstacle[END_REF] and [START_REF] Haustein | Discrete element modeling of deformable particles in yade[END_REF]. This type of damping acts on the velocity of the element. In YADE, different constitutive laws use this type of model; an example is the one implemented in the class Law2 ScGeom FrictViscoPhys CundallStrackVisco, which corresponds to the class Law2 ScGeom FrictPhys CundallStrack by incorporating the viscous damping model in the normal direction ( ŠMILAUER et al., 2015a).

Critical time step

Setting the value of ∆t and damping the equations of motion as described in the previous section are two key issues of DEM. The time step has to be less than a critical value for the second-order time-centered finite difference scheme of [START_REF] Cundall | A discrete numerical model for granular assemblies[END_REF] to produce a stable solution. The critical time step is related to the maximum eigenfrequency (i.e., minimum eigenperiod) of the total system. Its estimate follows the procedure of [START_REF] Hart | Formulation of a three-dimensional distinct element model -Part II. Mechanical calculations for motion and interaction of a system composed of many polyhedral blocks[END_REF]. For stability to the explicit integration scheme, it is possible to define an upper limit for ∆t by Eq. 3.39. ∆t cr = 2 ω max (3.39) where ∆t cr = critical time step;

ω max = the maximum eigenfrequency of the system.

Considering the spring contact model described above (i.e. mass-spring system), it is possible to rewrite Eq. 3.26 as: The solution of harmonic oscillation for the second-order linear differential equation (i.e., Eq. 3.40) is: By substituting the second time derivative of Eq. 3.41 into the differential equation obtained directly from Eq. 3.40 results in an expression for angular frequency as a function of stiffness and mass (i.e., Eq. 3.44).

m • ü = -K • u (3.
u(t) = A • cos (ω • t + ϕ) ( 3
d 2 u dt = - K m • u (3.42) -ω 2 • A • cos (ω • t + ϕ) = - K m • A • cos (ω • t + ϕ) (3.43) ω = K m (3.44)
Elements can move at the same time; therefore, the highest system frequency occurs when there is opposite motion between two elements. For elements with the same velocities and displacements equal to ∆u, there is an apparent stiffness K a defined as follows:

∆F = -K • [∆u -(-∆u)] = -2 • K • ∆u (3.45) K a = 2 • K (3.46)
where ∆F = force increment; ∆u = element displacement increment;

K a = apparent stiffness of the system.

The maximum system eigenfrequency ω max is the maximum value of all frequencies calculated in the system, and its calculation corresponds to replacing K by apparent stiffness values in Eq. 3.44. Given this and considering Eq. 3.39, the global critical time step is the minimum value among those calculated for each element by applying Eq. 3.47 separately for each degree of freedom, assuming them not coupled [START_REF] Hart | Formulation of a three-dimensional distinct element model -Part II. Mechanical calculations for motion and interaction of a system composed of many polyhedral blocks[END_REF]CHAREYRE;[START_REF] Chareyre | Dynamic spar elements and discrete element methods in two dimensions for the modeling of soil-inclusion problems[END_REF].

∆t cr = 2 ω max = 2 • m K a = √ 2 • m K (3.47)
It is important to note that K represents the stiffness matrix determined by the different stiffnesses of each interaction of a given element. For greater computation efficiency, YADE's engine GlobalStiffnessTimeStepper simplifies this approach by considering only translational terms of the stiffness matrix, which is usually the limiting factor, and there-fore estimates eigenvalues exclusively from diagonal terms ( ŠMILAUER et al., 2015a). As mentioned by [START_REF] Effeindzourou | Modelling of deformable structures in the general framework of the discrete element method[END_REF], in typical mechanical systems, this procedure provides very accurate estimates of the stability condition observed in benchmark tests.

In case there are no particle contacts, the maximum eigenfrequency of the system is zero (i.e., ω max = 0), and the critical time step tends to infinity. Therefore, to avoid numerical errors under these conditions, the engine GlobalStiffnessTimeStepper estimates in advance stiffness based on future interactions ( ŠMILAUER et al., 2015a).

Preliminary models and simulations

Development

The following sections present some of the modeling possibilities using YADE that may be useful in future work related to the topic of this thesis. These sections cover the preliminary simulations that, based on the DEM concepts studied (CUNDALL; STRACK, 1979; MUNJIZA, 2004; O'SULLIVAN, 2011), served to gain familiarity with YADE. All the simulations involved important fundamentals that allowed to manipulate YADE consistently in the works of this thesis. The models, here characterized as examples, require general notions of Python programming language (LUTZ, 2009) and YADE framework ( ŠMILAUER et al., 2015c) for the definition of simulation scripts.

The characteristics and parameters of the materials are purely hypothetical. Their determination is due solely to the purpose of the preliminary models and simulations. Therefore, this part does not cover the evaluation of these materials, focusing only on the ideas of numerical models. It is important to emphasize that the results presented in the following sections, including those with calculated and graphed values, have no other purpose than to understand the functioning of the codes and how to access some simulation data. In this regard, focusing on qualitative analysis, the examples omitted numerically calculated result units. Thus, the examples performed were:

• Contact between spheres;

• Spheres gravity deposition;

• Oedometer test;

• Simple shear;

• Triaxial test;

• Grid pullout test;

• Sample porosity control;

• Handling interactions.

Contact between spheres

The first simulation reproduces a vertically falling sphere that collides with another sphere of the same diameter fixed at a lower position in space and perfectly aligned (i.e., support sphere), as shown in Figure 3.8. Figures 3.8a and 3.8b show the initial and final state of the simulation, respectively. This simulation provided important aspects, such as fixing the position of an element in space and applying gravitational force to the system. The simulation also assumed energy dissipation obtained considering the numerical damping described earlier. Thus, a sequence of collisions between the two spheres occurred until the energy balance of the system. The yade.plot module contains utility functions for plotting values directly in YADE (i.e., Figure 3.8c). These functions have plot formatting limitations, although they can be improved. Figure 3.8c shows the evolution of the vertical displacement of the falling sphere (i.e., the upper sphere) during the simulation. The rebound height that decreases during the simulation is a function of the damping coefficient value.

Spheres gravity deposition

Gravity deposition comprised spheres of varying diameter arranged within a box, as shown in Figure 3.9. This simulation made it possible to explore dispatchers and functors responsible for determining the collision geometry, its mechanical properties, and the contact law for force application. It consisted of randomly arranging spheres of different diameters inside a box and then imposing gravitational force on the system. A calculated value of the unbalanced force determined the end of the simulation, that is, a target value for which the simulation cycles ceased. The ratio between the mean or maximum force acting on bodies (i.e., particles) and the average force acting on interactions give the unbalanced force. Depending on the boolean parameter useMaxForce of the unbalancedForce() function of YADE's module yade.utils, the numerator of this ratio will be the average value or maximum value. This ratio provides a measure of the system's static equilibrium condition. As the simulation stabilizes, the unbalanced force tends to zero. Sufficiently small value can be approximately 0.01 or even less, depending on how much balance is required. This model also included the calculation of work associated with gravity force, kinetic energy, elastic potential, and dissipation energies caused by non-viscous damping and plastic dissipation. The chart in Figure 3.9c shows these values as gravWork, kinetic, elastPotential, nonviscDamp, and plastDissip, respectively. This chart also includes the change in unbalanced force (i.e., unbalanced term in Figure 3.9c) relative to the interaction number throughout the simulation.

In Figure 3.9c, the unbalanced force fluctuates due to particles bouncing during interactions. Gravitational work (i.e., work against gravity) is maximum at the beginning of the simulation and tends to a negative value at the end of gravity deposition. The negative value indicates that the spheres have lost gravitational potential energy relative to the final arrangement in the box. As a result, the kinetic energy increases during particle drop and tends to zero on the particle resting condition at the end of the simulation. The elastic potential (i.e., elastPotential) is the elastic energy stored in contacts during interactions. The value is low during simulation because the forces that maintain active contacts are not high enough to cause significant overlap between elements. Non-viscous damping (i.e., nonviscDamp) is the damping dissipation. Plastic dissipation (i.e., plastDissip) is the total energy dissipated when frictional contact between elements slips.

Oedometer test

The model consisted of simulating the gravity deposition described above and, from an equilibrium condition, applying compression and decompression stages by imposing vertical displacement on a horizontal plate at the top of the sample. This simulation emphasized the use of real-time simulation and the number of iterations to control the periodicity of model functions, the constant velocity imposition on a given element, the contact force calculation for use in a displacement checker function (e.g., function to control the compression plate), and also the YADE batch mode. As its name implies, the batch mode allows sequential processing of a simulation with different parameters, facilitating parametric studies.

In this model, a threshold value for the contact force on the loading plate defined the maximum allowable compression and the beginning of the unloading stage (i.e., decompression stage). Figure 3.10 shows the initial condition and the compression and decompression stages of the oedometer test simulation, as well as the evolution of the calculated forces. Increasing the vertical force during loading implies compacting the numerical sample. Figure 3.10d shows that decreasing the vertical force during unloading results in a nonreversible vertical strain (i.e., plate displacement). Besides, it is possible to verify a hysteresis effect on the loading and unloading curves. The intensity of unbalanced forces is a function of the loading and unloading velocity applied to the sample.

An approximately zero value for the contact force on the loading plate defined the end of the simulation. When the loading plate practically loses contact with the sample, the sample dissipates the residual energy stored in the contacts, which results in the increased unbalanced forces observed at the end of the simulation (i.e., Figure 3.10d). In this case, the end of the simulation was before the total residual energy dissipation of the sample. Thus, the unbalanced forces are not equal to or close to zero at the end of the oedometer test.

Simple shear

The simple shear model used the YADE's periodic space feature to initially impose an isotropic normal strain and, subsequently, a constant distortion rate in the periodic cell created in the model. The repetition of a parallelepiped-shaped cell characterizes a periodic space. The concept is to satisfy periodicity conditions to avoid boundary effects in three-dimensional DEM simulations. YADE C ++ class yade.wrapper.Cell ( ŠMILAUER et al., 2015b) comprises the implementation of this functionality.

There are sphere packing generators in YADE that are useful for modeling volumes or samples. An example is the regularHexa (i.e., YADE's module yade.pack), which has a higher density from a hexagonal close packing arrangement (CONWAY; SLOANE, 1999). A regular hexagonal packing of spheres generated from a periodic random cloud of spherical elements characterized the sample in this model.

After normal isotropic deformation, a constant distortion rate applied to the periodic cell imposed the simple shear condition on the sample. Figure 3.11 shows the sheared numerical sample with the colored spheres based on their amount of rotation. Qualitatively, the green spheres had the highest rotations, and the dark blue spheres had the lowest rotations. Figure 3.11a shows the particle rotations for the final state of the sample (i.e., particle rotation at the end of the simple shear simulation). Figure 3.11b presents the shear stress curve as a function of shear strain. The forces that cause shear are a function of contact friction forces between the spheres as well as the particle entanglement mechanisms (i.e., particle imbrication). In Figure 3.11b, the contact friction influences the maximum shear stress that represents the imminence of the particle rolling mechanisms. During shear, the particles may roll over each other and modify the initial state of imbrication. Subsequent peaks and valleys of the stress-strain curve in Figure 3.11b illustrate this behavior. The peaks are due to the difficulty of rolling between particles, and the valleys represent the range in which particle imbrication changes.

Triaxial test

The triaxial test simulation provided an understanding of the YADE's engine PeriTri-axController ( ŠMILAUER et al., 2015b) that controls stresses and deformations in periodic samples. This simulation consisted of imposing two strain stages on a periodically defined sample (i.e., periodic boundary). The first stage involved isotropic compression of the sample to a predefined stress value. The second stage comprised an imposition of uniform vertical strain, keeping lateral stresses constant, until reaching 20 % of deformation concerning the initial height of the sample.

In this model, rigid aggregates (i.e., clumps) created by three spherical elements characterized each particle in the sample. Particle modeling from aggregate elements makes it possible to create complex geometric shapes. Figure 3.12 illustrates the triaxial test and some measured values. Figure 3.12b presents curves of deviator stress versus vertical strain and variation of sample volume versus vertical strain. Importantly, the maximum deviator stress value is a function of the density of the granular assembly, microscopic friction, and particle shape.

Among the many methods available, triaxial compression testing is reliable and widely used in engineering. This test is useful for measuring the mechanical properties of materials such as shear strength and stress-strain behavior. The interest of simulating a triaxial test is to calibrate the microscopic parameters to obtain the real macroscopic parameters of the material (e.g., friction angle and dilatancy).

Grid pullout test

This simulation comprised a pullout test of a grid created from flexible cylinders (i.e., Figure 3.13), based on the studies of [START_REF] Chareyre | Dynamic spar elements and discrete element methods in two dimensions for the modeling of soil-inclusion problems[END_REF], [START_REF] Bourrier | Discrete modeling of granular soils reinforcement by plant roots[END_REF], and [START_REF] Effeindzourou | Modelling of deformable structures in the general framework of the discrete element method[END_REF]. Nodes and connections between them characterize the grid of this model. In short, grid modeling followed two steps, that is, first the definition of the mesh of nodes and later the connections of these nodes by cylindrical elements. The part dedicated to developing the numerical model of this thesis describes and provides details of this type of element. The grid between two layers composed of spheres of the same diameter configured the model sample. A plate imposed a confining condition on the sample by vertical loading to a predetermined threshold value for the calculated force on the loading plate (i.e., a procedure similar to the Section 3.3.4 oedometer test). Imposing a constant velocity on one end of the grid initiated the pullout test. Vertical color bands applied to sphere groups made it easy to visualize particle movements throughout the simulation. The simulation contemplated the pullout force measurement until the total grid extraction from the sample. Figure 3.13 presents simulation steps and the evolution of the force at the end of the grid. Later, the numerical analysis part of this thesis, which includes geogrid force-displacement curves, details this type of simulation.

Sample porosity control

This simulation aimed to model a layer of rigid aggregates of different shapes and sizes and specific porosity. It is similar to the gravity deposition of Section 3.3.3, but with particles formed by rigid aggregates and with porosity control. Two measurements defined the range of possible porosities for the sample. The first one is the initial porosity calculated after stabilization of the gravity deposition stage, which defined the maximum possible porosity of the arrangement. The second one is the minimum porosity calculated after sample consolidation by the self-weight of the particles, imposing a zero friction angle between them. Thus, the porosity specified for the final state of the sample considered this value range.

Spheres and three different geometries of rigid aggregates characterized the model particles, which also varied in size. In YADE, it is possible to calculate the porosity by the utils.porosity() function of the yade.utils module. This function calculates the total volume by assuming an axis-aligned bounding box for the sample bounds and, for simulations involving a periodic boundary, by the periodic cell volume itself. From the total volume and particle volume, it is possible to calculate the volume of void-space, and hence the porosity. Thus, the final state of the simulated sample comprised the reduction of inter-particle friction angle, including in pre-existing contacts, waiting for the system forces to balance, and calculating porosity in each simulation cycle until the specified value.

Figure 3.15 shows the initial condition of the simulation and the sample states im-mediately after particle gravity deposition and for the specified porosity, as well as the evolution of energy potentials, energy dissipations, unbalanced force, and sample consolidation. In the chart (i.e., Figure 3.15d), it is possible to note that porosity n ranged from 0.76 (i.e., initial porosity) to 0.20 (i.e., specified porosity) throughout the simulation.

Handling interactions

This simulation focuses on how to model different types of interactions and make the model capable of handling the contact possibilities between elements of different classes. The created model is fictitious and only aimed at simulating contacts covering all possibilities of interaction between spheres, cylinders (CHAREYRE; The model comprised two grids similar to that modeled in Section 3.3.7 and three free-falling spheres colliding with the grids. The grids were placed spaced one above the other with the upper one fixed in space by the shorter sides and the lower one fully fixed (i.e., nodes and connections). Particle facets element sealed two openings in the central region of the lower grid and one of the upper grid. This arrangement made it possible to verify during the simulation sphere-sphere, sphere-cylinder, sphere-facet, cylindercylinder, cylinder-facet, and facet-facet contacts. Figure 1 illustrates this simulation.

Introduction

The preparation of numerical samples has particular importance in the final results of simulations. As in experimental studies, samples also need to be carefully prepared or assembled in numerical studies, observing their state and incorporated properties. Therefore, the procedure involves assembly aspects such as mechanical equilibrium and homogeneity. In the case of modeling geosynthetic reinforcements (e.g., geogrid), the mechanical behavior of these structures and how they interact with other elements also need to be evaluated to ensure their correct use in the numerical model.

In DEM analyses, especially in quasi-static conditions that hold a strong influence on the initial state of the specimen, the sample preparation procedure should be even more meticulous. As discussed by (RADJA Ï; DUBOIS, 2011), there are numerical recipes to built a well-controlled, homogeneous, and representative sample. These recipes involve the evaluation of numerical and mechanical parameters considering the micro and macro domain of the sample, as well as the type of boundary condition of the model (e.g., non-periodic and periodic boundary).

In the case of geogrid pullout simulations of this thesis, geogrid-reinforced granular assemblies characterize the numerical samples. Thus, it is important to obtain adequate sample conditions (e.g., homogeneity and equilibrium) before pulling out the geogrid in the simulations. The following sections detail the discrete numerical approach to geogrid and granular material and also describe the techniques used in preparing samples for the simulations. It is important to emphasize that the auxiliary tool for discrete modeling used here is the open-source framework YADE ( ŠMILAUER et al., 2015c).

Geogrid modeling in the DEM framework 4.2.1 DEM discretization

Geogrids are geosynthetic materials with the function of reinforcing or improving the mechanical behavior of granular layers. They are generally made of polymers such as polyester, polyethylene, polypropylene, polyvinyl alcohol and aramid, which provide a feature of high tensile strength at low elongation. Furthermore, the grid frame adds a very important mechanism in the reinforced layer which is the interlocking in the region of the geogrid apertures. The granular material of proper size can pass through the apertures and the two materials interlock together resulting in a composite behavior.

Geogrids vary in type depending on the manufacturing method and may even have different geometric characteristics (e.g., extruded, welded, and knitted or woven geogrids). Although a realistic numerical representation is required, a modeling that contemplates all the minutiae of a particular type of geogrid ends up restricting its application to very specific cases, in addition to resulting in a high cost of computational processing. Hence, the approach described in this section was developed to capture the main attributes of real generic geogrids such as stretching, bending, shearing, and twisting due to their interaction with granular media (i.e., interlocking, friction, and abutment).

Stretching usually depends on the tensile stiffness of the geogrid longitudinal yarns, whereas shearing depends essentially on the shear stiffness of the nodes or crossing points. The bending and twisting effects result from both the tensile stiffness of the longitudinal yarns and the stiffness of the crossing points. Pullout mechanisms can induce the abutment (i.e., a bearing or buttress action) between the transverse yarns of the geogrid and the granular particles blocked and interlocked within the apertures. In this case, it depends on the relative value between the thickness of the yarns and the particle sizes. Friction on the surface of the geogrid elements and interlocking mechanism allow transmitting the internal forces of the granular media skeleton to the reinforcement. Figure 4.1 illustrates these effects.

The discrete geogrid modeling was implemented in the YADE code using the cylinders and particle facets described respectively in the works of [START_REF] Bourrier | Discrete modeling of granular soils reinforcement by plant roots[END_REF] and [START_REF] Effeindzourou | Modelling of deformable structures in the general framework of the discrete element method[END_REF]. Cylinders are the product of the morphological dilation of a sphere and a segment. Similarly, the particle facet element, so-called PFacet in YADE, is geometrically created by the sum of a triangular facet and a sphere. Figure 4.2 shows these two element types. In both elements, the constitutive relations are defined by the information stored in the GridNodes, which are centers of the spheres that compose the boundaries of the element.

In YADE, cylinder consists of two GridNodes with a connection between them (i.e., The basis for the equations governing the external and internal interactions of cylinder and PFacet elements follows the most common law in three-dimensional discrete models, which is the law defined for interaction between two spheres. The sphere-sphere contact model can be used because of a virtual sphere idea, which is described in more detail later, as well as the presence of GridNodes. Figure 4.3 illustrates the final aspect of the numerical geogrid. Both longitudinal and transverse members on geogrid are modeled using cylinder and PFacet elements. In the enlarged detail (i.e., Figure 4.3) it is possible to notice the GridNodes, GridConnections, and facets that constitute the members of the geogrid. This discrete modeling approach for geogrids results in a flattened shape representative of real geogrids and has geometric attributes considered essential to the rolling mechanisms present in the interaction of this type of reinforcement with granular media.

Contact model

Contact models need to be defining to take into account the surface interaction between two elements or bodies, for example between two granular particles, between a granular particle and geogrid elements (i.e., GridNode, GridConnection and PFacet), and possibly between two geogrid elements involving the same folded structure or distinct structures. In a usual DEM model, contact between two spheres exists as soon as their volumes overlap. That is, the distance between the two centers is less than the sum of the radii of the spheres. The magnitude of the overlap can be related directly to the contact force by a force-displacement law. This is suitable for the case of interaction with other elements (i.e., external interactions). Figure 4.4 illustrates a feasible spring scheme for two spheres in contact. For element deformations (i.e., internal interactions), the magnitude of the internal forces is related to the relative displacements between GridNodes belonging to the element. A basic contact law that takes into account elastic and frictional forces is quite reasonable to describe interaction of the geogrid with other elements. The law used for this purpose implements the classical linear elastic-plastic model of [START_REF] Cundall | A discrete numerical model for granular assemblies[END_REF] and is written in the YADE C ++ class Law2 ScGeom FrictPhys CundallStrack ( ŠMILAUER et al., 2015b). This law computes contact forces between two elements based on relative displacements resulting from the integration of the relative velocity components with respect to time. Relative velocities are calculated using reference points in each element, usually the points on the line connecting the centers and crossing the boundaries of the elements in interaction. By combining these relative quantities with the previous law, the normal and shear forces increments can be expressed as:

dF n = k c n • un • dt, (4.1
) In the case of spheres, the contact can be assumed as two springs arranged in series with lengths equal to their respective radii (i.e., Figure 4.4). Thus, the normal stiffness in the contact is calculated by the harmonic mean of the individual rigidities of these springs. Since k c s is directly related to the k c n through Poisson's ratio, the equations for the normal and shear stiffness are:

dF s = k c s • us • dt, ( 4
k c n = 2 • E 1 • R 1 • E 2 • R 2 E 1 • R 1 + E 2 • R 2 , ( 4.3)
k c s = 2 • E 1 • R 1 • ν 1 • E 2 • R 2 • ν 2 E 1 • R 1 • ν 1 + E 2 • R 2 • ν 2 , ( 4.4)
where

E 1 , E 2 = Young's moduli of each sphere in contact; R 1 , R 2 = spheres radii; ν 1 , ν 2 = Poisson's ratios.
For each time step, normal and shear forces are calculated by the respective sum of all forces determined for previous time steps (i.e., prior normal and shear forces) plus the increment values given by Eqs. 4.1 and 4.2, respectively. In the calculation of the shear force F s , a threshold value is considered using a Coulomb-type friction law without cohesion in our case and inter-particle friction angle φ µ (i.e., Eq. 4.5). If an absolute value of F s is larger than the threshold value, the shear force modulus is set equal to F s lim .

F s lim = F n • tan (φ µ ) , (4.5)
where F s lim = the limit value of the shear force;

F n = normal force;

φ µ = inter-particle friction angle (i.e., the smaller value between the two interacting elements).

The force concentrated on the contact F c is calculated by the sum of the normal and shear forces described above. F c is applied to each interacting elements, which can generate moments in addition to forces in the course of the calculation step. These resultant moments and forces are then used to integrate motion equations for each element separately, obtaining position and orientation data. Eq. ( 4.6) gives the expression for the moment of force (i.e., torque) due to the contact force in each sphere.

τ = R i - δ 2 • n × F c , (4.6)
where τ = moment of force; R i = sphere radius; δ = the penetration depth of the contact (i.e., overlap measure); n = normal unit vector (i.e., parallel to the axis of interaction that passes through the centers of the elements);

F c = the concentrated force applied to the contact.

The feasibility of using the spheres contact law for geogrid interactions lies in the idea of defining a virtual sphere within the cylinder or PFacet. The diameter of this sphere is assumed equal to the thickness of the element. Virtual sphere is inserted in such a way that its center is in the projection of the contact point (i.e., contact between cylinder or PFacet and a spherical element) on the plane that passes through the GridNodes of the respective element in interaction. From this, all virtual sphere status update (i.e., translational and rotational velocities) is interpolated and the forces and moments are then distributed linearly on the GridNodes. Figure 4.5 illustrates the concept of virtual spheres in sphere-cylinder and cylinder-cylinder interactions. In the case of PFacet, coordinates of the virtual sphere center can normally be defined from the concept of barycentric coordinates. For the virtual sphere within the boundary of the element facet area, the coordinates of its center are calculated as weighting functions in relation to the reference points (i.e., the corresponding GridNodes) by Eqs. 4.7 to 4.9. YADE uses the method detailed by [START_REF] Ericson | Real-time collision detection[END_REF] for a more efficient implementation of these equations.

p 1 = Ψ P N 2 N 3 Ψ N 1 N 2 N 3 , (4.7)
p 2 = Ψ N 1 P N 3 Ψ N 1 N 2 N 3 , (4.8
)

p 3 = 1 -p 1 -p 2 , (4.9)
where p 1 , p 2 , p 3 = barycentric coordinates of the virtual sphere center;

N 1 , N 2 , N 3 = GridNodes; P = position of the center of the virtual sphere;

Ψ N 1 N 2 N 3 = area of the triangle formed by the GridNodes; Ψ P N 2 N 3 , Ψ N 1 P N 3 = areas of two sub-triangles formed by P and GridNodes.

Translational and rotational velocities of the PFacet's virtual sphere can be calculated from the respective values obtained at the nodes N 1 , N 2 , and N 3 as written in Eqs. 4.10 and 4.11. As regards the contact force calculated for the virtual sphere, it is distributed to the GridNodes by the application of Eq. 4.12. Figure 4.6 shows the virtual sphere in a sphere-PFacet interaction.

υ vs = 3 i=1 p i • υ N i , (4.10
)

ω vs = 3 i=1 p i • ω N i , (4.11)
F N i = p i • F vs , (4.12)
where υ vs = the translational velocity of the PFacet virtual sphere;

υ N i = the translational velocity of the node N i ;

ω vs = the rotational velocity of the PFacet virtual sphere;

ω N i = the rotational velocity of the node N i ;

F N i = distributed contact force at N i with i ∈ {1, 2, 3};

p i = the virtual sphere barycentric coordinate with i ∈ {1, 2, 3};

F vs = contact force on the virtual sphere.

Internal force model

The equations for calculating the internal forces in cylinders and PFacet elements are implemented in the YADE C ++ class Law2 ScGeom6D CohFrictPhys CohesionMoment ( ŠMILAUER et al., 2015b). This law is similar to the previous law (i.e., Law2 ScGeom -FrictPhys CundallStrack), but incorporates in its formulation adhesion as well as bending and twisting moments. Likewise, internal forces are calculated from relative displacements between the element GridNodes (i.e., Figure 4.7a). Internal moments are obtained analogously using linear laws (i.e., Figure 4.7b) which comprise relative rotations of the element nodes by the following equations: It is important to note that, since calculations are made in relation to the GridNodes, the deformation of cylinder and PFacet elements is governed by the connections of their nodes, known as GridConnections. This means that the parameters of these connections control the mechanical behavior of the element. Constitutive law parameters can be defined individually by the user or correlated to the mechanical and geometrical properties of the GridConnections using the following equations:

M r = k r • Ω r , (4.13) M tw = k tw • Ω tw , ( 4
k n = E c • π • R c 2 L c , (4.15
)

k s = 12 • E c • π•Rc 4 4 L c 3 , (4.16) k r = E c • π•Rc 4 4 L c , (4.17
)

k tw = G c • π•Rc 4 2 L c , (4.18)
where k n = normal stiffness; k s = tangential (i.e., shear) stiffness;

k r = bending stiffness;

k tw = torsional stiffness;

E c = Young's modulus of the GridConnection (i.e., the same Young's modulus of the first GridNode that characterizes the GridConnection);

G c = shear modulus of the GridConnection, that is, E c / [2 • (1 + ν n )]
, where ν n is the Poisson's ratio of the first GridNode that characterizes the GridConnection; R c = radius of the GridConnection;

L c = length of the GridConnection.

A failure condition is normally established by specifying elastic limits for the internal forces. These maximum values are related to specific adhesion values (i.e., Eqs. 4.19 and 4.20) that can be correlated in our case to the mechanical characteristics of the geogrid. It is important to emphasize that normal force is limited only under tensile conditions, since the compression at the contact is considered purely elastic. Normal and shear adhesion values are defined from the normal cohesion and shear cohesion parameters, respectively, which in turn can be used to modify the strength of the element. Thus, if the maximum tensile or shear force are reached, the cohesive links break and, consequently, the adhesion values become zero. In this situation, the connection between nodes will be lost and there will be no more internal interaction between them.

Plasticity condition is also possible in the model. This can be done by keeping the adhesion values. In this case, the behavior is perfectly plastic and the internal forces and also the moments are kept constant. The plasticity condition is defined in place of the failure condition by changing the boolean parameter fragile in the yade.wrapper.CohFrictMat and yade.wrapper.CohFrictPhys classes ( ŠMILAUER et al., 2015b). Regarding moments, the maximum bending and twisting values are defined in relation to the tensile strength of the element. The following conditions provide the limits for both conditions:

F n ≤ a n , (4.19) F s ≤ a s , (4.20) M r ≤ 4 • σ max n • I r R , (4.21)
M tw ≤ 2 • σ max n • I tw R , (4.22)
where Time-dependent behavior is implemented in Law2 ScGeom6D CohFrictPhys Cohe-sionMoment only for shear and twisting components. The creep behavior can be activated separately for each of these components. When both are activated, there is a viscosity incorporated in the model and the evolution of the elastic parts of shear displacement and relative twisting rotation is given by: FrictMat is used to characterize frictional materials (i.e., with contact friction). In this way, it is possible to define several material behaviors using different sets of parameters in the numerical modeling. Materials are assigned to geogrid elements and the parameters are used in the interactions to define the properties of the contact. According to the previous section, the internal behavior of the structure is calculated from the GridNodes. Therefore, in the case of the geogrid model, it is important to note that the contact parameters are applied only to GridConnections and PFacets.

F n =
du s,e dt = - F s ν s , ( 4
Contact parameters comprises measurable properties, whose values are used for interactions with bodies external to the geogrid (e.g., soil particles, gravel, and ballast). This set of parameters includes mainly the contact friction angle, but the interaction between soil and geogrid depends also on the dimension of the geogrid defined in the constructive process (i.e., length, width, thickness, aperture area, and cross-section of the longitudinal and transverse members).

Internal strength parameters

Internal parameters of the geogrid control its mechanical behavior. A cohesive-frictional material (i.e, with contact friction and cohesion) is defined by using the YADE C ++ class yade.wrapper.CohFrictMat and assigned to the GridNodes. In this case, the set of parameters comprises Young's modulus, Poisson's ratio, normal cohesion, and shear cohesion. Normal and shear cohesion were described in the previous section as parameters that respectively define the normal and shear strengths. Bending and twisting moments (i.e.,Eqs. 4.13 and 4.14) can be disregarded in interactions between elements defined with cohesive-frictional material. This is possible by changing the boolean parameter momentRotationLaw in the yade.wrapper.CohFrictMat class. This is convenient when geogrid connections must be completely hinged and can rotate freely with respect to each other.

Parameters are assigned to each body or element added to the model. Therefore, it is the combination of the values that globally characterizes the modeled geogrid. Another point to note is that the numerical model allows the insertion of different parameters for the bodies that compose the geogrid, which makes this approach very flexible to study different types of problems.

Geogrid simulations 4.4.1 Modeling remarks

The simulations were designed not only to demonstrate the numerical geogrid's response to classical load cases, but also to provide an insight of the many possibilities of study that this approach offers for the interaction of geogrids with granular materials. Concerning the numerical model, Figure 4.8 shows two possibilities for the numerical construction of the geogrid. The first one dicretizes the geogrid in segments with five GridNodes and four PFactes. The second one consists in using four GridNodes and two PFacets in each segment. Compared with two-PFacet configuration, four-PFacet configuration has six more elements per geogrid segment (i.e., one GridNode, two GridConnections, and two PFacets). This means that the number of elements in the model grows linearly with the number of segments that make up the modeled structure. In DEM, contact detection is one of the tasks with the greatest computational demand. It is possible to evaluate this demand by analyzing the complexity of the algorithm used in the model. In the case of YADE, the complexity of the contact detection algorithm is approximately O (n log n), by the Big Omicron notation [START_REF] Knuth | Big omicron and big omega and big theta[END_REF], where n is the number of model elements. For the four-PFacet structure, n is equal to the number of elements of the two-PFacet structure plus five times the number of segments considered for the geogrid. The larger the number of model segments, the greater the complexity of the algorithm, which can greatly increase the computational cost.

Despite the use of two PFacets per segment results in a lower computational cost, the four-PFacet configuration exhibits a better symmetry condition. The simulations in this work involved a geogrid modeled with fifty apertures totaling two hundred ninetysix segments (i.e., Figure 4.9). The simulations of this section also do not involve the representation of layers of granular material, thus the greater number of elements of the four-PFacet configuration does not result in a significant impact on the calculations durations of the cases studied. Therefore, the two geogrid segment configurations with the four PFacets and two PFacets for the numerical modeling were evaluated and compared in simple cases, so that, at the present time the comparison between the two is very minimalist. As regards parameters, fictitious values were assumed to simplify calculations, to intensify deformations for a better visualization in graphical outputs, and to increase computational efficiency during the simulations. Table 4.1 lists the parameters used in the analyses. It is important to note that the values used do not necessarily correspond to those of real geogrids and that the use of these parameters is not related to some limitation of the model.

Traction and shearing tests

Traction and shearing analyses were performed to verify how the geogrid behaves and to ensure that the combinations of the different types of elements that compose it do not result in unexpected behavior. These tests consisted of replicating a linear elastic behavior in simulations involving traction and shearing of the modeled structure. The results were compared with analytical solutions to validate the numerical model response.

= K N /K T ) 0.30 - Normal cohesion ∞ N/m 2 Shear cohesion ∞ N/m 2
The traction test consisted of pulling the geogrid axially at a constant velocity through its left end and calculating the axial force relative to the deformation produced by the traction. The shearing test allowed investigating the shear strain of the geogrid (i.e., shear tensile simulation). In this case, the reinforcement was pulled from the left end at a constant velocity in the transverse direction. As regards the boundary conditions, the structure was placed in a space without any gravitational force and fixed only at its right end. ∆L n and ∆L s are respectively the normal and shear displacements applied in the analyses and L o is the geogrid initial length. For a tensile elastic model with no cross section reduction phenomenon, the axial force can be calculated by:

F a = E • A cs • ε grid , (4.25)
where F a = axial force; E = Young's modulus;

A cs = the total yarn cross-sectional area;

ε grid = geogrid strain, that is, ∆L n /L o for normal tensile simulation and L o 2 + ∆L s 2 0.5 -L o /L o for shear tensile simulation.

Figure 4.10 shows the classical geogrid width reduction during the tensile test. For comparison with the analytical formulation, it is important to suppress this effect in the simulations, since it is not contemplated by Eq. 4.25. Numerically, this can be achieved by defining the stiffnesses of the diagonal GridConnections (i.e., connections required for PFacets) equal to zero or a negligible value. According to Eqs. 4.15 to 4.18, stiffnesses are directly proportional to the Young's modulus and to the radius. As the calculations are performed in relation to the GridNodes, the reduction of stiffnesses by the Young's modulus is not convenient in this case, since it would also affect the non-diagonal connections. Thus, a suitable way is to reduce the radius of the desired GridConnections. Considering geogrid modeling configuration with two PFacets per segment (i.e., Figure 4.8), it is possible to define the radius of the diagonal GridConnections equal to zero so that there is no influence of the stiffness of these elements in the model. Differently, in the configuration with four PFacets per segment, the internal nodes of the geogrid members make the model numerically non-convergent for zero-radius diagonal connections. Therefore, radius values close to zero can be used to obtain negligible stiffness.

Different effects are obtained by assuming a shear deformation of the geogrid without reducing the radius of diagonal GridConnections. Figure 4.11 shows a tendency of rotation of the transverse members during the imposition of transverse displacement at the left end of the geogrid. This effect, as in Figure 4.10, is caused by the constraint caused by the diagonal connections to the deformations imposed on the structure. Therefore, it can also be suppressed by reducing the stiffness of these connections Before analyzing the results of the numerical simulations (i.e., tensile and shearing tests), it is important to highlight the possibilities of the numerical model to describe complex geogrid behavior (e.g., tensile force and bending moments at the nodes). Since the nodes, longitudinal members, and transverse members of the geogrid are structures composed of GridNodes, GridConnections and PFacets, it is necessary to gather data from these elements to obtain localized information, such as forces and bending moments in nodes or longitudinal strain and transverse strain within particular member of the geogrid. For example (i.e., Figure 4.11), longitudinal and transverse strain of geogrid node can be determined by calculating the respective elongation and dividing it by the initial length. Since the elongation may be different in each node contour (i.e., GridConnections), the average elongation in the middle of the node can be calculated by averaging the values of the contours with the same direction. Node rotation can be calculated as a function of an axis passing through its center. Concerning the resultant force and bending moment acting on the node, they can be obtained by considering respectively the forces and moments for the central GridNode of the node structure. In order to exemplify this, graphs containing strain, resultant force, and bending moment, calculated for the detail node of Figure 4.11 during the shearing test, were plotted respectively in Figures 4.12,4.13,and 4.14. By convention, bending moment and rotation are positive clockwise. As demonstrated on these figures, without reducing the radius of diagonal GridConnections, bending moments, transverse strain, and transverse tensile force take place within the crossing nodes. This can be realistic but not suitable for comparison with analytical results. Figure 4.15 shows the graphical output of the normal and shear tensile simulations for the geogrid. Both simulations were performed assuming absence of the width reduction effect. In YADE, tensile force is computed relative to the GridNodes. For the simulations, this vector quantity was obtained by summing the resultant forces on the nodes in which the traction was imposed (i.e., first row of GridNodes at the left end of the geogrid). The tensile force in the shear traction simulation is the resultant in the direction of the longitudinal elements, it is calculated by the sum of its horizontal and vertical components. Figures 4.16 and 4.17

Bending test

Bending simulation aimed to evaluate the membrane effect of the geogrid under uniform loading (i.e., self-weight). Thus, this simulation consisted of fixing the geogrid by its two smaller sides in a given space and then imposing the gravitational force on the model. Boundary conditions were only defined for the left and right ends of the geogrid by fixing it in space. Strain and maximum deflection due to the self-weight of the structure were calculated numerically. Figure 4.18 illustrates the initial and final states of the simulation.

Concerning the analytical solution, [START_REF] Delmas | Sols renforcés par géotextiles -premières études[END_REF] proposed equations, later explored by [START_REF] Hello | Renforcement par géosynthétiques des remblais sur inclusions rigides, étude expérimentale en vraie grandeur et analyse numérique[END_REF], for the study of the membrane effect on geosynthetic materials subjected to a uniform loading. Assuming a horizontal geosynthetic sheet with linear elastic behavior and a uniform vertical load, the previous authors showed that the geosynthetic after loading had a parabolic deformation profile. They presented equations for the calculation of the tensile forces and maximum deflection. The calculation of these values is done in relation to the parameter β determined by the equation:

q • L o J = 3 • β • 1 + β 2 + arcsinh (β) -2 • β 3 + β 2 , ( 4.26)
where q = uniform distributed load;

L o = geosynthetic initial length; J = geosynthetic tensile stiffness.

For the bending simulation, the uniform distributed load, q, is approximately equal to 0.30 N/m 2 (i.e., no infinite width of the sheet). The geogrid tensile stiffness, J, can be obtained by the result of the previous simulation and is given per unit width by the relationship between the tensile force and deformation force shown in Figure 4.16 (i.e., J = 3.02E+02 N/m). The solution of Eq. 4.26 gives the value of β (i.e., β = 14, 38E-02), which leads to the determination of the horizontal component of the tensile force (i.e., T h = 102.89E-02 N/m given by Eq. 4.27), the maximum tensile force (i.e., T max = 103.31E-02 N/m given by Eq. 4.28) and the geosynthetic maximum deflection (i.e., d max = 3.60E-02 m given by Eq. 4.29). The maximum strain (i.e., ε max = 3.42E-03) is deduced of the maximum tensile force by Eq. 4.30.

T h = q • L o 2 • β , (4.27) T max = q • L o • (1 + β 2 ) 1/2 2 • β , (4.28) d max = β • L o 4 , ( 4
.29) [START_REF] Giroud | Determination of geosynthetic strain due to deflection[END_REF] introduced an approximation for a smooth deflection curve derived from equations of parabolic profile using truncated Taylor series expansion (i.e., Eq. 4.31). Assuming this approximation, the maximum geosynthetic strain can be estimated analytically by replacing in it the previously calculated deflection (i.e., ε max = 3.45E-03).

ε max = T max J , ( 4
ε max = 8 3 • d max L o 2 , (4.31)
Figures 4.19 and 4.20 show the comparison of the results obtained numerically (i.e., during time) and analytically for the bending of the geogrid. Numerical results converge fairly quickly to a stable solution through the use of damping. In Figure 4.19, the values of the maximum relative deflection (i.e., d max /L o ) and maximum strain (i.e., ε max ) are indicated respectively on the left and right axes of the graph. In terms of relative deflection, the result obtained by the numerical model is slightly higher than that obtained by the analytical equations (i.e., the difference between the two is 4.71 %). The explanation for this lies in the numerical discretization of the geogrid yarns: non-infinite grid width which leads to no identical values of the vertical displacements in the grid width and a limited number of segments in the length (i.e., the greater the discretization, the smaller the difference between the results and the greater the computational cost). In the case of the maximum tensile force and strain of the geogrid, the values obtained by Eq. 4.28 (i.e., 103.31E-02 N/m) and Eqs. 4.30 and 4.31 (i.e., 3.42E-03 and 3.45E-03, respectively) almost overlapped the result obtained by the simulation. Thus, it can be concluded that the results are consistent and demonstrate that the numerical model behaves properly.

Interation number

0.00e+00 2.00e+04 4.00e+04 6.00e+04 8.00e+04 1.00e+05 1.20e+05 0.00e+00 2.00e+04 4.00e+04 6.00e+04 8.00e+04 1.00e+05 1.20e+05 The numerical model can also be used to access data in different sections of the geogrid or even in its individual members. To exemplify this, tensile forces in a section 0.45 m from the left end of the geogrid were obtained in the bending simulation of the previous subsection. Figure 4.21 defines this section as well as the geogrid yarns in which the tensile forces were calculated. Considering the case of geogrid deflection by self-weight, the tensile force in different sections of the geosynthetic can be calculated from the following equations:

z(x) = β • x 2 L o - L o 4 , (4.32) ∂z ∂x = 2 • β • x L o , (4.33) T ϕ (x) = T h • 1 + ∂z ∂x 2 , (4.34)
where z(x) = deflection of the geogrid for a given section;

x = the position of the section relative to the center of the geogrid;

T ϕ (x) = tensile force in the predefined section (i.e., per unit width). Tensile forces in section S-S were computed by the numerical model for each of the six geogrid longitudinal yarns (i.e., non-identical values of tensile forces on the grid width due to boundary effect, videlicet, non-infinite grid width). Figure 4.22 shows these values. Higher tensile forces occur in the yarns closest to the longitudinal axis of the geogrid. As expected, the yarns symmetrical to the longitudinal axis show forces equal in magnitude observed by the overlap of the curves in Figure 4.22. The total tensile force per unit width in section S-S is given by the sum of the values obtained in the same section for each yarn. This sum is approximately equal to 102.90E-02 N/m which is the value of T ϕ (x) calculated for the section under analysis (i.e., x = 0.05 m).

A final simulation is presented below to illustrate the ability of the numerical model to take into account interactions between geogrid members and granular particles. For this purpose, the bending simulation was again used with the inclusion of two spheres side by side dropping them in the center of the geogrid. In the initial state, the spheres were positioned slightly above the geogrid at a distance of 3.56E-03 m. Fictitious parameters were adopted merely to simplify calculations. Thus, the same values of Table 4.1 were used for the geosynthetic and for the spherical particles the Young's modulus, Poisson's ratio, density, and friction angle were equal to 1.00E+05 N/m 2 , 0.30, 26.00 kg/m 3 , and 0.35 rad, respectively. The numerical model makes accessible the localized behavior in both the reinforcement structure and the interacting particles. The load is transmitted between the geogrid and the particles through the contact law described in the previous section. As presented earlier, when particles come into contact with the surface of GridConnections and PFacets, the forces are interpolated to the nodes that constitute the geogrid skeleton, which in turn transmit the forces internally to the structure. Since it is possible to numerically access the data of each geogrid element, this scheme makes it possible to study the behavior of specific parts of the geosynthetic.

Granular assembly

The arrangement and interlocking of the granular assembly can significantly influence its mechanical behavior. For example, [START_REF] Matsushima | Discrete element modeling for irregularly-shaped sand grains[END_REF] found greater shear strength in arrangements of elements that have more angular or non-convex shapes. The geometrical characteristics of the elements can develop particle roll constraints, which in turn makes the arrangement have greater shearing resistance. The bending moment generated in the interaction between non-convex elements also contributes to increasing the overall strength of the granular assembly. Thus, classical discrete models based on independent spherical elements may not properly reproduce the shear strength of, for example, triaxial tests on granular soils [START_REF] Salot | Influence of relative density on granular materials behavior: DEM simulations of triaxial tests[END_REF].

There are numerical studies dedicated to more realistically mimicking grain rolling. They can involve two approaches. The first one comprises the definition of bodies with non-spherical geometries [START_REF] Cundall | Formulation of a three-dimensional distinct element model -Part I. A scheme to detect and represent contacts in a system composed of many polyhedral blocks[END_REF][START_REF] Jensen | DEM simulation of granular media-structure interface: effects of surface roughness and particle shape[END_REF]LU;MCDOWELL, 2006) and the other one focuses on the integration of contact laws with rotation restrictions (IWASHITA; ODA, 2000). Although the first approach requires a higher computational cost, especially for the contact detection task, it has the advantage of not artificially modifying contact laws.

In the model of this thesis, rigid aggregates of individual spheres numerically recreate the non-convexity feature of granular materials. A non-convex volume occurs when a line segment formed by two internal points crosses the volume domain (MORRIS; STARK, 2015). For this purpose, the model used clumps elements. Clumping spheres together, without relative displacement between these spheres during simulation, makes it possible to model different particle geometries. It is important to note that the intention is not to model any specific particle but to obtain non-convex elements that can capture the interlocking effect.

The number of spheres per clump can produce different macroscopic results in, for example, triaxial test simulations since the element shape is directly related to its rolling mechanisms (i.e., perfect symmetry for one-sphere clumps, rotational symmetry for twosphere clumps, and asymmetry for three-sphere and four-sphere clumps). From a certain asymmetry (e.g., three-sphere or four-sphere clumps), the difference in the packing mechanical behavior for a larger number of spheres becomes less significant [START_REF] Salot | Influence of relative density on granular materials behavior: DEM simulations of triaxial tests[END_REF]. Therefore, aiming at particles capable of reproducing the behavior of granular materials more realistically and with less influence of geometric aspects (e.g., symmetry), this thesis used four spheres per clump.

Considering a sphere of radius r, the ratios used for the spheres sizes of each clump were r, 0.83r, 0.67r and 0.67r (i.e., two spheres of equal size and two spheres of different sizes). The numerical granular sample randomly combined clumps of the same geometry but with different sizes to avoid a regular packing in terms of elements. Four clump sizes were modeled by increasing and also reducing the radii of the spheres by 20 % and 40 % so that the ratio between the largest and the smallest clumps was 2.5. Figure 4.24 shows the geometry of the typical clump used in the numerical analyses. It is possible to evaluate the concave aspect of the clumps by a defined coefficient of angularity based on the radii and the distances between the component spheres, according to the equation:

ang = 2! (n -2)! n! 1≤i<j≤n d i,j R i + R j (4.35)
where ang = angularity coefficient; n = number of spheres;

d i,j = distance between sphere centers; R i , R j = spheres radii.
Null angularity means that the centers of spheres coincide, and therefore the overlap is complete, resulting in a perfectly spherical particle. Values greater than unity represent clumps created with non-contacting spheres. Figure 1 illustrates the concept of the angularity coefficient for clumps composed of two spheres. The concave aspect can play a significant role in the interlocking between the elements, especially at high inter-particle friction angles. Although angularity can modify the rolling mechanisms of the elements, and therefore the packing porosity, [START_REF] Salot | Influence of relative density on granular materials behavior: DEM simulations of triaxial tests[END_REF] found that the difference between the maximum and minimum porosities of granular samples remains relatively constant from a very low angularity coefficient (e.g., greater than 0.2). Thus, it is possible to treat the difference between the minimum and maximum porosities of a sample independently of the geometry of its elements. Because of that, it is possible to use the concept of relative density (i.e., as a function of the minimum and maximum packing porosity) to characterize numerical samples. Table 4.2 shows the angularity coefficient for the four-sphere clump used in this work. TABLE 4.2 -Angularity coefficient for the typical 4-sphere clump.

Sample preparation

The sample preparation consisted of obtaining different geogrid-reinforced granular volumes for given confinement before the geogrid pullout. The methodology for obtaining controlled compacted samples consisted of the following three main steps:

• To introduce elements into the model;

• To manage porosity and consolidate the sample;

• To get ready for pullout simulation.

The first step comprises generating the geogrid and the granular sample and defining their respective properties, as described in the previous sections. The geogrid is initially arranged fully stretched in the center of the sample between two layers of granular material within a box. The model generates the box by utils.aabbWalls() function of the yade.utils module ( ŠMILAUER et al., 2015b). This function uses AABB to insert into the model six walls that surround the sample from all sides. The first two steps impose a static condition on the geogrid, making it impossible for any movement or deformation. Another important point to note is that the model does not take into account the interactions between the box and the geogrid, so the only constraint exerted on reinforcement is due to the granular layers. This is also why it is possible to pull out the geogrid through the wall during the simulation.

Radius Expansion -Friction Decrease -REFD is a procedure for obtaining a dense granular assembly with specific porosity (CHAREYRE; [START_REF] Chareyre | Discrete element modeling of curved geosynthetic anchorages with known macro-properties[END_REF][START_REF] Salot | Influence of relative density on granular materials behavior: DEM simulations of triaxial tests[END_REF]. Thus, the second step involves applying this method to the sample. Initially, it consists of increasing the particle radius of the granular material while controlling the pressure within the sample boundaries to a preset value, which results in an initial porosity value for the sample. Subsequently, the sample achieves the desired porosity during the application of consolidation steps by simultaneously reducing the inter-particle friction angle and increasing the radii to maintain the confining conditions. Importantly, there is no friction between the box and the granular assembly in the model of this thesis. Thus, the walls are only to ensure the confinement of the sample.

The third and last step comprises some final sample settings for geogrid pullout simulations. After obtaining the desired porosity for the confined sample, it is necessary to recover the inter-particle friction angle to the initial value defined for the model or to any other value in case of behavior studies based on this parameter. Removal of the static condition of the geogrid imposed in the first two steps is also necessary to ensure the degrees of freedom of geosynthetic reinforcement. Finally, the sample is ready for simulation when the system reaches a static equilibrium condition after imposing the gravitational force on the model. It is useful to apply vertical color bands to particle groups to improve motion visualization during the simulation. Figure 4.26 shows two examples of samples samples during preparation, one with six thousand spheres and one with six thousand clumps. Figure 4.27 provides the same samples ready for the geogrid pullout test. 5 Numerical Analyses and Results

Geogrid pullout test

Geogrid pullout tests involved granular assembly samples with both spheres and clumps. In these two cases, the samples comprised three thousand, six thousand, and nine thousand particles, totaling six models (i.e., three using spheres and three using clumps). Sample preparation took into account all the steps described in Section 4.6. In the case of sphere packing, the model assumed uniform particle size distribution by defining the mean radius and standard deviation equal to 1.20E-02 m and 6.00E-03 m, respectively. Regarding clump packing, the random arrangement of rigid aggregates within the box comprised equal probability between the different sizes defined in Section 4.5. In the context of this thesis, the analyzes comprise samples prepared strictly in the same way and with the same parameters. Therefore, different states of consolidation between samples result from the number and geometry of the particles in each sample.

The preparation procedure assumed a confining pressure of 1.00E+03 N/m 2 . The defined boundary conditions do not take into account the friction between the container walls and the particles of granular material (i.e., Section 4.6). Aiming at achieving maximum sample consolidation for the applied confining pressure, the smallest possible porosity for the geogrid-reinforced granular assembly was the target during the REFD method. It is important to emphasize that the samples reach the minimum porosity by reducing in stages the inter-particle friction angle to zero during the consolidation step. As mentioned in the previous chapter, there is the restoration of the original friction angle value and application of gravitational force before the start of geogrid pullout tests.

After applying the gravitational force to the model, it is necessary to wait for the new static equilibrium condition of the sample. During this process, the particle self-weight increases contact overlaps, which represents higher contact forces. Thus, there is stress relief on the upper wall of the box, and stress increases at the other boundaries. In the simulations of this thesis, the deformations of the contacts caused the granular material to lose contact with the upper wall of the box. Therefore, there is no pressure applied to the top of the box at the beginning of the pullout tests.

Because box volume is invariable, more spheres or clumps in the model result in smaller particles. Therefore, the minimum porosity mentioned above is different between samples. This fact provides different dense samples whose results were analyzed and compared with each other. Table 5.1 shows the porosity, the smallest, largest, and average particle size values, and the ratio of average particle size to geogrid mesh size for each sample (i.e., samples with 3000, 6000, and 9000 spheres and samples with 3000, 6000, and 9000 clumps). The samples have an initial porosity characteristic of the first step of the REFD method and a minimum porosity (i.e., for the imposed conditions) that was the desired porosity for the assemblies. Particle sizes (i.e., minimum, maximum, and average) provide the length to circumscribe the particle, which is the diameter itself for spherical particles and the longest length between the component elements for rigid aggregates.

The opening of the geogrid mesh is 6.40E-03 m 2 . Considering the average particle size (i.e., Table 5.1), the estimated number of particles within a geogrid aperture for samples with 3000, 6000, and 9000 spheres are 4.21, 6.61, and 8.59 particles, respectively. For samples comprising 3000, 6000, and 9000 clumps, these values are 2.51, 3.96, and 5.15, respectively. Particle and mesh size ratio 6.95E-01 5.53E-01 4.85E-01 m/m

The average particle size in samples with 6000 particles is about 20 % smaller than in samples with 3000 particles in both cases involving spheres and clumps. For 9000 particle samples, this value is 30 % lower than 3000 particle samples and approximately 12 % smaller than 6000 particle samples. These proportions are associated with the minimum porosity values obtained for each model. It is important to note that the initial porosity does not influence the particle sizes of the prepared samples, so there is no perspective of proportionality between these values. All samples comprised the numerical model developed for true biaxial generic geogrids (i.e., the same mechanical properties in both machine and cross-machine directions) and described in Sections 4.2 to 4.4. The numerical construction of the geogrid involved the two-PFacet configuration since, for the case under study, it has lower computational cost and similar behavior regarding the four-PFacet configuration. The geogrid modeled for the samples consisted of fifty openings (i.e., 10 openings × 5 openings) and two hundred and ninety-six segments.

The geogrid pullout test analyses in this work are not intended to reproduce a real or experimental case. The central point in this thesis is to increase understanding of the interaction behavior between granular material and a generic geogrid. Studies concern how grains roll in reinforcement and how forces are mobilized locally on the longitudinal and transverse members of the geogrid. Thus, the models focused on capturing the main attributes of real generic geogrids (i.e., stretching, bending, shearing, and twisting) during their interaction with granular material (i.e., interlocking, friction, and abutment). Considering this, Table 5.2 presents the assumed parameters.

It is possible to determine the macroscopic parameters by simulating triaxial tests on granular material samples. Thus, obtaining these parameters consisted of performing simulations similar to Section 3.3.6 (i.e., non-periodic samples in this case) in the samples of granular material without geogrid. These simulations consisted of applying vertical loading while maintaining the confining stress of 1.00E+03 N/m 2 to the samples prepared as described in Section 4.6, but without the presence of the geogrid, until the maximum measured axial stress value exceeded the current value by ten percent. This procedure ensured triaxial tests on granular samples with the same characteristics used in geogrid pullout simulations. Figure 5.1 shows each sample at the end of the triaxial test. Figures 5.2 to 5.7 show the deviator stress and volumetric strain as a function of axial strain for the triaxial tests performed. Table 5.3 presents the macroscopic parameters obtained for the granular samples.

From the densities of Table 5.3,gravitational acceleration (i.e.,9.81 m/s 2 ), and half of the box height (i.e., 0.25 m according to Table 5.2), it is possible to calculate the vertical stresses at the interface of granular material and geogrid. Considering the mean between the calculated values for the three samples comprising spheres and for the three samples comprising clumps, the vertical stresses at the geogrid interface are 3.60E+03 N/m 2 and 4.19E+03 N/m 2 , respectively. Geogrid pullout simulations essentially consisted of imposing constant velocity dis-placement on the left end of the reinforcement to partially pull it out of the prepared sample. Usually, pullout tests promote geogrid displacements until the first resistance peak is defined. However, the studies of this thesis included the partial pullout of almost half of the geogrid embedded in two layers of granular material. The simulations assumed 0.44 m (i.e., reinforcement length equal to 1.00 m) pullout from the sample, also intending to evaluate the behavior in large deformations. To avoid the occurrence of dynamic effects, the simulations assumed a low pullout rate of 4.69E-08 m per time step. The simulations also adopted a non-viscous damping coefficient (i.e., YADE C ++ class yade.wrapper.NewtonIntegrator) of 0.60 to dissipate kinetic energy between particles. The tests provided the pullout force on the reinforcement as well as its displacement. Functions developed for the numerical model also tracked stresses in longitudinal and transverse members of two geogrid openings during simulations. These particular openings are on the longitudinal axis and equidistant from the center point of the geogrid. Figure 5.8 shows these positions, identifies the respective members or elements evaluated, and also indicates the geogrid pullout direction. Cells 1 and 2 characterize the mentioned geogrid openings. In these cells, the elements identified by A and B refer to the longitudinal members and by C and D to the transverse members. Poisson's ratio 0.28 0.30 0.28 -FIGURE 5.8 -Openings and their elements evaluated in the geogrid pullout test.

Sample with spheres

This section covered three geogrid-reinforced samples prepared under the same confinement pressure and minimum porosity and containing 3000, 6000, and 9000 spheres, respectively. The application of constant displacement at the left end of the geogrid imposed a partial removal of the reinforcement embedded in granular material (Figures 5.9 to 5.11). Importantly, the model did not include interactions between the box and the geogrid, so that the walls did not restrict geosynthetic movement. Figures 5.9 to 5.11 show the final state of the geogrid pullout tests on spherical granular material. The strong and weak force chains characterize the force chain network structure (MC-DOWELL et al., 1996;[START_REF] Peters | Characterization of force chains in granular material[END_REF]. The strong force chain follows the main principal stress of a granular medium. The weak force chain forms around the strong chain and acts as auxiliary support. Loading conditions cause changes in the contact points of the particles, which in turn causes the chain force to rearrange to maintain system stability [START_REF] Wen | Grain size effect on the mechanical behavior of cohesionless coarse-grained soils with the discrete element method[END_REF]. Figures 5.12 In the force chain network, the lines connect the centers of the particles in contact and characterize a force distribution structure. The thickness of the lines in these diagrams is proportional to the magnitude of the total contact forces. Comparing Figures 5.12 to 5.14, it possible to verify that for smaller particle sizes (i.e., larger number of spheres in sample results in smaller particles), the distribution of contact forces in the sample is more uniform but with lower values (i.e., thinner lines of force chains).

In Figures 5.12 to 5.14, it is possible to notice the absence of a force chain of the granular material with the upper plate at the right end of the sample. Figures 5.9 to 5.11 also show this lack of contact in this part of the sample. With the displacement of the end of the geogrid within the sample, the particles rearrange themselves in the empty spaces, which results in this condition. This aspect and the force chain with the upper plate on the left side of the box (i.e., imposing side of the geogrid pullout) indicate that the vertical stresses at the geogrid level may vary over the tests. For the final state of samples with 3000, 6000, and 9000 spheres (i.e., Figures 5.9 to 5.11), the total normal force applied to the top plate of the box is 1.91E+03 N , 2.36E+03 N , and 2.43E+03 N , respectively.

Other information obtained through the simulations was the force-displacement curve. Figure 5.15 shows this type of curve for pullout force (i.e., distributed force at the left end of the reinforcement) and displacement relative to the end of the geogrid. Therefore, in the results, the forces initially grow with zero displacement due to the elastic deformation of the reinforcement. This deformation occurs in all three samples up to a pullout force of approximately 2.50E+03 N/m.

From the results, it is possible to distinguish the initial and secondary peaks in the curves. The peaks represent the mobilization of the reinforcement strength during its pullout of the sample. The initial resistance mobilization observed in the first peaks of the graph is mainly due to the friction between the granular material and the reinforcement. This difficulty in moving the geogrid promotes its extensional strain governed by its mechanical properties to a limit where the pullout force is greater than the static friction. As the geogrid moves, it is possible to observe subsequent peaks that essentially characterize the difficulty of movement imposed by the geogrid's transverse yarns (i.e., the transverse members can act as a buttress). The space between the geogrid transverse yarns is 0.08 m, which is approximately the value of the distance between the curve peaks.

The confinement provided by the granular material and the presence of particles within the mesh openings result in interlocking mechanisms, which in turn mobilize the reinforcement when there is relative displacement between elements. The pullout force is proportional to the length of the reinforcement embedded in the granular sample. During the test, the pullout force decreases (i.e., the tendency of lower values) because the length of the geogrid inside the box also decreases. When comparing the development of the force-displacement curve obtained for the samples, it is possible to notice that the initial peak decreases for a larger number of spheres in the sample (i.e., smaller particle and mesh size ratio), although the observed difference is relatively small. In this case, there is also a tendency to reduce the initial slope of the curve with decreasing sphere size. Secondary peaks are more pronounced in the sample with 3000 spheres. As the number of spheres increases, subsequent peaks tend to a constant value (i.e., residual pullout force). The number of particles within the geogrid mesh also influences the results. The approximate number of spheres within the geogrid apertures (i.e., based on average particle size) is 4.21, 6.61, and 8.59 for samples with 3000, 6000, and 9000 spheres, respectively.

Fewer spheres lead to larger diameters, hence the rolling constraint is greater on these larger mass particles, which in turn results in greater pullout resistance. Also, in wellconsolidated samples with interlocked particles, the smaller the particle, the less expansion between contacts as the particle moves. This small disturbance compared to that caused by larger particles in the same confining conditions results in less force mobilization in the geogrid. Furthermore, as the particle diameter increases, the shear band size also increases, which results in higher pullout force. This explains the observed behavior for the initial peak as well as the lower secondary peaks with the tendency to a constant pullout force as the particle size becomes smaller in the sample. Another important point in these results is that, although the minimum porosity obtained is smaller in the samples with a larger number of particles (i.e., Table 5.1), this increase in the number of spheres does not have a significant impact on the value of the initial peak strength.

Regarding the axial stresses in the longitudinal members of the geogrid openings (i.e., elements A and B in Figure 5.8), the values obtained in the longitudinal elements of cell 1 are always higher than in cell 2, as shown in Figures 5.16 and 5.17. The distribution of the anchoring forces explains this behavior. That is, the greatest efforts are concentrated in the pullout force application region and tend to zero at the opposite end of the geogrid. The axial stress values observed in cell 1 elements A and B for the 3000-sphere sample are slightly lower than those for 6000 and 9000 spheres, although cell 2 does not openly exhibit the same behavior. Except for this point, the stresses in elements A and B presented average values very close between the samples (i.e., Table 5.4). A specific arrangement of the spheres during geogrid pullout, which results in greater confinement of the reinforcement members, can explain inversions of the values as well as exceptional peaks obtained in the longitudinal elements during simulations. Figures 5.18 and 5.19 show for both cells 1 and 2 the axial stresses of the transverse elements identified respectively as B and C. In the transverse elements, differently from the longitudinal elements, the values of the axial stresses in cells 1 and 2 follow the same range of values. This behavior implies that the main interaction mechanism of these elements is similar to that of a buttress that restricts the movement in the longitudinal direction of the reinforcement since they work as obstacles to the rolling of particles that occupy the space between the geogrid openings. Moreover, the mesh formed by the existence of the transverse elements allows the interlocking of the particles in these openings, which intensifies the restriction of horizontal movement between the granular material and the reinforcement. Figures 5.20 to 5.22,which show the frontal view of the geogrid pullout, emphasize this behavior as well as suggest that the axial stresses in the transversal yarns are mainly due to bending mechanisms. It is possible to notice in the last transverse yarn that leaves the sample the members bending due to interactions with the particles.

The mean value calculated considering the six stress-displacement curves of the transverse elements indicated by the red line in Figures 5.18 and 5.19 (i.e.,2.05E+06 N/m 2 and 2.07E+06 N/m 2 for elements C and D, respectively) is smaller compared to those of the longitudinal elements, where the resistive component due to friction is higher. Neverthe-

Sample with clumps

Similar to the previous section, it was possible to evaluate the behavior in samples with 3000, 6000, and 9000 clumps. These models considered the same confining conditions and also the preparation under minimum porosity. The application of the same pullout rate used in the spherical particle samples imposed the equivalent partial geogrid removal. The boundary conditions were also the same in all simulations (i.e., Sections 4.6 and 5.1). Figures 5.23,5.24,and 5.25 show the final state of geogrid pullout tests on non-spherical granular material containing 3000, 6000, and 9000 clumps, respectively. Concerning the force-displacement curve shown in Figure 5.26 (i.e., distributed force at the left end of the reinforcement versus displacement relative to the end of the geogrid), forces initially increase with zero displacement due to the elastic deformation of the reinforcement. This occurs up to a pullout force of 2.50E+03 N/m for the 3000-clump sample and 3.20E+03 N/m for the 6000-clump and 9000-clump samples. In the obtained curves, it is also possible to distinguish initial and secondary peaks. The three samples characterized the initial peak with two slope sections of the curve. The displacement value of 0.03 m defines the end and start of the first and second slope sections, respectively. The first slope refers to the mobilization of the reinforcement strength by the friction between the granular material and the geogrid yarns. The results showed a minimal influence of particle size on this part of the curve. In this part, the force-displacement curve tends to a level that defines the maximum frictional mobilization. However, due to the particles in the geogrid mesh, the curve follows a new slope until it reaches the initial pullout peak. As mentioned previously, the particles within the mesh interact with the geogrid transverse elements and restrict the movement of the reinforcement. The results indicated that the initial peak is higher for the sample with the largest clump size (i.e., sample with 3000 clumps). However, samples with 6000 and 9000 clumps showed the initial peak with very close values (i.e., 2.00E+04 N/m). In this case, the 9000-clump sample requires greater geogrid displacement to mobilize the same strength as the 6000-clump sample. It is important to emphasize that the average value of the force-displacement curve (i.e., the average value for all points on the curve) obtained for the 6000-clump sample is lower compared to the other two samples (i.e., 1.73E+04 N/m for 3000 clumps, 1.54E+04 N/m for 6000 clumps, and 1.67E+04 N/m for 9000 clumps). Despite this unexpected behavior, taking into account the extremes (i.e., 3000-clump and 9000-clump samples), it is possible to observe the tendency of smaller pullout force peaks in samples with smaller particles or with a smaller particle and mesh size ratio (i.e., samples with larger amount of particles). Sample with 9000 clumps (i.e., fewer particles) showed lower secondary peaks, although the residual pullout force (i.e., curve tendency to a constant value) in this sample was higher than in the 6000-clump sample. As the size of the particles in the same sample varies, it is possible that this characteristic caused, for geogrid displacement greater than 0.10 m, an arrangement of the particles in the sample with 6000 clumps favorable to the rolling mechanisms, which resulted in lower stresses at the interface between the granular material and the reinforcement and thus the tendency for smaller residual pullout force compared to the sample with smaller particle size (i.e., 9000-clump sample). Regarding the peaks of the curves, it is also possible to notice that the distance between them is approximately 0.10 m, a value close to the distance between the transverse yarns of the geogrid (i.e., 0.08 m).

The force chain network created by the particle contacts plays a substantial role in the mechanical behavior of the samples. In the case of the geogrid pullout test, the distribution in the system of stresses caused by the imposition of reinforcement displacement occurs through the contacts. Therefore, the force chain network varies throughout the test, adjusting according to the arrangement of system elements. This adjustment is more uniform in samples with smaller particles, due to contacts networks that are also more uniform in these samples. The more uniformity, the more branches of the force chain, which makes the weak force chain more present in the system. A weak force chain, as its name implies, concentrates lower contact forces than a strong force chain. Consequently, the reinforcement force mobilization is lower in samples with smaller particles. Figure 5.26 indicates this behavior for samples with 3000 and 9000 clumps. Figures 5.27 to 5.29 show the force chain network for the samples composed of non-spherical granular material in the final state of the simulation.

In the case of the 6000-clump sample, the larger particle size results in less uniform force chain than in the 9000-clump sample, which justifies the more notable secondary and subsequent peaks rather than a tendency for a residual value. Nevertheless, its lower mean value of the force-displacement curve (i.e., 1.54E+04 N/m for 6000 clumps) suggests that the reasoning of contact forces is not straightforward in samples with non-spherical particles, indicating the importance of future work on this topic.

As additional information, the 9000-clump sample best characterized high and low force concentration zones. It is possible to see in Figure 5.29 the force lines forming at least two shells facing the geogrid pullout direction and centered on the reinforcement position. The convex side of these shells characterizes the zones of high concentration of forces while the concave side characterizes the zones of low pressure. This aspect of the force chain network is due to the geogrid transverse yarns that interact with the granular material to produce a shell-like force distribution. The concavity of these force distribution shells in the granular medium towards the pullout direction illustrates the idea of a parachute-like mechanism resisting the imposed motion. Figures 5.27 to 5.29,as well as Figures 5.23 to 5.25, show a lack of contact between the granular material and the upper plate at the right end of the sample. In contrast, there is a force chain with the upper plate on the left side of the box (i.e., imposing side of the geogrid pullout). Therefore, vertical stresses at the geogrid level may vary over the tests. For the final state of samples with 3000, 6000, and 9000 clumps (i.e., Figures 5.23 Concerning the axial stresses in the longitudinal members of the geogrid openings (i.e., cells 1 and 2), Figure 5.30 and 5.31 show the stress-displacement curves for elements A and B, respectively. The higher values obtained in the longitudinal elements of cell 1 compared to those of cell 2 are due to the distribution of the anchoring forces, which is greater the closer to the point of application of the pullout force (i.e., the left end of geogrid). relative to the residual stress. After the first peak of axial stress, the values tend to stabilize at the same level, except for cell 2 elements A and B of the 3000-clump sample that showed a slight curve ascendancy throughout the test. This behavior suggests that the interlocking caused by non-spherical particles keeps the confinement at the interface of granular material and geogrid, resulting in a maintenance of stresses even after large deformations.

The stress-displacement curves do not indicate an abrupt drop after the initial peak

The axial stress values observed in cell 1 elements A and B for the 6000-sphere sample are slightly lower than those for 3000 and 9000 spheres. This result is consistent with the behavior shown by the force-displacement curve of the geogrid pullout test (i.e., Figure 5.26), although the expectation was an intermediate stress-displacement curve compared to the other two samples.

Importantly, geogrid cell 2 elements A and B did not exhibit the same behavior observed in the cell 1 elements. The 3000-clump sample had the lowest initial peak, although the secondary and subsequent peaks are relatively higher than the other samples. For the samples with 6000 and 9000 clumps, the computed values were very close. This aspect again indicates the importance of future work involving contact forces in geogrid-reinforced layers with non-spherical particles. Table 5.5 shows the mean values for the stress-displacement curves of the three samples (i.e., elements A and B of cells 1 and 2). The 6000-clump sample presents slightly lower mean values compared to the other two samples, which is in agreement with the previously analyzed results. The sample with larger particles (i.e., 3000-clump sample) had the highest mean values for axial stress in the analyzed elements (i.e., 5.62E+07 N/m 2 and 5.59E+07 N/m 2 respectively for elements A and B of cell 1 and 2.59E+07 N/m 2 and 3.08E+07 N/m 2 respectively for elements A and B of cell 2). Despite that, the stresses in the elements of the respective cells indicated very close average values between the samples.

In the transverse elements, the values of the axial stresses in cells 1 and 2 have the same range of values. Figures 5.32 and 5.33 show for the two cells the axial stresses of the transverse elements B and C, respectively. This behavior is consistent with the interaction mechanism of the transverse elements that act as obstacles to the rolling of particles present inside the geogrid openings. Therefore, in this case, the stress values in the transverse elements of geogrid cells at different positions should not differ significantly. Figures 5.34 to 5.36 show the frontal view of the geogrid pullout tests. The bending of the transverse members of the yarn leaving the sample suggests that the axial stresses are mainly due to bending mechanisms.

The mean value calculated for the stress-displacement curves of the transverse elements (i.e., considering the three samples) indicated by the red line in Figures 5.32 and 5.33 (i.e.,3.60E+06 N/m 2 and 4.19E+06 N/m 2 for elements C and D, respectively) is smaller compared to those of the longitudinal elements. During the geogrid pullout, the friction mobilization in the transverse elements is lower than in the longitudinal elements, which favors the occurrence of also lower stresses in these elements. The stresses mobilized on the transverse elements are not negligible, including the first stress peak that is close to the average value calculated on the element throughout the test. For samples with clumps, the average value obtained for the stress-displacement curves of transverse element C (i.e., 3.60E+06 N/m 2 ) is very close to that of element D (i.e., 4.19E+06 N/m 2 ). The average value obtained together for the transverse elements C and D (i.e., 3.89E+06 N/m 2 ) is approximately 8 % and 14 % of the average value calculated for the longitudinal elements of cells 1 (i.e., 5.16E+07 N/m 2 ) and 2 (i.e., 2.71E+07 N/m 2 ), respectively. Importantly, stress peaks in the transverse elements during geogrid pullout can exceed three times the calculated average value.

Comparison of results

Figure 5.37 compares the force-displacement curves for samples with the same number of particles, but with different geometries (i.e., spheres and clumps). As the measured displacement is relative to the end of the geogrid, there is a record of forces at the beginning of these curves up to a limit value at which the end of the reinforcement begins to move. This behavior is due to the initial elastic deformation of the geogrid. In samples with spherical particles, this limit is 2.55E+03 N/m, 2.60E+03 N/m, and 2.42E+03 N/m for samples with 3000, 6000, and 9000 particles, respectively. In the samples comprising clumps, these forces are respectively 2.53E+03 N/m, 3.20E+03 N/m, and 3.21E+03 N/m. Except for samples with 3000 particles, samples containing clumps indicate greater elongation of the geogrid due to the greater initial pullout forces. Higher forces imply higher shear stresses at the interface of the geogrid with the granular material, which in turn can intensify the dilatancy (i.e., macroscopic behavior) of the medium during the pullout.

By comparing the results obtained for both sphere and clump samples (i.e., Figure 5.37), it is possible to verify that the maximum pullout force characterized by the initial peak of the curves occurs for a larger reinforcement displacement in the samples with clump particles. Although these samples indicate an initial pullout strength approximately twice that of spherical particle samples (i.e., the first force peak in Figure 5.37), a greater geogrid displacement is required to mobilize this condition.

The force-displacement curve defines two different slopes to the first peak in clumpparticle samples (i.e., the displacement value of 0.04 m defines approximately the end of and start of the first and second slope sections, respectively, for samples comprising 3000, 6000, and 9000 clumps), which indicates two types of force mobilization. The first one is mainly due to the friction between the granular material and the geogrid yarns. In this case, the non-spherical particles provided greater force in the geogrid, which can be verified by the greater slope of the initial portion of the curves involving clumps compared to those involving spheres. Possibly spherical particles allow fewer points of contact with the reinforcement, resulting in a lower frictional force component. The second one is a result of the imbrication (i.e., interlocking) of non-spherical particles that hinder the rolling mechanism and therefore require greater force during the geogrid pullout. In this type of force mobilization, different arrangements of the particles that interact with the reinforcement can occur and result in a change in the slope of the force-displacement curve. A similar phenomenon may have occurred in samples with 3000 and 9000 clumps, as Figures 5.37a As the number of particles increases (i.e., a decrease in the particle and mesh size ratio), the initial and subsequent peaks of force decrease, although the 6000-clump sample showed well-characterized subsequent peaks. It is also possible to verify in Figure 5.37 the tendency to a residual pullout force with the increasing number of particles in the model. In the case of 9000-particle samples, the clump-particle sample indicated a tendency to maintain pullout strength near the initial peak level rather than an abrupt drop in value. The interlocking effect in non-spherical particles seems to keep the same level of confinement at the interface of granular material and geogrid even after large deformations. Indeed, the difference between the initial and subsequent peaks in Figures 5.37a and 5.37b is smaller in clump-particle samples.

As discussed in Sections 5.2 and 5.3, there is a concentration of contact between the granular material and the upper plate in the first part of the box. This condition provides greater confinement of the sample in this region, which results in greater vertical stresses. Considering the final state of the samples, the total normal force applied to the top plate of the box is approximately 2.23E+03 N for samples with spheres (i.e., average between 1.91E+03 N , 2.36E+03 N , and 2.43E+03 N calculated for samples with 3000, 6000 and 9000 spheres, respectively) and 5.65E+03 N for samples with clumps (i.e., average between 5.85E+03 N , 5.19E+03 N , and 5.92E+03 N calculated for samples with 3000, 6000 and 9000 clumps, respectively). This difference in the value of approximately 2.50 times greater for non-spherical samples compared to spherical samples is due to the dilatancy that is greater in the clumps assembly.

Regarding the axial stresses in longitudinal elements A and B, Figures 5.38 and 5.39 compare their stress-displacement curves in both cells 1 and 2 for samples with the same number of particles (i.e., spheres and clumps). Although these results were previously presented separately in the previous sections, Figures 5.38 and 5.39 allow more direct visualization of important aspects of sample behaviors. One aspect is that the ratio between the averages of calculated stress for a specific element in cell 1 and cell 2 is very similar between samples. This ratio is close to 2.00, as shown in Table 5.6. Another aspect is that the axial stresses ratio between corresponding geogrid elements of clump-particle sample to sphere-particle sample resulted in approximately equal values in samples with the same number of particles (i.e., Figures 5.38 and 5.39). For samples with 3000, 6000, and 9000 particles, the average ratio is 2.75, 1.92, and 2.43, respectively. Table 5.7 presents the calculated ratios. Although further studies are needed, this linear principle of behavior may be worth studies of shape-matching factors for simulations with Except for Figure 5.38c, the results in Figures 5.38 and 5.39 show a better characterized initial peak in the sphere-particle samples. In samples with clumps after the first peak, the stresses on elements A and B show a reasonable propensity of a residual value equal to the peak stress. The results indicate that the higher uniformity of a medium with a larger number of particles (e.g., Figure 5.38c) or the interlocking of non-spherical particles (e.g., clump-particle samples) favored this behavior.

Finally, concerning the stress-displacement curves in the transverse elements, the ratio of the mean value calculated in all samples in non-spherical granular medium to spherical granular medium is 1.76 for element C and 2.02 for D element (i.e., Figures 5.18,5.19,5.32,and 5.33). Therefore, the increase in calculated stresses for samples with clumps is approximately double those calculated for samples with spheres. Importantly, the geogrid displacement required to mobilize the initial stress peak in these geogrid elements on both types of granular material is very close, requiring approximately 0.02 m of geogrid displacement.

When comparing the results, it is possible to verify that the increase in the number of particles minimally influences the interaction mechanism of the granular material with the geogrid transverse yarns (i.e., overlapping stress-displacement curves). On the other hand, the particle geometry of the granular medium plays a fundamental role in this mechanism of interaction. As assessed, the presence of non-spherical particles in the model approximately doubled the axial stresses in the geogrid elements C and D. Therefore, the angularity of particles is a significant point to consider in geogrid-reinforced layer solutions.

As a final aspect, it is possible to compare the ratio between the value calculated by the pullout force divided by the geogrid surface area (i.e., total area) and the vertical stress at the geogrid interface. Table 5.8 presents the data and values of this ratio for samples comprising spheres and clumps. The ratios calculated for samples with clumps are, on average, 65 % higher than those calculated for samples with spheres. This difference between the values is mainly due to the higher pullout forces obtained in the tests performed on the samples with non-spherical particles. The calculated ratios for the geogrid (i.e., Table 5.8) are greater than those estimated for geosynthetic sheets. In this latter case, the proportion is equal to 2 • tan (φ), where φ is the interface friction angle between the granular material and the sheet. Thus, assuming φ equal to 28.00 o , the ratio is equal to 1.06 for geosynthetic sheets. In the case of the geogrid, the ratio is greater (i.e., 2.81 to 5.49) because it integrates not only the friction between the granular material and the geogrid but also buttresses and shearing at the interface. The ratio between the division of pullout force by the surface area of the geogrid and the vertical stress at the geogrid interface for the pullout test samples. 6 Conclusion

Particles

Summary and reflection

A three-dimensional discrete modeling approach for geogrids has been provided to contribute to new studies related to the interface mechanisms between this type of geosynthetic reinforcement and granular materials. The use of DEM in this work was due to the current potential of this method for the analysis of problems involving interaction between elements with a relatively small number of parameters. Furthermore, exclusively discrete models do not have the difficulty of precisely defining the contact behavior between different numerical methods, such as numerical approaches involving, for example, the multi-domain coupling between DEM and FEM for particles and structures, respectively.

The geogrid model comprised deformable elements implemented in YADE, which is an auxiliary tool for DEM. The elements used were cylinders, characterized by two nodes and one connection, and particle facets (i.e., PFacets), characterized by three nodes, three connections, and two facets. The basis for the force-displacement law used in external and internal interactions of the model followed the classical linear elastic-plastic model of [START_REF] Cundall | A discrete numerical model for granular assemblies[END_REF]. Elastic limits for the internal forces defined the reinforcement failure condition.

The developed model has two possibilities for the numerical construction of the geogrid. The first possibility involved the discretization of the geogrid in segments with five nodes and four PFacets and the second in segments with four nodes and two PFacets. The first alternative has been shown to have a better symmetry condition, although the second one resulted in a lower computational cost due to the use of fewer discrete bodies for model construction. For the case under study, the two-PFacet configuration had similar behavior concerning the four-PFacet configuration. Therefore, it was the choice for the geogrid pullout simulations.

Geogrid traction and bending analyses made it possible to verify the numerical geogrid response. These analyses consisted of traction and deflection tests and subsequent comparison of results with analytical solutions. These analyzes also discussed the possible effect of reinforcement width reduction when under tensile conditions. It is possible to numerically suppress this effect by reducing the stiffness parameter of the diagonal connections that form the geogrid segments to zero or close to zero. Regarding the results obtained, the values computed numerically were consistent with those calculated analytically. In the case of the tensile test, the results overlapped and confirmed the expectation of matching the results. Turning to the deflection test, the relative deflection obtained by the numerical model was very close to the one obtained analytically from a parabolic curve for the deformation profile of the geogrid, which proved to be a good hypothesis. Concerning the calculated strain, it practically overlapped the simulation result, with the difference being negligible.

Tensile forces obtained in a specific section and different geogrid yarns exemplified that the numerical model allows analyzing particular parts of the structure. The inclusion of spherical particles in the geogrid model demonstrated the ability to take into account the interaction between the granular particles and the reinforcement. From the analyses, it was possible to conclude that the behavior of the numerical model was correct and consistent. The use of this numerical approach proved to be flexible and able to capture the continuous nature of the geogrid elements in discrete modeling. It also demonstrated to take into account attributes of actual geogrids, such as stretching, bending, shearing, and twisting, which are essential to properly describe their mechanical behavior.

The currently available geogrids can have different features inherent to the manufacturing method, including geometric particularities. A minimalist and very specific numerical representation restricts the application of the model and may also result in a high computational processing cost. The flat shape modeled for the geogrid members is representative and enables more realistic rolling mechanisms in interface simulations, which is an improvement considering current geogrid modeling referenced in the literature review of this thesis. Thus, the numerical approach focused on a generic representation of real geogrids involving attributes considered important for their interaction with granular materials.

Geogrid pullout simulations allowed to evaluate the mobilization of forces in the reinforcement during the tests. The analyzes comprised samples with different amounts of both spherical and non-spherical (i.e., clumps) particles. The results obtained were generally consistent with the expected behavior in regions of geosynthetic reinforcement. For example, the reinforcement region closest to the pullout force application concentrated the greatest axial stresses, as seen in the longitudinal elements of the two openings of the geogrid analyzed separately. The force-displacement curve of the geogrid pullout tests also showed an initial force peak followed by a tendency to a residual pullout force, especially when increasing the number of particles in the sample, which results in a smaller particle granular medium.

Interaction mechanisms between granular material and geogrid yarns indicated sensitivity to particle size. The peaks in force-displacement curves and stress-displacement curves are less pronounced in samples with smaller particle sizes. Despite less pronounced peaks, the geogrid transverse yarns showed curves with approximately the same range of values for samples with 3000, 6000, and 9000 particles, which means a smaller influence of both the number and geometry of the particles on the mechanisms of interaction with these yarns.

The simulations showed lower strength mobilization in the transverse members compared to the longitudinal members of the geogrid. For the geogrid openings individually analyzed in simulations comprising spherical particles, the average value of axial stresses in the transverse elements was 9 % of the average value in the longitudinal elements closest to the pullout region (i.e., cell 1) and 18 % of the average value in the longitudinal elements furthest from the pullout region (i.e., cell 2). In the case of clump-particle samples, these percentages were 8 % and 14 %, respectively. Although smaller than the stresses in the longitudinal elements, the axial stresses in the transverse elements of the geogrid are not negligible. The peaks of stress-displacement curves in these elements can exceed three times the calculated average value.

By comparing geogrid pullout tests between samples containing spherical and nonspherical particles (i.e., clumps), it was possible to analyze behaviors of the reinforced granular layer. The maximum pullout force (i.e., initial peak) occurs for a larger reinforcement displacement in the clump-particle samples. The values are also higher, and the pullout force is approximately twice as high in these samples for the presented cases.

Clump-particle samples characterized the first force peak of the pullout curve with sections of different slope, highlighting distinct mechanisms that occur at the interface between granular material and geogrid. These mechanisms are mainly due to the friction between granular material and geogrid yarns and the particle interlock that controls particle rolling at the interface with the reinforcement during the geogrid pullout.

The interlocking effect on non-spherical particles kept the same level of confinement at the interface of granular material and geogrid even after large geogrid displacement. The tendency to maintain pullout strength near the initial peak level rather than an abrupt drop to a residual value indicated this behavior.

Another aspect verified when comparing samples with spherical and non-spherical particles was that the ratio between the averages of calculated stress in corresponding elements of the geogrid resulted in approximately equal values for samples with the same number of elements. This behavior may indicate the possibility of using shape-matching factors in simulations that comprise spherical particles to mimic non-spherical particle behavior, which would reduce the computational cost of the model. Future studies devoted to this topic are needed to validate this hypothesis.

In the case of geogrid transverse yarns, the analyses indicated that the number of particles minimally influences the mechanism of interaction with the granular material. In contrast, the particle geometry of the granular medium is of great importance in the interaction with the transverse elements of the geogrid. The presence of non-spherical particles in the model caused axial stresses approximately twice as high in the transverse elements analyzed.

Regarding the properties and parameters that govern the behavior of the geogridreinforced layer numerical model, it is possible to separate the internal geogrid strength parameters and the geogrid and particle contact properties and parameters. Among the internal strength parameters, Young's modulus is one of the main, since it is related to the deformation behavior of the elements that make up the reinforcement. In the case of contact properties and parameters, the relative size between geogrid apertures and particles of granular material has a significant influence on the mobilization of reinforcement strength (i.e., especially in the longitudinal yarns of the geogrid).

In addition to the reinforcement mesh size, particle angularity also plays an essential role in the mechanical behavior of the geogrid reinforced sample. The geometrical characteristics of more angular or non-convex particles (i.e., clump-particle sample in this thesis) can restrict the rolling mechanism, which causes the interface between the geogrid and the granular material to have greater shear strength. Although the determination of contact parameters is 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.

Proposals for future work

The great advantage of a numerical approach is the numerous work possibilities it provides. Section 3.3 has presented preliminary models and simulations in which it is possible to know some of the many possibilities YADE offers to study the problem. For example, interface behavior analysis of geogrid-reinforced samples in simple shear tests. A triaxial cell model could also be used in geogrid pullout tests to ensure the initial confining conditions of the sample during the simulation.

One unexpected behavior in the pullout analyses of this thesis was the lower forcedisplacement curve values obtained for the 6000-clump sample compared to the other samples. As noted, the larger particle size of this sample resulted in a less uniform force chain network than in the 9000-clump sample. This less uniform force chain may explain the most pronounced force peaks in the force-displacement curves of the sample. However, it is unclear whether any particular arrangement of the particles in the geogrid interface region or another type of interaction mechanism caused the lowest average of the force-displacement curve found in the sample with 6000 clumps. Therefore, to better understand this occurrence, future studies could address clumps of different geometry, but with the same angularity (i.e., 0.79). A comparison between simulations involving the same amount of clumps may define whether or not this behavior is repeatable. In this case, the evaluation of contact forces in more members of the geogrid can help in understanding the mechanisms of interaction.

Another possibility for future work is to study the effect of geogrid multi-layers on the same sample. It would be very useful to understand the benefits of using multilayers solutions compared to just one geogrid layer of higher strength and stiffness. As mentioned in Section 2.1, the BS 8006-1 (BSI, 2010) highlights that this topic is not yet fully understood.

The yarn cross-overs modeling in the geogrid node region would be an addition or even an improvement to the numerical model proposed in this thesis. This modeling would allow representing more realistically two categories of geogrids that are the woven geogrids and bonded geogrids. Moreover, it would also make it possible to take into account node strength lower than that of yarns. Woven geogrids have the cross-overs joined by knitting or intertwining and a coating that protects the entire unit (e.g., bitumen, polyvinyl chloride, latex). Bonded geogrids comprise extruded strips of polyester or polypropylene welded together in a grid-like pattern. Modeling of cross-yarn nodes could be done using spheres that connect longitudinal and transverse yarns through a contact law that takes into account traction, compression, bending, and twisting (e.g., YADE C ++ class Law2 -ScGeom6D CohFrictPhys CohesionMoment). In the case of woven geogrids, it is possible to model yarns by juxtaposed cylindrical elements (i.e., side by side GridConnections) and, in bonded geogrids, by one of the geogrid segment configurations developed in this thesis (i.e., two-PFacet configuration and four-PFacet configuration).

On a final note, the use of numerical tools capable of dealing with the nature of granular materials and also with the complexity of material-reinforcement interaction can greatly assist the understanding of the behavior of the material-geogrid interface. In turn, this understanding is indispensable for an optimized design of reinforced structures. Therefore, the importance of developing numerical models that capture the attributes of real geogrids.

FIGURE 4 . 1 -

 41 FIGURE 4.1 -Important effects to be considered in geogrid modeling. . . . . . . .

FIGURE 4 . 2 -

 42 FIGURE 4.2 -Cylinder and particle facet elements. . . . . . . . . . . . . . . . . .

FIGURE 4 . 3 -

 43 FIGURE 4.3 -Numerical modeling of a generic geogrid in the DEM framework. . .

FIGURE 4 . 4 -

 44 FIGURE 4.4 -Contact scheme of two spheres. . . . . . . . . . . . . . . . . . . . .

FIGURE 4 . 5 -

 45 FIGURE 4.5 -Virtual sphere used in sphere-cylinder and cylinder-cylinder interactions. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

FIGURE 4 . 6 -

 46 FIGURE 4.6 -Depiction of the virtual sphere in a sphere-PFacet interaction. . . .

FIGURE 4 . 7 -

 47 FIGURE 4.7 -Types of induced reactions in members of the geogrid structure. . .

FIGURE 4 . 8 -

 48 FIGURE 4.8 -Geogrid segment configuration with four PFacets and two PFacets. .

FIGURE 4 . 9 -

 49 FIGURE 4.9 -Skeleton of the geogrid model. . . . . . . . . . . . . . . . . . . . . .

FIGURE 4 .

 4 FIGURE 4.10 -Length variation of the geogrid in the tensile direction and orthogonal to the tensile direction. . . . . . . . . . . . . . . . . . . . . . .

FIGURE 4 .

 4 FIGURE 4.11 -Rotation of transverse members in shear tensile simulation. . . . . .FIGURE 4.12 -Longitudinal and transverse strains of the chosen node of the geogrid.

FIGURE 4 .

 4 FIGURE 4.13 -Horizontal and vertical components of tensile force in the selected geogrid node. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

FIGURE 4 .

 4 FIGURE 4.14 -Moment-rotation curve of the selected geogrid node. . . . . . . . . .

FIGURE 4 .

 4 FIGURE 4.15 -Normal and shear tensile simulations for the geogrid. . . . . . . . .

FIGURE 4 .

 4 FIGURE 4.16 -Comparison of results for normal tensile simulation. . . . . . . . . .

FIGURE 4 .

 4 FIGURE 4.17 -Comparison of results for shear tensile simulation. . . . . . . . . . .

FIGURE 4 .

 4 FIGURE 4.18 -Initial and final state of the geogrid bending simulation. . . . . . . .

FIGURE 4 .

 4 FIGURE 4.19 -Comparison of relative deflection and strain for bending simulation.

FIGURE 4 .

 4 FIGURE 4.20 -Comparison of maximum tensile force for bending simulation. . . . .

FIGURE 4 .

 4 FIGURE 4.21 -Section S-S and yarns where the tensile forces were obtained. . . . .

FIGURE 4 .

 4 FIGURE 4.22 -Comparison of relative deflection and strain for bending simulation.

FIGURE 4 .

 4 FIGURE 4.23 -Views of the final state of the simulation involving interaction between geogrid and particle. . . . . . . . . . . . . . . . . . . . . . . .

FIGURE 4 .

 4 FIGURE 4.24 -Typical four-sphere clump. . . . . . . . . . . . . . . . . . . . . . . .

FIGURE 4 .

 4 FIGURE 4.25 -Angularity coefficient for two-sphere clump. . . . . . . . . . . . . . .

FIGURE 4 .

 4 FIGURE 4.26 -Geogrid-reinforced granular volumes during sample preparation. . .

FIGURE 4 .FIGURE 5 . 3 -FIGURE 5 . 4 -FIGURE 5 . 5 -

 4535455 FIGURE 4.27 -Samples prepared for geogrid pullout simulation. . . . . . . . . . . .

FIGURE 5 .

 5 FIGURE 5.9 -Geogrid pullout embedded in 3000 spheres (i.e., particle and mesh size ratio equal to 5.37E-01 m/m). . . . . . . . . . . . . . . . . . .

FIGURE

  FIGURE 5.10 -Geogrid pullout embedded in 6000 spheres (i.e., particle and mesh size ratio equal to 4.28E-01 m/m). . . . . . . . . . . . . . . . . . .

FIGURE 5 .

 5 FIGURE 5.35 -Front view of the pullout test on 6000-clump sample (i.e., particle and mesh size ratio equal to 4.28E-01 m/m). . . . . . . . . . . . . FIGURE 5.36 -Front view of the pullout test on 9000-clump sample (i.e., particle and mesh size ratio equal to 3.76E-01 m/m). . . . . . . . . . . . . FIGURE 5.37 -Comparison of the force-displacement curve between geogrid-reinforced samples composed of spheres and clumps. . . . . . . . . . . . . . . . FIGURE 5.38 -Comparison of the element A stress-displacement curve between geogrid-reinforced samples composed of spheres and clumps. . . . . FIGURE 5.39 -Comparison of the element B stress-displacement curve between geogrid-reinforced samples composed of spheres and clumps. . . . .

  FIGURE 1.1 -Ultimate limit states for basal reinforced embankments. Adapted from BSI (2010).

FIGURE 1

 1 FIGURE 1.2 -Serviceability limit states for basal reinforced embankments. Adapted from BSI (2010).

FIGURE 2

 2 FIGURE 2.1 -The sequence of a typical DEM modeling.

FIGURE 2

 2 FIGURE 2.2 -Two disks laterally compressed by rigid walls. Adapted from Cundall and Strack (1979).

Figure 2 .

 2 Figure 2.2 shows two negligible weight disks, represented by x and y, which are laterally compressed by a pair of rigid walls. The walls move toward each other at a constant velocity v. Initially, at time t0, the disks touch the walls with all contact forces equal to zero. After the time interval ∆t, the walls move at a distance equal to v • ∆t relative to the original position. Under the assumption that disturbance propagation can only occur through immediate neighboring elements, both disks remain in their positions until the time interval t1 (i.e., Figure 2.2).

FIGURE 2 . 3 -

 23 FIGURE 2.3 -Geogrid modeling from sphere aggregates (KONIETZKY et al., 2004).

  The method used by[START_REF] Chareyre | Dynamic spar elements and discrete element methods in two dimensions for the modeling of soil-inclusion problems[END_REF], Bourrier et al. (2013), and Effeindzourou et al. (2016) allows the modeling of deformable elements. Effeindzourou et al. (2016) demonstrated the versatility and capabilities of this approach and some potential applications to study soil-inclusion systems. This thesis highlights the cylinder presented by Bourrier et al. (2013) and the particle facets described by Effeindzourou et al. (2016). Effeindzourou et al. (2017) used these elements to model a damping system for rockfall protection. Damping systems can work as cushion layers for rockfall galleries. Effeindzourou et al. (2017) modeled a damping module composed of cylindrical wire mesh, two steel rings, a geotextile lining and granular filling material (i.e., Figure 2.4). The authors used cylinders to model the wire mesh and steel rings, particle facets to model the geotextile, and spheres to model the filling material.

  FIGURE 2.4 -Discrete modeling of damping module for rockfall protection (EFFEIND-ZOUROU et al., 2017).

FIGURE 2 . 5 -

 25 FIGURE 2.5 -Effect of top boundary condition on pullout test results (PALMERIA; MIL-LIGAN, 1989).

  FIGURE 2.6 -Influence of the initial geogrid position within the sample on the pullout test result (PALMEIRA, 1987; PALMEIRA, 2009).

FIGURE 2

 2 FIGURE 2.7 -Pullout test results for different grids. Adapted from Palmeira (2004).

  FIGURE 2.8 -The effect of confining pressure on the pullout response of a geogrid[START_REF] Farrag | Pull-out resistance of geogrid reinforcements[END_REF].

  Figure 3.1 summarizes a typical simulation cycle in YADE.

FIGURE 3 .

 3 FIGURE 3.1 -YADE's simulation cycle.

  (a) No collision due to lack of intersection between the set of points (b) No collision due to intersection on the x-axis only (c) Virtual collision by overlapping P i e P j (d) Exact collision detection between elements

FIGURE 3 .

 3 FIGURE 3.2 -SAP algorithm example involving two-dimensional collision detection scenarios. Adapted from Šmilauer et al. (2015a).

FIGURE 3 . 4 -

 34 FIGURE 3.4 -Representation of normal stiffness in the interaction between two spheres ( ŠMILAUER et al., 2015a).

  (a) 1 DOF for normal strain (b) 2 DOF for shearing (c) 1 DOF for twisting (d) 2 DOF for bending

FIGURE 3 . 5 -

 35 FIGURE 3.5 -Degrees of freedom for interaction between two spherical elements.

  position of the centers of the spheres at the current calculation time step; C o = interaction reference point at the current calculation time step; d 1 , d 2 = interaction radii (i.e., reduced or expanded radius of the sphere connecting its center to the point of contact); d 0 = the sum of interaction radii (i.e., d 0 = d 1 + d 2 ).

FIGURE 3 .

 3 FIGURE 3.6 -Contact created between two spheres, highlighting the overlapping between them. Adapted from Šmilauer et al. (2015a).

  14) where ∆u T 1 , ∆u T 2 = first part increments (i.e., interaction movement); u - T = shear displacement at the previous calculation time step; n -= the unit vector in the direction of the interaction axis at the previous calculation time step; ω 1 , ω 2 = respective angular velocities of elements with index 1 and 2 defined for the previous time step; ∆t = time step.

FIGURE 3 .

 3 FIGURE 3.7 -Shear component update due to displacement between elements. Adapted from Šmilauer et al. (2015a).

  acceleration at the current calculation time step; F = force acting on the element; m = mass of the element.

  3.38) where uo ξ = velocity components of the element at the current calculation time step; üo ξ = acceleration components of the element at the current calculation time step. YADE's engine NewtonIntegrator includes the implementation of Eq. 3.38 ( ŠMILAUER

  (a) Free-falling upper sphere and lower fixed sphere (b) Collision between the spheres (c) Evolution of the vertical displacement of the falling sphere

FIGURE 3 .

 3 FIGURE 3.8 -Collision between a free-falling sphere and a fixed sphere in the space.

  (a) Initial condition showing randomly arranged spheres (b) Final condition showing the spheres in equilibrium (c) Evolution of energy potentials, energy dissipations, and unbalanced force

FIGURE 3 .

 3 FIGURE 3.9 -Gravity deposition of spheres in a box.

  FIGURE 3.10 -Oedometer test with compression and decompression stages.

  (a) Shear applied towards the right side (b) Shear stress versus shear strain curve

FIGURE 3 .

 3 FIGURE 3.11 -Regular hexagonal packing of spheres under simple shear.

  (a) Deformed sample at the end of the test (b) Deviator stress and sample volume as a function of vertical strain

FIGURE 3 .

 3 FIGURE 3.12 -Triaxial test of a sample composed of clump particles.

FIGURE 3 .

 3 FIGURE 3.13 -Pullout test of a grid confined by spherical element layers.

Figure 3 .

 3 14 shows the four different geometries for the particles. The first clump comprises three aligned spheres with partial overlap (i.e., Figure 3.14b). Two aligned spheres and one out of alignment of different diameter in the middle of the other two assembled the second clump geometry (i.e., Figure 3.14c). The third geometry involves a cubic clump formed with four spheres (i.e., Figure 3.14d).

FIGURE 3 .

 3 FIGURE 3.14 -Particle geometry of the numerical sample.

  (a) Initial condition (b) Sample with initial porosity (c) Sample with specified porosity (d) Evolution of unbalanced force, energy potentials, energy dissipations, and porosity

FIGURE 3 .

 3 FIGURE 3.15 -Sample simulation of rigid agglomerates with specified porosity.

  Sphere-cylinder, sphere-facet, cylinder-cylinder, and cylinderfacet contacts (d) Sphere-cylinder, sphere-facet, cylinder-cylinder, cylinder-facet, and facet-facet contacts

FIGURE 3 .

 3 FIGURE 3.16 -Possible interactions involving spheres, cylinders, and particle facets.

FIGURE 4 .

 4 FIGURE 4.1 -Important effects to be considered in geogrid modeling.

FIGURE 4 .

 4 FIGURE 4.2 -Cylinder and particle facet elements.

FIGURE 4 . 3 -

 43 FIGURE 4.3 -Numerical modeling of a generic geogrid in the DEM framework.

FIGURE 4 . 4 -

 44 FIGURE 4.4 -Contact scheme of two spheres.

FIGURE 4 . 5 -

 45 FIGURE 4.5 -Virtual sphere used in sphere-cylinder and cylinder-cylinder interactions.

FIGURE 4 .

 4 FIGURE 4.6 -Depiction of the virtual sphere in a sphere-PFacet interaction.

M

  tw = twisting moment; k r = stiffness associated with the bending moment; k tw = stiffness associated with the twisting moment; Ω r = bending components of the relative rotations; Ω tw = twisting components of the relative rotations.

FIGURE 4 . 7 -

 47 FIGURE 4.7 -Types of induced reactions in members of the geogrid structure.

  internal normal force; F s = internal shear force; M r = internal bending moment; M tw = internal twisting moment; a n = normal adhesion; a s = shear adhesion; σ max n = tensile strength of the element; I r = bending moment of inertia (i.e., π • R 4 /4); I tw = polar moment of inertia (i.e., π • R 4 /2); R = internal reference radius (i.e., minimum radius between the interacting GridNodes).

  e = elastic shear displacement; Ω tw,e = the elastic part of the twisting component of the relative rotation; ν s = viscosity associated with shear; ν tw = viscosity associated with twisting; parameters YADE C ++ class yade.wrapper.

FIGURE 4 .

 4 FIGURE 4.8 -Geogrid segment configuration with four PFacets and two PFacets.

FIGURE 4 .

 4 FIGURE 4.9 -Skeleton of the geogrid model.

FIGURE 4 .

 4 FIGURE 4.10 -Length variation of the geogrid in the tensile direction and orthogonal to the tensile direction.

FIGURE 4 .

 4 FIGURE 4.11 -Rotation of transverse members in shear tensile simulation.

  Figure 4.15 shows the graphical output of the normal and shear tensile simulations for the geogrid. Both simulations were performed assuming absence of the width reduction effect. In YADE, tensile force is computed relative to the GridNodes. For the simulations, this vector quantity was obtained by summing the resultant forces on the nodes in which the traction was imposed (i.e., first row of GridNodes at the left end of the geogrid). The tensile force in the shear traction simulation is the resultant in the direction of the longitudinal elements, it is calculated by the sum of its horizontal and vertical components.Figures 4.16 and 4.17 present the comparison of tensile forces obtained by both Eq. 4.25 and the simulations. As can be seen, assuming the absence of the width reduction effect, the values overlap in this comparison, which confirms the expectation of matching the results.

FIGURE 4 .

 4 FIGURE 4.15 -Normal and shear tensile simulations for the geogrid.

  FIGURE 4.18 -Initial and final state of the geogrid bending simulation.

FIGURE 4 .

 4 FIGURE 4.20 -Comparison of maximum tensile force for bending simulation.

FIGURE 4 .

 4 FIGURE 4.21 -Section S-S and yarns where the tensile forces were obtained.

FIGURE 4 .

 4 FIGURE 4.22 -Comparison of relative deflection and strain for bending simulation.

Figure 4 .

 4 Figure 4.23 shows the final state of the simulation comprising the interaction between the geogrid and spherical particles. In comparison with the previous simulation, the maximum relative geosynthetic deflection was 0.07 m/m (i.e., 3.60E-02 m/m previously) and the tensile force per unit width in the same section S-S (i.e., 0.45 m from the left end of the geogrid) was 1.72 N/m (i.e., 102.90E-02 N/m previously). The magnitudes of the resultant forces on the spheres A and B (i.e., Figure 4.23) were respectively 0.22 N and 0.23 N (i.e., weight of the spheres).

FIGURE 4 .

 4 FIGURE 4.23 -Views of the final state of the simulation involving interaction between geogrid and particle.

FIGURE 4 .

 4 FIGURE 4.24 -Typical four-sphere clump.

FIGURE 4 .

 4 FIGURE 4.25 -Angularity coefficient for two-sphere clump.

  (a) Particles modeled as spheres (b) Particles modeled as clumps

FIGURE 4 .

 4 FIGURE 4.26 -Geogrid-reinforced granular volumes during sample preparation.

  (a) Sample using spheres (b) Sample using clumps

FIGURE 4 .

 4 FIGURE 4.27 -Samples prepared for geogrid pullout simulation.

FIGURE 5 . 1 -

 51 FIGURE 5.1 -Samples at the end of the triaxial test.

FIGURE 5 .

 5 FIGURE 5.7 -Deviator stress and volumetric strain versus axial strain for the 9000-clump sample.

FIGURE 5 .

 5 FIGURE 5.9 -Geogrid pullout embedded in 3000 spheres (i.e., particle and mesh size ratio equal to 5.37E-01 m/m).

FIGURE 5 .

 5 FIGURE 5.10 -Geogrid pullout embedded in 6000 spheres (i.e., particle and mesh size ratio equal to 4.28E-01 m/m).

FIGURE 5 .

 5 FIGURE 5.11 -Geogrid pullout embedded in 9000 spheres (i.e., particle and mesh size ratio equal to 3.76E-01 m/m).

  to 5.14 show the force chain network for the final state of the samples composed of spherical granular material.

FIGURE 5 .

 5 FIGURE 5.12 -Force chain network for the final state of the geogrid pullout test for the 3000-sphere sample.

FIGURE 5 .

 5 FIGURE 5.13 -Force chain network for the final state of the geogrid pullout test for the 6000-sphere sample.

FIGURE 5 .

 5 FIGURE 5.14 -Force chain network for the final state of the geogrid pullout test for the 9000-sphere sample.

FIGURE 5 .

 5 FIGURE 5.16 -Element A stress-displacement curve for geogrid pullout test on samples composed of spherical particles.

FIGURE 5 .

 5 FIGURE 5.20 -Front view of the pullout test on 3000-sphere sample (i.e., particle and mesh size ratio equal to 5.37E-01 m/m).

FIGURE 5 .

 5 FIGURE 5.21 -Front view of the pullout test on 6000-sphere sample (i.e., particle and mesh size ratio equal to 4.28E-01 m/m).

FIGURE 5 .

 5 FIGURE 5.22 -Front view of the pullout test on 9000-sphere sample (i.e., particle and mesh size ratio equal to 3.76E-01 m/m).

FIGURE 5 .

 5 FIGURE 5.23 -Geogrid pullout embedded in 3000 clumps (i.e., particle and mesh size ratio equal to 6.95E-01 m/m).

FIGURE 5 .

 5 FIGURE 5.24 -Geogrid pullout embedded in 6000 clumps (i.e., particle and mesh size ratio equal to 5.53E-01 m/m).

FIGURE 5 .

 5 FIGURE 5.25 -Geogrid pullout embedded in 9000 clumps (i.e., particle and mesh size ratio equal to 4.85E-01 m/m).

FIGURE 5 .

 5 FIGURE 5.27 -Force chain network for the final state of the geogrid pullout test for the 3000-clump sample.

FIGURE 5 .

 5 FIGURE 5.28 -Force chain network for the final state of the geogrid pullout test for the 6000-clump sample.

FIGURE 5 .

 5 FIGURE 5.29 -Force chain network for the final state of the geogrid pullout test for the 9000-clump sample.

  Figures 5.27 to 5.29, as well as Figures 5.23 to 5.25, show a lack of contact between the granular material and the upper plate at the right end of the sample. In contrast, there is a force chain with the upper plate on the left side of the box (i.e., imposing side of the geogrid pullout). Therefore, vertical stresses at the geogrid level may vary over the tests. For the final state of samples with 3000, 6000, and 9000 clumps (i.e.,Figures 5.23 to 5.25), the total normal force applied to the top plate of the box is 5.85E+03 N , 5.19E+03 N , and 5.92E+03 N , respectively.

FIGURE 5 .

 5 FIGURE 5.34 -Front view of the pullout test on 3000-clump sample (i.e., particle and mesh size ratio equal to 5.37E-01 m/m).

FIGURE 5 .

 5 FIGURE 5.35 -Front view of the pullout test on 6000-clump sample (i.e., particle and mesh size ratio equal to 4.28E-01 m/m).

FIGURE 5 .

 5 FIGURE 5.36 -Front view of the pullout test on 9000-clump sample (i.e., particle and mesh size ratio equal to 3.76E-01 m/m).

  and 5.37c indicate a third slope change in approximately 0.06 m of geogrid displacement.

FIGURE 5 .

 5 FIGURE 5.37 -Comparison of the force-displacement curve between geogrid-reinforced samples composed of spheres and clumps.

  Géogrille, Méthode des Eléments Discrets, Renforcement des Sols, Essai d'Extraction List of Figures FIGURE 1.1 -Ultimate limit states for basal reinforced embankments. Adapted from BSI (2010). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Pullout test results for different grids. Adapted from Palmeira (2004). FIGURE 2.8 -The effect of confining pressure on the pullout response of a geogrid (FARRAG et al., 1993). . . . . . . . . . . . . . . . . . . . . . . . . . . FIGURE 3.1 -YADE's simulation cycle. . . . . . . . . . . . . . . . . . . . . . . . . FIGURE 3.2 -SAP algorithm example involving two-dimensional collision detection scenarios. Adapted from Šmilauer et al. (2015a). . . . . . . . . 3.5 -Degrees of freedom for interaction between two spherical elements. . FIGURE 3.6 -Contact created between two spheres, highlighting the overlapping between them. Adapted from Šmilauer et al. (2015a). . . . . . . . . FIGURE 3.7 -Shear component update due to displacement between elements. Adapted from Šmilauer et al. (2015a). . . . . . . . . . . . . . . . . . FIGURE 3.8 -Collision between a free-falling sphere and a fixed sphere in the space. FIGURE 3.9 -Gravity deposition of spheres in a box. . . . . . . . . . . . . . . . . FIGURE 3.10 -Oedometer test with compression and decompression stages. . . . . FIGURE 3.11 -Regular hexagonal packing of spheres under simple shear. . . . . . .

FIGURE 1.2 -Serviceability limit states for basal reinforced embankments. Adapted from BSI (2010). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . FIGURE 2.1 -The sequence of a typical DEM modeling. . . . . . . . . . . . . . . . FIGURE 2.2 -Two disks laterally compressed by rigid walls. Adapted from Cundall and Strack (1979). . . . . . . . . . . . . . . . . . . . . . . . . . FIGURE 2.3 -Geogrid modeling from sphere aggregates (KONIETZKY et al., 2004). FIGURE 2.4 -Discrete modeling of damping module for rockfall protection (EF-FEINDZOUROU et al., 2017). . . . . . . . . . . . . . . . . . . . . . . . FIGURE 2.5 -Effect of top boundary condition on pullout test results (PALMERIA; MILLIGAN, 1989). . . . . . . . . . . . . . . . . . . . . . . . . . . . . FIGURE 2.6 -Influence of the initial geogrid position within the sample on the pullout test result (PALMEIRA, 1987; PALMEIRA, 2009). . . . . . . . FIGURE 2.7 -FIGURE 3.3 -Simplified representation of the contact between two spheres. Adapted from Coetzee and Els (2009). . . . . . . . . . . . . . . . . . . . . . . FIGURE 3.4 -Representation of normal stiffness in the interaction between two spheres ( ŠMILAUER et al., 2015a). . . . . . . . . . . . . . . . . . . . . FIGURE FIGURE 3.12 -Triaxial test of a sample composed of clump particles. . . . . . . . . FIGURE 3.13 -Pullout test of a grid confined by spherical element layers. . . . . . . FIGURE 3.14 -Particle geometry of the numerical sample. . . . . . . . . . . . . . . FIGURE 3.15 -Sample simulation of rigid agglomerates with specified porosity. . . . FIGURE 3.16 -Possible interactions involving spheres, cylinders, and particle facets.

  5.10 -Geogrid pullout embedded in 6000 spheres (i.e., particle and mesh size ratio equal to 4.28E-01 m/m). . . . . . . . . . . . . . . . . . . embedded in 9000 spheres (i.e., particle and mesh size ratio equal to 3.76E-01 m/m). . . . . . . . . . . . . . . . . . . FIGURE 5.12 -Force chain network for the final state of the geogrid pullout test for the 3000-sphere sample. . . . . . . . . . . . . . . . . . . . . . . . . . FIGURE 5.13 -Force chain network for the final state of the geogrid pullout test for the 6000-sphere sample. . . . . . . . . . . . . . . . . . . . . . . . . . FIGURE 5.14 -Force chain network for the final state of the geogrid pullout test for the 9000-sphere sample. . . . . . . . . . . . . . . . . . . . . . . . . . 5.22 -Front view of the pullout test on 9000-sphere sample (i.e., particle and mesh size ratio equal to 3.76E-01 m/m). . . . . . . . . . . . . FIGURE 5.23 -Geogrid pullout embedded in 3000 clumps (i.e., particle and mesh size ratio equal to 6.95E-01 m/m). . . . . . . . . . . . . . . . . . . FIGURE 5.24 -Geogrid pullout embedded in 6000 clumps (i.e., particle and mesh size ratio equal to 5.53E-01 m/m). . . . . . . . . . . . . . . . . . . 5.27 -Force chain network for the final state of the geogrid pullout test for the 3000-clump sample. . . . . . . . . . . . . . . . . . . . . . . . . . FIGURE 5.28 -Force chain network for the final state of the geogrid pullout test for the 6000-clump sample. . . . . . . . . . . . . . . . . . . . . . . . . . FIGURE 5.29 -Force chain network for the final state of the geogrid pullout test for the 9000-clump sample. . . . . . . . . . . . . . . . . . . . . . . . . . FIGURE 5.30 -Element A stress-displacement curve for geogrid pullout test on samples composed of clump particles. . . . . . . . . . . . . . . . . . . . FIGURE 5.31 -Element B stress-displacement curve for geogrid pullout test on samples composed of clump particles. . . . . . . . . . . . . . . . . . . .

curve for geogrid pullout test on samples composed of spherical particles. . . . . . . . . . . . . . . . . . . . . . . . FIGURE 5.16 -Element A stress-displacement curve for geogrid pullout test on samples composed of spherical particles. . . . . . . . . . . . . . . . . . . FIGURE 5.17 -Element B stress-displacement curve for geogrid pullout test on samples composed of spherical particles. . . . . . . . . . . . . . . . . . . FIGURE 5.18 -Element C Stress-displacement curve for geogrid pullout test on samples composed of spherical particles. . . . . . . . . . . . . . . . . . . FIGURE 5.19 -Element D stress-displacement curve for geogrid pullout test on samples composed of spherical particles. . . . . . . . . . . . . . . . . . . FIGURE 5.20 -Front view of the pullout test on 3000-sphere sample (i.e., particle and mesh size ratio equal to 5.37E-01 m/m). . . . . . . . . . . . . FIGURE 5.21 -Front view of the pullout test on 6000-sphere sample (i.e., particle and mesh size ratio equal to 4.28E-01 m/m). . . . . . . . . . . . . FIGURE FIGURE 5.25 -Geogrid pullout embedded in 9000 clumps (i.e., particle and mesh size ratio equal to 4.85E-01 m/m). . . . . . . . . . . . . . . . . . . FIGURE 5.26 -Force-displacement curve for geogrid pullout testin on samples composed of clump particles. . . . . . . . . . . . . . . . . . . . . . . . . FIGURE FIGURE 5.32 -Element C Stress-displacement curve for geogrid pullout test on samples composed of clump particles. . . . . . . . . . . . . . . . . . . . FIGURE 5.33 -Element D stress-displacement curve for geogrid pullout test on samples composed of clump particles. . . . . . . . . . . . . . . . . . . . FIGURE 5.34 -Front view of the pullout test on 3000-clump sample (i.e., particle and mesh size ratio equal to 5.37E-01 m/m). . . . . . . . . . . . .

TABLE 4 .

 4 1 -Simulation parameters. . . . . . . . . . . . . . . . . . . . . . . . . .

TABLE 4 .

 4 2 -Angularity coefficient for the typical 4-sphere clump. . . . . . . . . .

TABLE 5 .

 5

1 -Porosity and particle size of samples. . . . . . . . . . . . . . . . . .

TABLE 5.2 -Geogrid pullout test model parameters. . . . . . . . . . . . . . . . . TABLE 5.3 -Macroscopic parameters of granular material samples. . . . . . . . . TABLE 5.4 -Mean values of the stress-displacement curve in elements A and B for samples composed of spherical particles. . . . . . . . . . . . . . . TABLE 5.5 -Mean values of the stress-displacement curve in elements A and B for samples composed of clump particles. . . . . . . . . . . . . . . . TABLE 5.6 -Ratio of cell 1 to cell 2 for mean values of stress-displacement curves in elements A and B. . . . . . . . . . . . . . . . . . . . . . . . . . . TABLE

TABLE 4 .

 4 1 -Simulation parameters.

	Category	Description	Value	Unit
		Length (L o )	1.00	m
		Width	0.50	m
		Thickness	4.00E-03	m
	Contact	Aperture area	6.85E-03	m 2
		Yarn cross-section	2.51E-05	m 2
		Density	1.00E+02 kg/m 3
		Young's modulus (E c )	2.00E+06 N/m 2
	Internal	Poisson's ratio (ν n

  Longitudinal and transverse strains of the chosen node of the geogrid.

		0.00e+00		1.00e+04	2.00e+04	3.00e+04	4.00e+04	5.00e+04
			Longitudinal strain
			Transverse strain
		0.00e+00		1.00e+04	2.00e+04	3.00e+04	4.00e+04	5.00e+04
					Interaction number
	-5.00 -2.00 -1.00 0.00 1.00 2.00 0.00e+00 0.00e+00 FIGURE 4.12 -Tensile resultant force (N) -4.00 -3.00	1.00e+04 1.00e+04 x-axis component 2.00e+04 2.00e+04 y-axis component	3.00e+04 3.00e+04	4.00e+04 4.00e+04	5.00e+04 -5.00 -2.00 -1.00 0.00 1.00 2.00 5.00e+04 -3.00 -4.00
					Interaction number
	FIGURE 4.13 -Horizontal and vertical components of tensile force in the selected geogrid
	node.
	Bending moment (N.m)	-2.00e-05 0.00e+00 2.00e-05 4.00e-05 6.00e-05 0.00 8.00e-05 0.00	0.10 0.10	0.20 0.20	0.30 0.30	0.40 0.40	0.50 -2.00e-05 0.00e+00 2.00e-05 4.00e-05 6.00e-05 8.00e-05 0.50
					Rotation (rad)

FIGURE 4.

14 -Moment-rotation curve of the selected geogrid node.

TABLE 5 .

 5 1 -Porosity and particle size of samples.

	Particles	Description	3000	6000	9000	Unit
		Initial porosity	5.15E-01 5.02E-01 4.96E-01	-
		Minimum porosity	4.54E-01 4.44E-01 4.38E-01	-
		Minimum particle size	2.20E-02 1.76E-02 1.54E-02	m
	Spheres
		Maximum particle size	6.60E-02 5.27E-02 4.62E-02	m
		Average particle size	4.40E-02 3.51E-02 3.08E-02	m
		Particle and mesh size ratio 5.37E-01 4.28E-01 3.76E-01 m/m
		Initial porosity	4.97E-01 4.99E-01 4.99E-01	-
		Minimum porosity	3.58E-01 3.49E-01 3.45E-01	-
		Minimum particle size	3.42E-02 2.72E-02 2.38E-02	m
	Clumps
		Maximum particle size	7.97E-02 6.36E-02 5.56E-02	m
		Average particle size	5.70E-02 4.54E-02 3.98E-02	m

TABLE 5 .

 5 2 -Geogrid pullout test model parameters.

	Model	Parameter	Value	Unit
		Length	1.00	m
	Box	Width	0.50	m
		Height	0.50	m
		Length	1.00	m
		Width	0.50	m
		Thickness	4.00E-03	m
		Mesh size	0.08 × 0.08	m
		Yarn cross-section	2.51E-05	m 2
	True Biaxial	Density	1.00E+02 kg/m 3
	Geogrid	Contact friction angle	20.00	o
		Young's modulus	4.00E+09 N/m 2
		Poisson's ratio (K N /K T )	0.30	-
		Normal cohesion	∞	N/m 2
		Shear cohesion	∞	N/m 2
		Density	2.60E+03 kg/m 3
	Granular Material	Contact friction angle	30.00	o
	(microscopic parameters)	Young's modulus	1.00E+06 N/m 2
		Poisson's ratio (K N /K T )	0.50	-

TABLE 5 .

 5

		3 -Macroscopic parameters of granular material samples.
	Particles	Description	3000	6000	9000	Unit
		Density	1.45E+03 1.47E+03 1.48E+03 kg/m 3
		Friction angle	31.00	31.66	31.98	o
	Spheres
		Young's modulus 9.51E+04 8.53E+04 8.45E+04 N/m 2
		Poisson's ratio	0.44	0.40	0.38	-
		Density	1.70E+03 1.71E+03 1.72E+03 kg/m 3
		Friction angle	50.48	50.43	50.07
	Clumps

o Young's modulus 1.51E+05 1.70E+05 1.60E+05 N/m 2

TABLE 5 .

 5 4 -Mean values of the stress-displacement curve in elements A and B for samples composed of spherical particles. Element C Stress-displacement curve for geogrid pullout test on samples composed of spherical particles.

	Element	Number of spheres Average value Unit
		3000	2.07E+07	N/m 2
	Element A -cell 1	6000	2.51E+07	N/m 2
		9000	2.31E+07	N/m 2
		3000	1.80E+07	N/m 2
	Element B -cell 1	6000	2.34E+07	N/m 2
		9000	2.30E+07	N/m 2
		3000	1.18E+07	N/m 2
	Element A -cell 2	6000	1.24E+07	N/m 2
		9000	1.21E+07	N/m 2
		3000	1.03E+07	N/m 2
	Element B -cell 2	6000	1.17E+07	N/m 2
		9000	1.10E+07	N/m 2
	less, the stresses mobilized in the transverse elements are not negligible, their average value
	(i.e., 2.06E+06 N/m 2 ) is approximately 9 % and 18 % of the average value obtained for
	the longitudinal elements of cells 1 (i.e., 2.22E+07 N/m 2 ) and 2 (i.e., 1.16E+07 N/m 2 ),
	respectively. Also, it is important to note that stress peaks in the transverse elements
	during geogrid pullout can exceed four times the calculated average value.

Element B stress (N/m²)

	1.00e+08 0.00	0.10	0.20	0.30	0.40	1.00e+08 0.50
	8.00e+07					8.00e+07
	6.00e+07					6.00e+07
	4.00e+07					4.00e+07
	2.00e+07					2.00e+07
	0.00e+00 0.00	0.10	0.20	0.30	0.40	0.00e+00 0.50
			Displacement (m)
	3000 clumps (cell1)	6000 clumps (cell1)	9000 clumps (cell1)
	3000 clumps (cell2)	6000 clumps (cell2)	9000 clumps (cell2)
	FIGURE 5.31 -Element B stress-displacement curve for geogrid pullout test on samples
	composed of clump particles.

TABLE 5 .

 5 5 -Mean values of the stress-displacement curve in elements A and B for samples composed of clump particles.

	Element	Number of clumps Average value Unit
		3000	5.62E+07	N/m 2
	Element A -cell 1	6000	4.71E+07	N/m 2
		9000	5.28E+07	N/m 2
		3000	5.59E+07	N/m 2
	Element B -cell 1	6000	4.39E+07	N/m 2
		9000	5.37E+07	N/m 2
		3000	2.59E+07	N/m 2
	Element A -cell 2	6000	2.45E+07	N/m 2
		9000	2.94E+07	N/m 2
		3000	3.08E+07	N/m 2
	Element B -cell 2	6000	2.29E+07	N/m 2
		9000	2.94E+07	N/m 2

Element C stress (N/m²)

	3.00e+07 0.00	0.10	0.20	0.30	0.40	3.00e+07 0.50
		MEAN VALUE
	2.00e+07					2.00e+07
	1.00e+07					1.00e+07
	0.00e+00 0.00	0.10	0.20	0.30	0.40	0.00e+00 0.50
			Displacement (m)
	3000 clumps (cell1)	6000 clumps (cell1)	9000 clumps (cell1)
	3000 clumps (cell2)	6000 clumps (cell2)	9000 clumps (cell2)
	FIGURE 5.33 -Element D stress-displacement curve for geogrid pullout test on samples
	composed of clump particles.
	3.00e+07 0.00	0.10	0.20	0.30	0.40	3.00e+07 0.50
		MEAN VALUE
	2.00e+07					2.00e+07
	1.00e+07					1.00e+07
	0.00e+00 0.00	0.10	0.20	0.30	0.40	0.00e+00 0.50
			Displacement (m)
	3000 clumps (cell1)	6000 clumps (cell1)	9000 clumps (cell1)
	3000 clumps (cell2)	6000 clumps (cell2)	9000 clumps (cell2)
	FIGURE 5.32 -Element C Stress-displacement curve for geogrid pullout test on samples
	composed of clump particles.

Element D stress (N/m²)

TABLE 5 .

 5 6 -Ratio of cell 1 to cell 2 for mean values of stress-displacement curves in elements A and B.

	Particle type Element Number of particles Cell 1 to cell 2 ratio
			3000	1.75
	Sphere	A	6000	2.02
			9000	1.91
			3000	1.74
	Sphere	B	6000	2.00
			9000	2.10
			3000	2.17
	Clump	A	6000	1.92
			9000	1.80
			3000	1.82
	Clump	B	6000	1.92
			9000	1.82

TABLE 5 .

 5 7 -Ratio of clump-particle sample to sphere-particle sample for mean values of stress-displacement curves in corresponding elements.

	Element A stress (N/m²)	0.00e+00 2.00e+07 4.00e+07 6.00e+07 8.00e+07 1.00e+08

(c) Element A stress-displacement curve for geogrid pullout test on 9000-particle samples FIGURE 5.38 -Comparison of the element A stress-displacement curve between geogridreinforced samples composed of spheres and clumps.

(c) Element B stress-displacement curve for geogrid pullout test on 9000-particle samples

FIGURE 5.39 -Comparison of the element B stress-displacement curve between geogridreinforced samples composed of spheres and clumps. simpler geometry particles.

The open-source code is available at: <https://yade-dem.org/>

Acknowledgments

The Convention de Cotutelle de Thèse (joint PhD agreement) between Université Grenoble Alpes (ComUE -UGA -3SR) -France and Instituto Tecnológico de Aeronáutica (ITA) -Brazil gave rise to the work presented in this thesis. The integrated experience and access to the two research environments for the benefit of this thesis was a great privilege.

This study was financed in part by the Coordenação de Aperfeiçoamento de Pessoal de Nível Superior -Brasil (CAPES) -Finance Code 001. Laboratoire 3SR (CNRS -Grenoble INP -UGA) also supported scientific training courses and events. I am very grateful for the financial support received that made this work possible.

Many people have helped me and made the path possible for my studies and research. Their contributions are sincerely appreciated and gratefully acknowledged. In particular, I would like to express my appreciation to:

Wilson Cabral for the support provided by the Engenharia de Infraestrutura Aeronáutica Graduate Program -ITA.

Jean-Pierre Gourc for introducing me to Laboratoire 3SR and making the Convention de Cotutelle de Thèse possible.

Delma Vidal for paving the way for the French research environment through her professionally developed partnerships, her kindness to me, and the hours she spent discussing my work.