Simulation of the Mechanical Behavior of Fiber Reinforced Sand using the Discrete Element Method

. The general characteristics of granular soils reinforced with fibres have been reported in previous studies and have shown that fibre inclusion provides an increase in material strength and ductility and that the composite behaviour is governed by fibre content, as well as the mechanical and geometrical properties of the fibre. The present work presents a numerical procedure to incorporate fiber elements into an existing discrete element code (GeoDEM). The fiber elements are represented by linear elastic-plastic segments that connect two neighbor contacts where the fiber is located. These elements are characterized by an axial stiffness, tensile strength and length. The effect of the addition of fibers was evaluated numerically by comparing the stress-strain behavior of a pure sand with and without fibers. These simulations showed that the addition of fibers provides a significant increase in strength for the mixture in comparison with strength of the pure sand.


Introduction
Past research has demonstrated that inclusion of fibers significantly improves the engineering response of soils. Gray & Ohashi (1983) studied the mechanics of fiber-reinforcement in cohesionless soils and showed that inclusion of fibers increased peak shear strength and ductility of soils under static loads. A number of factors such as fiber content, orientation of fibers with respect to the shear surface, and the elastic modulus of the fiber were each found to influence the contribution of fibers to the shear strength. Later work (e.g. Maher & Gray, 1990;Maher & Ho, 1993;Morel & Gourc, 1997;Consoli et al., 1998Consoli et al., , 1999Consoli et al., , 2002Consoli et al., , 2004Consoli et al., , 2005Consoli et al., , 2007aConsoli et al., , 2007bConsoli et al., , 2009Zornberg, 2002;Vendruscolo, 2003;Michalowski & Cermák, 2003;Heineck et al., 2005;Casagrande, 2005 andCasagrande et al., 2006) has improved understanding of the mechanisms involved and the parameters affecting the behavior of fiber reinforced soils under static loading conditions. Within this context, the Discrete Element Method (DEM) was used in the present work, as a modeling tool for analysis of micro-scale reinforced soil, to understand how these materials interact. In particular evaluation was made of the behavior of soil-polypropylene fibers mixtures, when compared to a non-reinforced soil.
DEM was introduced by Cundall (1971) for the analysis of rock mechanics problems and later specialized to granular media by Cundall & Strack (1979). Since then, the DEM has been utilized to represent the mechanical behavior of pure sands subjected to different loading conditions Plassiard et al., 2009 andZhao &Evans, 2009). However, little work has been reported on the behavior of sand mixed with fibers using DEM. One of these studies, Maeda & Ibraim (2008) used DEM for the study of granular soils reinforced with fibers. In the approach adopted by these authors, fibers were modeled as particles themselves connected by contact bonds possessing tensile and shear strength but no restriction to rotation, behaving as hinge connections.
A two dimensional DEM simulation of the mixture of granular soil and fibers was carried out by implementing a special procedure to take into account the fiber-soil interaction. This implementation was carried out in the numerical code GeoDEM (Velloso, 2010) which is a general purpose discrete element code. Numerical results under biaxial compression are presented and discussed.
DEM is a method that carries relatively high computational costs. Simulation times for a fair sized problem may take hours or even days to process. However, the proposed formulation for inclusion of fibers inside the granular mass does not alter in a significant way the demands in computing power.
The main objective of the present work is to describe and to evaluate the proposed procedure for the inclusion of fibers into the DEM.

The discrete element method
Motion of particles and their interaction are modeled in the present work by the Discrete Element Method (DEM) implemented in the GeoDEM code (Velloso, 2010). DEM was introduced by Cundall (1971) for the analysis of rock mechanics problems and later specialized to granular media by Cundall & Strack (1979). In this method, Newton's second law is used to describe the motion of an individual particle. The equation that governs the translation motion of an individual particle i is: where m i e v i are respectively mass and velocity of particle i, nc i is the number of contacts of particle i. Forces acting on particle are gravitational force (F gi ), contact forces between particles i and j (F cij ), forces due to fibers (F fi ). The equation that governs the rotational motion of particle i is: where w i is the angular velocity of particle i, I i is the moment of inertia of particle i, and T ij is the torque generated by the contact forces between particles i and j. Motion Eqs. 1 and 2 are integrated using a central difference scheme. Details of expressions for accelerations, velocities and displacements can be found in Cundall & Strack (1979). The mechanical (macroscopic) behavior of a particular material is simulated in DEM through constitutive models associated to each contact (through the force-displacement relationship). By using the adopted constitutive model, the contact forces between two particles can be determined. In the present work, a linear elastic-perfectly plastic (sliding) behavior, whose graphical representation is shown in Fig. 1, was adopted.
The contact force F c can be decomposed in normal F c N and shear forces F c S : Normal force is given by: where K N is contact normal stiffness, n is the normal vector to the contact plane and U N is the superposition between particles.
being R 1 e R 2 particle radii in contact and d the distance between particles centroid.
The shear force is determined in an incremental form. When a contact is formed, the shear force is initially zero and the subsequent shear displacements result in increments of this force: being v c s the relative velocity in the shear direction at the contact, k s the contact shear stiffness and Dt the time interval. In DEM application to sands, tensile forces at the contact are not allowed and therefore: The sliding model between particles is defined by the friction coefficient (m) that limits the shear force to a maximum value given by: If F c S c s F ³ max , sliding occurs making: After calculations of the contact forces, these are transferred to the particles and used in Eqs. 1 and 2. Equations 1 and 2 correspond to dynamic conditions. Natural dynamic systems contain some degree of damping otherwise the system would oscillate indefinitely. Various types of damping have been proposed in the solution of problems using the DEM. For the quasi-static problems, such as the ones included in the present work, local, non-viscous damping is indicated. Details of characteristics of different forms of damping can be found in Figueiredo (1991) and in the user's manual of the code PFC2D (Itasca, 2002).

