Discrete Element Analysis of the Strength Anisotropy of Fiber-Reinforced Sands Subjected to Direct Shear Load

Abstract

Soil reinforcement with natural or synthetic fibers enhances its mechanical behavior in various applications. Fiber-reinforced sands (FRS) can be relatively anisotropic because of the fiber self-weight and the compaction technique. However, the microscopic mechanisms underlying the anisotropy are still poorly understood. This study used a discrete element approach to analyze the microscopic mechanisms underlying the strength anisotropy of FRS due to fiber orientation. Analysis of contact networks revealed that the optimum fiber orientation angle is perpendicular to the main direction of strong contact force in direct shear testing. These fibers produced the largest increase in shear zone thickness, normal force around the fiber body, effective contact area, tensile force along fibers, and energy storage/dissipation. This study is valuable for further understanding of the mechanical behaviors of FRS.

1. Introduction

In geotechnical engineering, various soil reinforcement techniques have been commonly used to enhance the engineering properties of soil. Soil reinforcement by fiber materials is considered an effective ground improvement method because of its cost-effectiveness, easy adaptability, and reproducibility [1]. Of particular concern is the strength of fiber-reinforced soils (FRS), which varies according to different fiber and soil properties. This indicates a need to understand the reinforcement mechanisms responsible for these differences by exploration with different perspectives and methods.

Factors that influence the shear strength of fiber-reinforced soil have been explored in a variety of laboratory experiments using, for example, triaxial tests and direct shear tests [2,3,4,5,6,7,8,9,10]. The existing body of research suggests that different soil conditions (i.e., the density, particle shape, and size grading) and fiber traits (i.e., type, fraction, dimension, and orientation) present significant effects on the reinforcement of fibers [11]. Several theories on the reinforcing mechanisms have also been proposed, such as force-equilibrium/mechanistic, energy dissipation, statistical, and superposition of the effects of soil and fibers [11,12,13,14,15,16]. In addition, different authors have investigated different problems in fiber reinforced soil using numerical methods, including continuum models [17] and discrete element method (DEM) models [18,19,20]. In a comparison of the two numerical approaches, Mao et al. [21] quantified and evaluated the effect of root reinforcement by simulating direct shear tests.

The principal advantages of fiber-reinforced soil are generally considered to be the relative isotropy of the composite material, and many researchers have assumed isotropy for simplicity during analysis [22]. However, Michalowski [23] and Diambra et al. [24] suggested that anisotropy would be a more realistic assumption. They argue that anisotropy is brought about by the limitations of the construction equipment and the method by which the fibers are incorporated into the soil and then compacted. The design shear strength of the fiber-soil composites should account for the preferential orientation of fibers. So far, however, there has been little discussion about the influence of fiber orientation on reinforcement [11,23,25,26,27,28,29]. In addition, these previous studies can only obtain macroscopic responses and far too little attention has been paid to the micromechanical mechanism that governs the reinforcement. These two aspects motivated this study. This article places emphasis on the analysis of the cause of the shear strength anisotropy of fiber reinforced sands (FRS) induced by a preferential orientation of fibers with DEM.

2. DEM Simulation of Direct Shear Test

In the current study, a direct shear test was modeled using the commercial code PFC3D produced by Itasca (Minneapolis, USA). PFC3D is based on rigid body and soft contact approaches. In the PFC3D formulation, the spheres are rigid and cannot change their shape and size. The “soft contact” approach assumes that interaction of spheres is represented by assigning each sphere a stiffness and allowing them to overlap. The extent of overlap is used in conjunction with a contact force law to give instantaneous forces from knowledge of the current positions, orientations, velocities, and spins of the particles [30,31]. The calculation cycle of a DEM simulation are schematically represented in Figure 1.

2.1. Development of PFC3D Model

