Discrete Element Analysis of the Load Transfer Mechanism of Geogrid-Ballast Interface under Pull-Out Load

Geogrids have been extensively used in subgrade construction for stabilization purposes of unconfined ballast. Based on wellcalibrated microparameters, a series of geogrid-reinforced ballast models with different geogrid sizes and particular structures were developed to reproduce the mechanical behavior of the geogrid under pull-out load in this paper. And the rationality of the DEM model is verified by comparing the evolution law pull-out force measured by laboratory tests and numerical simulations under comparable conditions. Moreover, themacro pull-out force and the internal force distribution of the geogrid were analyzed, and the contact force statistical zones of the particle system were divided accurately according to the results. Meanwhile, both the force transfer mechanism in the geogrid-ballast interface and the sectionalized strain of the geogrid were discussed. And results unveil that the pull-out load is transmitted along the longitudinal ribs to the transverse ribs, and nearly 90% of the load is transmitted to the contact network (in statistical zone 1) in front of the first transverse rib, resulting in strong interlocking between the particles occurs in statistical zone 1. And the second transverse rib is the strength dividing line between strong and weak contact forces. +en, additional pull-out tests on the control groups were conducted, and the sectionalized strain of the geogrid and the peak pull-out force, as well as the energy dissipation were systematically analyzed. In addition, the proposed method used in simulation holds much promise for better understanding of the reinforcement mechanism and further optimizing the performance of geogrid-reinforced structures.