The fiber element
In this work, it is proposed the addition of fiber elements in the code GeoDEM (Velloso, 2010). The implementation is carried out by the addition of fiber elements in GeoDEM through the use of linear elastic-plastic segments connecting two neighboring contacts as shown in Fig. 2 202 Soils and Rocks, São Paulo, 35 (2) (Ferreira, 2010). These fiber elements are characterized by axial stiffness, tensile strength and fiber length. The relative displacement of these contacts generates a force that is transmitted to each pair of particles that forms the contact. The force due to the fiber is calculated as: where F f is the fiber related force to be transmitted to the particles, K f is the fiber stiffness, L f is the present length of the fiber segment, L f0 is the initial length (at the moment the fiber was generated) of the fiber segment and n f is the unit vector in the direction that connects the two neighboring contacts. If the fiber segment shortens, (L f < L f0 ), the fiber force is zero, that is, the fiber only transfers forces to the particles when subjected to tensile forces. If the force calculated by Eq. 10 in a fiber segment is larger than the tensile strength of the fiber, the fiber force becomes zero, this fiber segment is disabled, indicating failure of the fiber at this point. With the addition of the fiber elements, the calculation cycle of DEM is shown in Fig. 3. The generation of fibers is carried out after the granular medium reaches an equilibrium configuration. Fibers can be generated stochastically, whose distribution are described by average and standard deviation values for length and direction. Fibers can be generated in the whole domain or in specific and prescribe regions.

Results
A granular sample 5 mm wide and 10 mm high formed by circular particles having an average diameter of 0.25 mm and a 20% porosity was generated and is shown in Fig. 4. The micromechanical parameters, shown in Table 1, were calibrated in a way to obtain the stress strain behavior of a sand having an effective friction angle of 30°, and Young's modulus of approximately 10 MPa for a confining stress of 100 kPa. One hundred fiber elements having a length equal to five times the average diameter of the particles were generated in random directions of the sample.   The axial stiffness of the fibers is 100 MN/m and their tensile strength is 10 kN. Figure 5 shows the stress-strain curves obtained from the numerical simulations of a sample of pure sand and one of pure sand with the inclusion of fibers. One is able to observe a significant increase in strength with the addition of fibers. This increase in strength is also reflected in in the strength envelope represented in Fig. 6. Figure 7 shows the influence of the fiber stiffness in the mechanical behavior of the tested samples. One is able to note that an increase in stiffness of the fibers, as implemented in the present work, has the effect of increasing both strength and stiffness of samples (reflected by their Young's modulus).
Experimental work (Casagrande, 2005) has shown that the addition of fibers increase both strength and ductility of sandy soils. The increase in ductility was not observed in the numerical simulations carried out. A possible reason for this discrepancy is related to the elasto-plastic model adopted for the fibers. This model incorporates soft-ening, that is, fibers fail completely when reaching their tensile strength. 204 Soils and Rocks, São Paulo, 35 (2)

Conclusions
This work presented a simple procedure for the inclusion of fiber elements in the discrete element formulation. In the proposed procedure, the fiber element is represented by a virtual segment connecting particles.
Numerical simulations of biaxial compression tests were carried out in samples of sand and sand mixed with fibers. These simulations showed that the addition of fibers provides a significant increase in strength for the mixture in comparison with strength of the pure sand.
However, the increase in ductility due to the addition of fibers as reported in the literature was not fully observed in the numerical simulations, indicating that the constitutive model adopted for the fibers has to be improved in order to reproduce all features observed in the experimental work A 2D simulation has some limitations to the system, such as simulating the actual interaction between the granular material and the fibers and their continuity. All phenomena implemented for the 2D case can be extended to 3D case; however, the 2D model has shown, in a satisfactory case, the increase in soil strength with fibers insertion.