The materials, dimensions, and boundary conditions in the numerical model were based on direct shear tests conducted by Darvishi et al. [32]. Figure 2a shows the DEM model of this direct shear box for FRS with random fibers. The size of the direct shear box was chosen as 25 mm in height, 60 mm in length, and 60 mm in width (Figure 2c). The virtual shear box was enclosed by 10 rigid wall elements and divided horizontally into 2 equal halves (upper and lower). The lid of the shear box placed on top of the particles was a servo-controlled wall, and normal force was exerted on the particle assembly within the box. The rigid load plate moved up and down during compression and shearing. In addition to the 10 walls, 2 more walls (1 on the left and 1 on the right) were used for closing the horizontal surfaces on either side of the shear plane. This is a necessary step to avoid particles escaping from the box during shearing.

Sand particles were modeled as rigid balls with conventional contact elements (Figure 2b). The grain size was slightly modified and scaled up to simulate the problem effectively. In the sand matrix, a normal distribution was used for a particle radius range from 1.0 to 1.4 mm. A fiber was modeled as a series of linked particles with a bond contact algorithm, where the strength of the bond is fairly high. Each ball element in a fiber had the same diameter. The ball diameter and number of balls could be varied to match the diameter and length of the fiber being simulated. Any arbitrary fiber shape also could be constructed by changing the arrangement of the balls. In this study, we used a straight virtual fiber (Figure 2d). The fibers could be located at horizontal, vertical, and oblique directions with respect to the shear plane. Considering the test situation, fibers were randomly oriented in a direct shear box during the calibration process (Figure 2c).

The material friction coefficients could be adjusted to large and small values to achieve loose and dense setups, respectively. A very low friction coefficient was adopted to generate a dense assemblage. The assembly porosity of the test sample at the initial condition was controlled at 0.4, similar to the samples tested in the laboratory [32]. The sand–fiber system was then brought to a force equilibrium state by iteration. The friction coefficients were then reset to a specific value. A vertical displacement was then applied slowly to the top wall until the desired vertical normal loading was reached. This process formed a test specimen containing a homogeneous and well-connected granular assembly with a predetermined normal load.

When the expected normal pressure (100 or 200 kPa) was achieved and the system reached equilibrium again, the shearing stage was started. Shearing was induced by moving the upper part of the shear box at a constant velocity in a positive x-direction, while the lower section was fixed. Note that the loading speed (0.01 m/s) needed to be extremely slow to guarantee that the test was conducted under quasi-static loading conditions. During the whole shearing stage, the velocity of the top wall was adjusted using Itasca’s built-in servo-controlled functions to maintain a constant normal pressure. A horizontal velocity was applied to the upper box until the sample failed by shearing or the horizontal displacement reached 10 mm. The physical quantities, such as shear stress and horizontal displacements, were measured and recorded at the corresponding boundaries.

2.2. Modeling Contact in the Samples

In DEM simulations of granular material, the macroscopic behavior of the whole assembly may be highly dependent on the inter-particle contact law. The contact types of binary mixtures are more complicated than the single contact type in a single substance. For FRS, there are four types of contacts: within fibers, between fibers, between fibers and the sand matrix, and within the sand matrix (Figure 3). Considering the computation cost and acceptable performance, sand particles were modeled as spheres in this study. As far as we know, the particle shape has a substantial effect on the mechanical behavior of a granular solid. As an alternative approach to account for particle shape effects, the rolling resistance linear model was integrated into the DEM model to simulate sand–sand interaction. The linear parallel-bond model was applied to contacts between spheres within virtual fibers, which can reasonably capture the compression, stretching, bending, and shearing behaviors of fibers. This linear parallel-bond model has been successfully used to simulate the behavior of various fibers [33,34] and a polypropylene geogrid [35,36,37]. The linear contact model was applied to describe the interaction between a fiber particle and the sand matrix at contact. Contacts between fiber particles were also modeled using a linear contact model. The details of each contact method are described in the following subsections.

2.2.1. Linear Model

The linear model describes the behavior of an infinitesimal interface that does not resist relative rotation so that the contact moment $M_c$ equals zero. As shown in Figure 4, the contact force is resolved into linear force $F^l$ and dashpot force $F^d$. The linear force can be updated as:

$$F_n^l=k_n\times \Delta {\delta }_n,$$

(1)

$$F_s^l=k_s\times \Delta {\delta }_s,$$

(2)