Introduction
Due to its excellent toughness and unique mesh structure, geogrids are widely used in geotechnical engineering [1][2][3]. Generally, granular geomaterials always possess considerable compressive strength and shear strength, but low tensile strength. e overall stability of the soil can be improved by combining the outstanding tensile strength of the geogrid with the compressive strength of the soil [4][5][6]. Especially, for granular ballasted track bed, a considerable part of the internal tensile stress is transmitted through the interaction of the geogrid and the aggregates so that the stress in the track bed is efficiently diffused [1,7]. en, the lateral deformation of the track bed can be effectively controlled by embedding the geogrid in the ballast.
With the extensive application of geogrids, the specific structure of geogrids has also attracted wide attention. For the geogrid-stabilized granular aggregates structure, the reinforcement performance is not only provided by the interface friction between the ribs and aggregates but also mainly by the interlock between the grid and the particles [8]. erefore, the reinforcement performance is related not only to the friction characteristics of the geogrid-aggregates interface but also to the grid area ratio in the reinforcedplane and the thickness of the ribs. e corresponding parameters of the geogrid can be characterized as aperture size, rib thickness, node height, and aperture geometry. e design criteria of reinforced soil generally stipulate the tensile strength and pull-out strength of geogrids, as well as the friction strength of the geogrid-soil interface. In addition, Hashash and Yoe [9] and Chen et al. [10] demonstrated that pull-out test is one of the most effective methods to study geogrid-soil interaction behavior and evaluate the reinforced performance of geogrid. Many researchers have conducted a lot of valuable experimental [11][12][13][14] and numerical [15][16][17][18] studies on these characteristics. In all the mentioned works, the achieved results confirmed the fundamental role of geogrid in enhancing the stability of the granular materials and found that the geogrid-soil interaction mechanism is highly complex. Among the experimental approaches, Moraci and Recalcati [14] reported that when the confining stress of the pull-out test is greater than 25 kPa, the strain hardening behavior occurs. Indraratna et al. [7] presented experimental evidence based on process simulation test that the deformation and the breakage of the ballast can be alleviated by embedding the geogrid in a suitable placement location. Moreover, the reinforced performance is influenced by the geogrid type. In terms of numerical simulation, Hussein and Meguid [17] adopted a three-dimensional finite-element method to study the pull-out behavior of biaxial geogrid embedded in granular backfill material and conducted that the bearing resistance contributes more to the total pull-out capacity than the frictional resistance. Based on the pull-out test by the two-dimensional discrete element method, Chen et al. [19] analyzed the influence of the reinforcement length and thickness of the geogrid on the reinforcement performance and presented that the displacement of the geogrids in the pull-out tests showed a nonlinear decrease along the longitudinal ribs.
All these research studies have promoted the understanding of the geogrid-soil interface behavior. Nevertheless, the influence of the specific structure of geogrids on the coarse-grained soil system and interaction of mesomechanism between coarse aggregate and various geogrid inclusions remain ambiguous. e contact between the particles of the granular material is extremely sensitive to the induced load and particle movement. Even slight changes of the external load can cause great changes in the force chain and the contact network [20]. And the involvement of geogrids makes the contact behavior of the reinforced system more complicated. However, research studies at the experimental level cannot accurately monitor the contact behavior of granular materials, and it is difficult for the finite-element method to reflect the mechanical behavior of discontinuous media. As a discontinuum-based numerical method, the three-dimensional DEM has been proven and widely used to look insight into the mesoscale mechanical response of geogrid-ballast interface [10,[21][22][23][24].
In this study, based on the well-established constitutive models in DEM, the geogrid is modeled by regularly arranging unit spheres through parallel bonding bonds to reproduce its particular mesh geometry and mechanical behavior. And 4 sizes of rectangular geogrid and 1 size of triangular geogrid models were established to simulate the pull-out behavior of the geogrid embedded in the ballast. According to the evolution of the axial force of the geogrid along the longitudinal ribs, the contacts in the shear zone surrounding the geogrid is reasonably distinguished.
Meanwhile, the contact force fabric and the load transfer mechanism of the geogrid-ballast interface under the pullout load are analyzed. In addition, in light of the development of peak pull-out force and the evolution of frictional energy consumption, the influence of the joint protuberance and aperture size on geogrid reinforcement performance are also emphasized.

DEM Modeling
2.1. Modeling of Ballast and Geogrid. Generally, the particle size of the ballast is 10-60 mm [1]. e interaction between geogrid and coarse grains is significantly affected by the angular characteristics of the particles [22,23]. To balance the computational efficiency of numerical simulation and the angular characteristics of the particles, a tri-ball particle developed by Miao et al. [23] being proven to reflect the ballast rotation and interparticle occlusion behavior in PFC3D was used in this study. Specifically, this ballast particle model was generated by three equal-sized unit spheres through "clump" logic, as presented in Figure 1(a). Besides, the interlocking effect between geogrid and aggregates is strongly influenced by the relationship between the ballast size and the aperture size [25][26][27]. Considering the influence of the geogrid-ballast size relationship, the current model adopts the ballast particle size distribution suggested by Indraratna et al. [28]. e grading curve for the ballast is shown in Figure 1 e modeling method of the geogrid proposed by Ngo et al. [15] was adopted in the present study. e size and geometry of the geogrid were adjusted based on this method. Correspondingly, the microparameters of the ballast and geogrid, well calibrated by Ngo et al. [15], are also used in the current model (as summarized in Table 1). Specifically, the geogrid model consists of spherical elements, 2 mm in radius, bonded together by parallel bonding, and the geometry of the aperture can be flexibly set by specifying the arrangement direction of the sphere elements ( Figure 2). In addition, the nodes of the biaxial and triaxial geogrids were generated by the bonding of 5 and 7 clump elements, respectively. e specific structure of the clump and the joints of the geogrid are shown in Figure 2. Table 2 shows the size and geometry of the 5 geogrid numerical specimens in this study. For each geogrid specimen, a pull-out test was conducted at 20, 30, 40, and 50 kPa normal stress, respectively. Except for the different geogrid specimens, the steps for establishing the numerical pull-out test are essentially the same in all cases. erefore, only the specimen BG1 under 20 kPa normal stress is used as an example to briefly describe the modeling procedures of numerical tests.

Setup of Numerical Pull-Out Test in DEM.
In order to avoid damage to the geogrid caused by the precompression of the initial ballast specimen, the initial ballast particles were divided into two layers to be generated (each layer is 200 mm high, 300 mm long, and 300 mm wide). en, the geogrid model was generated in the 6 mm high space between the two ballast layers. In the plane view, the geogrid is centered in the model, and the distance to the model boundary is 86 mm, 142 mm, 86 mm, and 50 mm, respectively (Figure 3(a)). Since the initial sample ballast particles partially overlap, the two servo mechanisms were used to control the upper and lower walls of the two ballast layers to precompress the ballast sample to achieve equilibrium, respectively ( Figure 3(b)). e two walls above and below the geogrid were deleted so that contacts between the geogrid and the ballast were allowed, and the servo mechanisms of the walls were redefined to apply 20 kPa normal stress to the reinforced specimen. In the positive direction of the yaxis (the pull-out direction), a constant velocity of 0.375 m/s is applied to the clamped end, as presented in Figure 3(c). e key indices are recorded during the pullout process until the displacement of the geogrid end reaches 90 mm.
It is noteworthy that the inherent properties of the material modeled depend on the microparameters not relating to the arrangement of the unit spheres (i.e., the aperture geometry and size of the geogrid) of the model in PFC3D. Miao et al. [23] also demonstrated that the stress-strain curves of geogrids with the same microparameters but a different aperture geometry (the biaxial and triaxial geogrid with 40 mm aperture) are similar through the single-width tensile tests.

Pull-Out Behavior and Axial Force Distribution of Geogrid.
e results of the numerical test for BG1 and TG1 are first analyzed in this section to demonstrate the general pull-out behavior of geogrid embedded in ballast. e pullout force is one of the macroindicators to evaluate the pullout behavior of the geogrids. Figure 4 reveals the pull-out force against the pull-out displacement of BG1 and TG1 under different normal stress. Clearly, before the displacement reaches 15 mm, an approximately linear increase in each curve is observed, and its growth rate is positively correlated with the normal stress. It also can be seen that the higher normal stress corresponds to a higher peak pullout force, and the displacement corresponding to the peak    is relatively lagging. Moreover, the curve of TG1 under the same normal stress fluctuates more than the curve of BG1, and strain hardening of the pull-out curve of TG1 under 50 kPa normal stress (within 90 mm of pull-out displacement) occurs. In the pull-out test of geogrid-reinforced coarse soil, if the normal stress is not large enough (greater than 25 kPa), the relatively stable peak of the pull-out curve may not be observed. As the normal stress increases, the displacement required for the peak pull-out force to appear is greater. Stahl et al. [29] conducted a series of laboratory pull-out tests of geogrid (biaxial geogrid was pulled out by 30 mm). Figure 5 presents the comparison of the laboratory tests by Stahl et al. [29] and the numerical simulation results in this study. Correspondingly, the evolutionary trend of the pull-out force with displacement is consistent with the previous analysis. e correlation between normal stress and pull-out behavior may be the cause of the differences in the curves. In general, the numerical results are basically consistent with the laboratory test results and the comparison with the experimental results also further validates the rationality of the current model.
In order to explore the internal force transmission of the geogrid under pull-out load, Figure 6 reveals the axial force distribution of BG1 under different normal stresses. It should be noted that because the external geogrid of the reinforced specimen is the clamping end and there is no contact with the ballast particles, the figure only shows the axial force distribution and deformation of the geogrid embedded in the ballast layer. e axial force is transmitted along the longitudinal ribs from the clamping end to the depth of the ballast (negative direction of the y-axis), and a stepwise drop occurs when passing the transverse rib. Moreover, when the pull-out force reaches the peak value, the contribution of the transverse ribs at different positions to the pull-out resistance is different. Specifically, the axial force at the first transverse rib decreases the most. Previous studies reported [23,27] that the axial force distribution of the triaxial geogrid is roughly the same as the biaxial geogrid in Figure 6. Based on the results shown in   Figure 6, the geogrid is divided into 4 segments (the length of each segment of the longitudinal rib is 36 mm) with the transverse ribs as the limit. And in order to facilitate understanding, these segments of geogrid are numbered 1, 2, 3, and 4, respectively to indicate their position. To illustrate the contribution of the transverse ribs at different positions to the pullout resistance, F i /F is defined as the percentage of axial force,where F is the peak value of the pull-out force and F i is the sum of the axial forces of the longitudinal ribs in the segment I, as shown in Figure 7. It can be seen that even if the normal stresses are different, the axial force of the segment 1 of the geogrid accounts for about 88%-90% in all conditions, and the segment 2 to 4 decrease to 23-37%, 12%-19%, and 4%-8% in sequence. When the pull-out force is close to the peak value, the transmission method of the axial force in the geogrid is not affected by the normal stress, but the axial force increases significantly with the increase of the normal stress level.

Contact Forces and Fabric Analysis.
e contact force chains of granular materials are extremely sensitive to the change of local force, and the slight fluctuation of external load is enough to change the force chain structure greatly. e complex contact network is formed by the continuous contact among the particles, and the contact force is transmitted along with the special structure of the chain path (i.e., the force chain) in the contact network [20]. e analysis of contact force evolution in granular systems is an important method to investigate the interaction mechanism of geogrid-soil interface. Figure 8 shows the projection of the normal contact force chain on the y-o-z plane for BG1 and TG1 under 20 kPa normal stress. When the pull-out displacement is 0, the distribution of contact force chains is relatively uniform. As the geogrid is gradually pulled out, the force chains develop in front of each transverse rib, and the strength of the force chains at the front end of the first transverse rib are greatest. e results suggest that the force chain distribution agrees well with the axial force distribution of the geogrid. e interlocking among the particles mainly occurs between the ballast surrounding the transverse ribs. Since each transverse rib of BG1 is perpendicular to the pull-out direction, the bar-shaped strong force chains in Figure 8(d) are separated by apertures. And according to the division of the geogrid segments in the previous section, the contact force of the main shear bond is also divided into four statistical zones, which correspond to the geogrid segment number, as revealed in Figure 8   Advances in Civil Engineering e Fourier series approximation (FSA) method proposed by Rothenburg et al. [30] was used to describe the anisotropy of contact force distribution quantitatively. e specific fitting function is as follows: where f n is the distribution of average normal force density; f 0 is the average normal contact force in statistical each zone; a n is the coefficient of normal contact force anisotropy; and θ and θ n are the direction of normal contact force and the principle direction of anisotropy, respectively. Figure 9 presents polar plots and FSA of projection of normal contact force on the y-o-z plane in four statistical zones of BG1. Each polar plot is divided into 36 angular bins of 10°a nd the mean value of the normal force in each bin is used to plot. e number of contacts in each statistical zone is basically distributed at the same level (728-798). With the change of the statistical zone (from Figure 9(a) to Figure 9(b), i.e., moves along the negative direction of the y-axis), the amplitude of the contact force gradually decreases, and the anisotropy coefficient a n is 0.76, 0.68, 0.36, and 0.45, respectively. e principle directions of anisotropy from statistical zone 1 to statistical zone 4 are 41.34°, 44.34°, 29.01°, and 0.01°, respectively. e main transmission path of the pull-out load is the longitudinal ribs-transverse ribs-surrounding ballast contact network. Besides, the evolution of the contact network is accompanied by friction between particles.
In order to distinguish the strength of the normal contact force chain, the percentage of the average contact force in each statistical zone P f is calculated as follows: where f a is the total average of the normal contact forces in the four statistical zones and f 0 i is the mean value of the normal contact forces in the statistical zone i. e strong contact force    chain is defined as the set of forces greater than f a , otherwise, it is the weak contact force chain. Figure 10 shows the variation of the contact force strength of the particle system and the axial force of the geogrid with the statistical zones. Overall, the average normal contact force chain strength changes from a strong contact network to a weak contact network near the second transverse rib (i.e., statistical zone 2). Under different normal stresses, the curves of the percentage of axial force and that of normal contact force are approximately the same, and the differentiation occurs after the second transverse rib. is is because the pull-out resistance also includes friction resistance.
It should be noted that the values of the pull-out displacement in all conditions are the same as in Figure 6, and the corresponding pull-out force is close to the peak value rather than the actual peak value of the pull-out force. Even so, the axial force of the longitudinal ribs and the normal contact force of the particles are believed to represent the pull-out force and the interlocking behavior between the particles within a reasonable error range.   Advances in Civil Engineering

Influence of Pull-Out Direction of Rectangular Geogrids.
Dong et al. [18] proved that the tensile performance of the biaxial geogrid is closely related to the tensile direction and reported that the highest tensile strength occurs in the direction of the ribs. Considering that the number of transverse ribs and the number of longitudinal ribs in BG1 or BG2 are different, it is not suitable for comparison of different pull-out directions. en, the pull-out test model for BG3 was adjusted. Specifically, the inside of the reinforced sample is exactly the same, and the pull-out direction is adjusted from the positive y-axis direction to the positive x-axis direction. e adjusted specimen is described as BG3′ for facilitate expression. e sectionalized strain of BG3 and BG3' (i.e., the strain from the first transverse rib to the last transverse rib) inside the test box is obtained. Figure 11 shows the sectionalized strain of geogrid versus the pull-out displacement extracted in simulations. In the initial stage of pull-out process (0 < pull-out displacement < 15 mm), the normal stress and sectionalized strain increment Δε negative correlation. However, the overall strain of the geogrid is positively correlated with the applied normal stress under pull-out loads, as demonstrated by Alagiyawanna et al. [31] and Miao et al. [23]. It can be inferred that the main strain occurs before the first transverse rib. As the pull-out displacement changes, the interaction between the geogrid and the ballast is a dynamically changing process. High normal stress promotes the geogrid to mobilize the particles to form a stronger interlocking effect. In addition, after the pull-out displacement reaches 20 mm, the strain curves of the geogrid under various normal stresses fluctuate. Combined with the results of the   previous sections, these fluctuations are inferred as the stages where the rearrangement of ballast and reorganization of the bearing structure of the particle system happens. Corresponding to the load transfer mechanism, the strain of the geogrid is also gradually transferred from the pull-out end to the anchored end. Due to the different side lengths of the longitudinal and transverse ribs, pull-out along the direction of the shorter rib will cause a small reduction in the strain of the anchored segment.

Influence of Joint Protuberance and Aperture
Size on Geogrid Reinforcement Performance. As an important structure for connecting geogrid ribs, the strength and thickness of geogrid joint are generally greater than ribs. In the reinforcement of fine-grained materials, the protuberance of the geogrid nodes can provide part of the frictional resistance and uniformize the internal force distribution of the ribs [32,33]. In geogrid-reinforced ballast systems, the interlocking effect of geogrids and particles is dominant compared to the interface friction. In this section, the influence of the node protuberance on the reinforcement performance of the geogrid was evaluated using the peak pull-out force as an indicator.
Two additional models, described as BG1 * and TG1 * , were established to simulate the geogrids without node protuberance by replacing the clump nodes of BG1 and TG1 with spherical particles (4 mm in diameter). And the comparison of peak pull-out force with and without node protrusions of BG1 and TG1 are shown in Figure 12. It can be observed that the peak pull-out force of the geogrid with or without node protuberance under 20 kPa normal stress is very close. At 30 kPa normal stress, the peak pull-out force of BG1 is significantly higher than BG1 * , but the peak pull-out force of TG1 is lower than TG1 * , which may be due to differences in mesh geometry and aperture size.
Furthermore, compared with the geogrid without node protuberance, the peak value of the biaxial geogrid (BG1) under 40 kPa and 50 kPa normal stress increased by 13.6% and 17.1%, respectively, and the triaxial geogrid (TG1) increased by 8.8% and 8.6%, respectively. Overall, taking peak pull-out force as the evaluation index, the node protuberance under high normal stress is more effective to improve the reinforcement performance of the geogrid.
In view of what have been discussed above, Figure 13 reveals the frictional energy dissipation of the particle system against the pull-out displacement under 50 kPa normal stress. e friction energy evolution of the particle system can be used as a method to quantify the ability to mobilize ballast particles of the geogrid. e geogrid mobilizes the ballast particles sequentially from the pull-out end to the  anchored end. Meanwhile, the relative displacement between the particles, including rotation and movement, can increase the friction energy consumption. Due to the large particle size of the ballast particles (10-60 mm), the number of nodes has little effect on frictional energy consumption. And the mesh area (effective area of geogrid) has a significant influence on the friction energy. In the case of the same reinforcement area (13824 mm 2 ), the friction energy consumption of BG4 is greater than others. And until the pullout process is terminated, the friction energy consumption of the triaxial geogrid (TG1) with 40 mm aperture is minimal. In general, corresponding to the macropull-out force, the friction energy consumption in the initial stage (0 < pullout displacement < 15 mm) also increases roughly and linearly. e joint protuberance may help to improve the plane rigidity of the geogrid to maintain its original geometry.

Conclusions
In this paper, based on well-calibrated microparameters, a series of geogrid-reinforced ballast models with different geogrid sizes and specific structures were developed to reproduce the mechanical behavior of the geogrid under pull-out load. e macro pull-out force and the internal force distribution of the geogrid were analyzed, and the contact force statistical zones of the particle system were divided accurately according to the results. Meanwhile, both the force transfer mechanism in the geogrid-ballast interface and the sectionalized strain of the geogrid were discussed. Furthermore, the influence of node protuberance and aperture size on the reinforcement performance of the geogrids was systematically analyzed based on the evolution of peak pull-out force and frictional energy consumption of the particle system. e main conclusions of this study can be summarized as follows: (1) e main transmission path of pull-out load is the longitudinal ribs-transverse ribs-surrounding ballast contact network. Specifically, the pull-out load is transmitted along the longitudinal ribs to the transverse ribs, and nearly 90% of the load is transmitted to the contact network in front of the first transverse rib, resulting in strong interlocking between the particles occurs in statistical zone 1. Moreover, the second transverse rib is the strength dividing line between strong and weak contact forces. e strong contact force, defined by the average normal contact force, is mainly distributed in front of the second transverse rib. (2) In the initial stage of pull-out process (0 < pull-out displacement < 15 mm), the overall strain of the biaxial geogrid growth as the normal stress increases, but the normal stress and sectionalized strain (reinforced area with complete mesh in ballast layer) of the geogrid increment Δε are negatively correlated. e development of strain increment corresponds to the force transmission mechanism of geogrid under pull-out load. (3) Based on the analysis of peak pull-out force and friction energy consumption of the particle system, the contribution of the node protuberance to the peak pull-out force of the geogrid is inferred that the node protuberance is beneficial to maintain the plane rigidity of geogrid. Appropriate plane stiffness enables the interlocking effect between the geogrid and the ballast particles to be fully utilized to resist the induced loads. And the aperture size and reinforced-plane area are the main controlling factors for reinforcement performance.

Data Availability
e data used to support the findings of the study are available from the corresponding author upon request.

Conflicts of Interest
e authors declare that they have no conflicts of interest.