where $k_n$ and $k_s$ are normal and shear stiffness, respectively, and $\Delta {\delta }_n$ and $\Delta {\delta }_s$ are the increment displacement during time step $\Delta t$ at contact normal and shear directions, respectively.

A coulomb limit is also enforced at the contact shear direction using a friction coefficient $\mu$:

$$F_s^l\le \mu \times F_n^l.$$

(3)

The dashpot normal and shear contact force between particles can be calculated as:

$$F_n^d=2{\beta }_n\sqrt{m_c{k}_n}{\dot{\delta }}_n,$$

(4)

$$F_s^d=2{\beta }_s\sqrt{m_c{k}_s}{\dot{\delta }}_s,$$

(5)

$$m_c=\frac{m_1\times m_2}{m_1+m_2},$$

(6)

where $m_1$ and $m_2$ are the masses of the contact particles and ${\beta }_n$ and ${\beta }_s$ are damping ratios in the normal and shear directions, respectively.

2.2.2. Rolling Resistance Linear Model

The rolling resistance contact model in PFC3D is based on the linear model, with the addition of a rolling resistance mechanism. The linear and dashpot forces are confirmed as in the linear model, while the rolling resistance moment is updated in the following steps:

$$M^T:=M^T-k_r\Delta {\theta }_b,$$

(7)

where $\Delta {\theta }_b$ is the relative bend–rotation increment and $k_r$ is the rolling resistance stiffness defined as:

$$k_r=k_s{\overline{R}}^2,$$

(8)

and $\overline{R}$ is the contact effective radius defined as:

$$\frac{1}{R}=\frac{1}{R^{(1)}}+\frac{1}{R^{(2)}}.$$

(9)

In Figure 4, $R^{(1)}$ and $R^{(2)}$ are the radii of contact particles. If one side of the contact is a wall, the corresponding radius $R^{(2)}$ will be infinite.

The magnitude of the updated rolling resistance moment is then checked against a threshold limit:

$$M^T=\{\begin{array}{l}{M}^T,\Vert M^T\Vert \le M^{*}\\ M^{*}(M^T/\Vert M^T\Vert ),otherwise\end{array}.$$

(10)

The limiting torque is defined as:

$$M^{*}={\mu }_r\overline{R}{F}_n^l,$$

(11)

where ${\mu }_r$ is the rolling resistance coefficient and $F_n^l$ is the normal linear force.

2.2.3. Linear Parallel-Bond Model

The linear parallel bond can act together with the linear contact model, unless the maximum stresses from the moment and force overcome the bond strength that cancels the bonding and leaves only the linear contact model between two contact particles. A graphical representation of this behavior is shown in Figure 5. The contact force is resolved into linear force $F^l$, dashpot force $F^d$, and the parallel-bond force. The contact moment $M_c$ is the parallel-bond moment $\overline{M}$. The linear and dashpot forces are updated as in the linear model. The parallel–bond force $\overline{F}$ is resolved into normal force ${\overline{F}}_n$ and shear force ${\overline{F}}_s$, and the parallel-bond moment $\overline{M}$ is resolved into torsional moment ${\overline{M}}_t$ and bending moment ${\overline{M}}_b$:

$$\overline{F}=-{\overline{F}}_n{\widehat{n}}_c+{\overline{F}}_s,$$

(12)

$$\overline{M}={\overline{M}}_t{\widehat{n}}_c+{\overline{M}}_b.$$

(13)

The processes by which contact forces, particle motion, tensile forces, and rotational moments are determined are presented next.

The normal force is updated:

$${\overline{F}}_n={\left({\overline{F}}_n\right)}_0+{\overline{k}}_n\overline{A}\Delta {\delta }_n,$$

(14)

where ${\left({\overline{F}}_n\right)}_0$ is the parallel component normal force at the beginning of the time step, $\overline{A}$ is the area of the bond, ${\overline{k}}_n$ is the normal stiffness, and $\Delta {\delta }_n$ is the relative normal displacement increment.

The shear force is updated:

$${\overline{F}}_s={\left({\overline{F}}_s\right)}_0-{\overline{k}}_s\overline{A}\Delta {\delta }_s,$$

(15)

where ${\left({\overline{F}}_s\right)}_0$ is the parallel component shear force at the beginning of the time step, ${\overline{k}}_s$ is the shear stiffness, and $\Delta {\delta }_s$ is the shear displacement increment.

The torsional moment is updated:

$${\overline{M}}_t={\left({\overline{M}}_t\right)}_0-{\overline{k}}_s\overline{J}\Delta {\theta }_t,$$

(16)

where ${\left({\overline{M}}_t\right)}_0$ is the twisting moment at the beginning of the time step, $\overline{J}$ is the polar moment of inertia of the bond, and $\Delta {\theta }_t$ is the relative twist rotation.

The bending moment is updated:

$${\overline{M}}_b={\left({\overline{M}}_b\right)}_0-{\overline{k}}_n\overline{I}\Delta {\theta }_b,$$

(17)

where ${\left({\overline{M}}_b\right)}_0$ is the bending moment at the beginning of the time step, $\overline{I}$ is the moment of inertia of the bond, and $\Delta {\theta }_b$ is the bend–rotation increment.

2.3. Calibration of PFC3D Model

The values of the micromechanical parameters of the FRS sample must be determined through a numerical calibration process, in which the macro-properties of the model are compared to the relevant measured response in the laboratory. The laboratory experimental data of Darvishi et al. [32] were used to calibrate DEM parameters in PFC3D. In the experiments, sand particles consisting mainly of quartz were clean and had mainly semi-circular, semi-angular shapes. Polypropylene fibers with specific gravity of 910 kg/m³ were used in this study, and their tensile strength was 570–660 MPa. The fibers were short (19 mm) and mixed into the dry sand at a concentration of 0.5% per dry weight of sand.

Figure 6 shows comparisons of the shear stress–horizontal displacement curves of unreinforced and reinforced sands at different normal stress levels obtained from DEM simulations and laboratory measurements. It is seen that the predicted DEM responses match well with the laboratory data, thus indicating that the proposed DEM model can capture the shear behavior of FRS. The set of micromechanical parameters selected for DEM simulation in the current analysis is shown in

Table 1. Microscopic parameters for DEM simulation.

Parameters	Value
Density of sand particle/(kg/m3)	2690
Effective modulus of sand particle/Pa	5e7
Normal-to-shear stiffness ratio of sand particle	1
Friction coefficient of sand particle	0.2
Rolling resistance coefficient of sand particle	0.25
Density of fiber particle/(kg/m3)	910
Normal stiffness of fiber particle/(N/m)	5e4
Normal-to-shear stiffness ratio of fiber particle	1
Parallel bond normal stiffness of fiber particle/(N/m)	1e11
Parallel bond stiffness ratio of fiber particle	1
Tensile strength of fiber particle/Pa	1e12
Cohesion of fiber particle/Pa	1e12
Friction coefficient of sand–fiber contacts	0.90
Friction coefficient of fiber–fiber contacts	0.90
Damping coefficient	0.7

.

3. Results and Discussion

3.1. Shear Stress–Strain Analysis

To evaluate the effect of different fiber orientations on the mechanical reinforcement, a series of numerical direct shear tests were performed on FRS using PFC3D. As shown in Figure 7, the angular orientations of fibers in the simulations were set to 0°, 30°, 60°, and 90° against the shearing surface.

Figure 8 shows the shear stress–horizontal displacement curves obtained for samples with different fiber orientations under the applied normal stresses of 200 kPa in DEM simulations. The shear resistance of FRS was not significantly affected by fiber orientation at small horizontal displacements (less than 2 mm). However, a preferential orientation of fibers induced significant shear strength anisotropy, in particular at large horizontal displacements. The highest shear stress is observed at a fiber orientation angle of 60°, and then gradually decreasing values are observed for fiber orientation angles of 90° and 30°. The specimen with a fiber orientation angle of 0° produced the weakest reinforcement effect. This result is consistent with that of Gray and Ohashi [11], who also found that an angle of 60° produced the greatest strength increase.

As shown, the fiber orientation not only affects the ultimate shear strength level, but it also changes the shapes of the shear stress–total displacement curves. Specimens with fibers having orientations of 0° and 30° showed peak shear strength followed by softening behavior. However, those with fibers having orientations of 90° and 60° showed continuous hardening behavior, indicating greater ductility.

3.2. Microscopic Observations

3.2.1. Characterization of Contact Networks and Force Transmission

Considering that intergranular forces are strongly dependent on the orientation of contacts, a complete description of load transfer in granular assemblies requires information on the distribution of contact orientations. The morphological characteristics of force chain networks (strong and weak) of FRS with different fiber orientations before and after shearing provide insights into the paths transmitting contact force in FRS. Figure 9 presents polar histograms in the XZ orthogonal view, showing the distribution of contact normals before and after shearing in specimens with different fiber orientations under 200 kPa of normal stress.

As shown in Figure 9, a similar fabric of strong normal contacts is observed in the four samples, which indicates that fiber orientations have little effect on strong contact distribution. In the initial state, strong force chains are found with their major orientation aligned in the horizontal direction, and the histogram of normal contacts has a spindle shape. In the final state, the principal direction of anisotropy of the contact normals for strong contacts is about 150°, which signifies that the strong contacts are apt to be orientated in the diagonal direction going from upper-left to lower-right of the box. The contact normals in the direction of 60° mainly decrease, so the distribution diagram transforms to a peanut-shaped distribution. Interestingly, the principal direction of anisotropy of the contact normals for strong contacts was observed to be perpendicular to the optimum orientation for maximum fiber contribution to shear strength. When the orientation of fibers is perpendicular to the main direction of strong contact force, more fibers are being subjected to radial compression, which is conducive to the formation of axial tensile forces of fibers and can contribute considerably to the strength of the composite. These results may help us to understand the strength anisotropy of FRS induced by fiber orientation.

The results indicated that the fiber orientation plays a key role in the anisotropy of the contact normals for weak contacts (see Figure 9). At the initial state, the rose diagrams of the contact normals for weak contacts are approximately spindle-shaped. The principal directions of anisotropy of the contact normals for weak contacts varied with the initial fiber angle. At the final state, the contact normals for weak contacts oriented in the direction of maximum extension decreased to a certain extent. In addition, the contact normals for weak contacts in horizontal orientations increased greatly under horizontal load.

Furthermore, the stress state also influenced the evolution of anisotropy and the force chain network [38]. A representative stress for any sub-region of the granular material can be determined by summing together the stresses in the material cells within the chosen measurement region. We defined the measurement region in the middle section along the y direction of the generated FRS sample. Figure 10 shows the measured major principal stress field in FRS with various fiber orientations under a vertical pressure of 200 kPa. It is noticed that the distribution of principal stress magnitude is very close. In the initial state, a uniform distribution of major principal stress was observed, as shown in Figure 10. In the final state, the major principal stress mainly concentrated in the direction from upper left to lower right. Cui et al. [39] and Wang et al. [40] have shown that the major principal stresses have an angle of 60° relative to the vertical in direct shear tests with the DEM. This confirms the conclusion that the maximum reinforcement can be achieved when the major principal stress direction is perpendicular to the preferred fiber-orientation plane, which is in agreement with the result by Gao et al. [41]. In principle, maximum forces are carried by contacts with orientations close to the direction of maximum principal stress [42]. This conclusion supports the finding that the achieved optimum orientation angle of fibers is perpendicular to the main direction of strong contact force in the direct shear testing.

3.2.2. Particles Displacement Vector

The particle displacement vectors can be accurately and easily captured in the DEM simulation. Figure 11 presents the total particle displacement vectors in the four samples with different fiber orientations at the end of the direct shear test. As shown in Figure 11, the direct shearing results in an arc-shaped narrow zone, which is called the shear zone in this study. Large shear distortion occurs in the shear zone, while the particles in the far field undergo limited shear distortion. Shewbridge and Sitar (1985) pointed out that the thickness of the shear zone is an important parameter that influences the fiber-reinforced shear strength. It is interesting to observe that the angle of fiber direction likely influences the shear zone thickness. As shown in Figure 11, the shear zone thickness of FRS with horizontal fibers is the smallest. The specimen with fibers inclined at an angle tends to form a thicker localized zone, and the shear zone thickness of FRS with fibers having an orientation of 60° is the largest. This can be explained by fibers inclined at an angle interfering with the development of the shear zone by transferring the shear stresses to zones farther away in which they are entangled. A larger shear zone thickness means that more particles are mobilized to impede the shearing motion, accounting for a higher strength reinforcement.

3.2.3. Average Normal Contact Force around Fibers

Identifying the causes of the FRS strength anisotropy induced by fiber orientation requires detailed analysis of the mechanisms governing interactions between sand and fibers. The interfacial shear resistance of FRS is related to friction, which could significantly prevent the fibers’ relative movement (slide) in a sand matrix and resist fiber pull-out. The fiber–sand interfacial friction depends primarily on the normal stress around the body of the fiber, efficient contact area of the fiber interface, and fiber surface roughness.

The normal force on the fibers has an important influence on their reinforcing effect. The stresses required to deform the matrix translate into an additional normal force on the matrix–fiber interface, and therefore into additional frictional and pull-out resistance. To understand the mechanism of interface interaction between fibers and sand particles, the effect of fiber orientation on the normal force of sand–fiber contact during direct shear testing is described in Figure 12. This figure makes it apparent that the normal force on the sand–fiber contact with fibers having an orientation of 60° is obviously larger than that with fibers having other orientations, contributing to the largest macro-scale shear strength. This discrepancy could be attributed to the phenomenon that the main direction of strong contact force is perpendicular to the orientation of 60°, which increased the normal force on the sand–fiber interface to the greatest extent.

3.2.4. Coordination Number for Sand–Fiber Contact

The magnitude of the interfacial friction was affected directly by the effective contact area of fibers and sand particles. Fiber reinforcement effects are more prominent for specimens with increased fiber/sand contact area [43]. An alternative parameter to the contact density is the coordination number (CN). The larger CN for sand–fiber contact indicates that larger amounts of sand particles will contact with fibers and increase the effective contact area of sand and fibers. The CN for sand–fiber contact $Z_{s-f}$ is calculated as:

$$Z_{s-f}=\frac{2C_{s-f}}{N},$$

(18)

where $C_{s-r}$ is the number of sand–fiber contacts and N is the total number of particles.

Figure 13 shows the change in sand–fiber CN that accompanies shear loading in FRS with different fiber orientations, which can provide an explanation of differences in the reinforcement effects. The results show that variation of CN for sand–fiber contacts is observed due to fiber orientation. The CN of sand–fiber contacts in the sample with 60° fibers was the highest, with successively lower CNs for 90°, 30°, and lastly 0°. This result indicates that the 60° orientated fibers give rise to the largest effective interfacial contact area of sand and fibers, which explains why the optimum orientation of fiber reinforcement is 60°.

3.2.5. Fiber Tensile Force

Fibers perform their mechanical reinforcement function by working as tension-carrying members that transfer the shear stresses in the sand matrix into tensile resistance via the interface friction along their surface. Thanks to the bond that develops between the particles, the parallel-bond contact model allows for the development of tensile forces. The tensile force along the fibers was calculated by the total tensile force at parallel-bonded contacts. Figure 14 displays the variation of the tensile force computed for all bond contacts in the fibers within the FRS samples. Figure 14 clearly shows that the mobilization of tensile force inside the reinforcing fibers is strongly influenced by shear displacement. A gradual increase of the tensile force in the fibers was recorded as the shear displacement increases up to 6 mm, which indicates the contribution of fibers progressively mobilized during the process. However, above a shear displacement of 6 mm, the tensile force is relatively stable and stays high. These results further support the idea that mobilization of the fiber reinforcement generally requires a certain horizontal displacement in the shear stress–displacement curve of the FRS.

The modeling results also clearly show that significant differences of fiber tensile force exist in FRS with different fiber orientations. Obviously, the highest tensile force was mobilized when fibers were placed at an angle of 60°, with successively lower mobilization at 90°, 30°, and lastly 0°. Interestingly, the tensile force in the specimen with horizontal fibers is nearly zero, which suggests that horizontal fibers have no effect on reinforcement. The change laws of mobilized tensile force following the change of fiber orientation is in accordance with laws of stresses and displacements in the FRS composite. This could confirm that fiber reinforcement was actually influenced by the fiber tensioning process.

Contribution of fibers to the shear strength is effective when FRS are subjected to tension, although the fibers play no role in compression as they may buckle or kink [44]. Figure 15 shows the spatial distributions of the contact force in fibers for different shear displacements. The figure shows that fibers may be subject to tension (yellow) or pressure (blue) in FRS with different fiber orientations and horizontal displacements. To analyze the difference quantitatively, Figure 16 shows the proportion of numbers of tensile contacts to the total number of contacts in fibers. Obviously, most of the fibers in different samples are subjected to tension at the initial state, although the value is very small. The number of tensile contacts in fibers with fiber orientations of 0° and 30° first decreased and then increased to a stable value. However, the number of tensile contacts in FRS with fiber orientations of 90° and 60° always stayed high during the whole progression. The FRS with fiber orientation of 60° had the highest numbers of tensile contacts, which contributes to the understanding that the fiber initial orientations of 60° mobilized the largest tensile force.

3.2.6. Energy Storage and Dissipation

The energy dissipation mechanism provides an important way to explain the micro-mechanical mechanisms in granular materials. To study the effect of fiber orientation on reinforcement, it is necessary to study the evolution of individual energy components in the sand–fiber system under direct shear load. Figure 17 shows the evolution of four energy components, namely bond strain energy, friction energy, strain energy, and kinetic energy during the shear process under normal stress of 200 KPa. Fiber deformation leads to the storage of strain energy in the parallel bonds that join the fiber balls to one another. As shown in Figure 17a, the effect of bond strain energy storage in the fibers gradually appears with the development of direct shear, which also means that the mechanical properties of fibers have been developed gradually. In addition, the bond strain energy stored in fibers having an orientation of 60° is the largest, which means that the mechanical properties of fibers having an orientation of 60° can be mobilized to the greatest extent. Figure 17b indicates that the friction energies increase monotonically with the increase of horizontal displacement. Furthermore, the friction energy dissipation in the sample with fibers having an orientation of 60° is higher than that in samples with other orientations. The stored strain energy in the FRS material increases first and then remains stable during the shearing process (see Figure 17c). The stored strain energy in the sample with fibers having an orientation of 60° is the largest. The enhanced capacity of the assembly to store strain energy could also be a factor in enhancing its shear strength. Because of the low loading rate, the kinetic energy is very small and almost unchanged during the shearing process (see Figure 17d).

The energy dissipation analysis can explain the mechanisms of optimizing the reinforcement performance in the presence of different fiber orientations. Numerical results show that friction energy and strain energy play a significant role in energy storage and dissipation. The high reinforcement effects of fibers having an orientation of 60° could be attributed to an increase in the energy storage/dissipation potential, which enables the particle system to transfer and bear loads in a relatively stable state.

4. Conclusions

In this study, a DEM approach was used to analyze the microscopic mechanisms underlying the strength anisotropy of FRS subjected to direct shear load. The following conclusions can be drawn based on the present simulation. When the orientation of fibers is perpendicular to the strong contact force, the reinforcing effect of the fibers to the composite strength is most efficient. This preferred direction is also perpendicular to the direction of the major principal stress. The orientation of the fibers can also affect the particle displacement field. The fibers with an orientation of 60° also cause the largest increase in normal stress around the fiber body and the effective contact area, which leads to an increase in sliding friction between fibers and the composite matrix. The modeling results also clearly showed that significant differences in fiber tensile force exist for different fiber orientations. The highest tensile stress was mobilized in 60° fibers, which can confirm that fiber reinforcement was mainly due to the tension process. The high reinforcement effects of fibers having an orientation of 60° can be attributed to an increased energy storage/dissipation potential, which enables the particle system to transfer and bear loads in a relatively stable state. This work will generate fresh insight into the fiber-reinforcing mechanism in FRS.