Abstract
The discrete element method (DEM) was used to study the behavior of crushable angular gravel in the cyclic soil-structure interface test. Two shapes of agglomerates were simulated by filling two scanned angular gravels with spheres connected by bonds that were given the shear and normal strength complying with Gaussian distribution to simulate random flaws. The proportion of these two shapes to constitute a numerical sample was named composite pattern. Good agreement in terms of macromechanical behavior between DEM simulation and laboratory test results has been attained. Agglomerate breakage is deeply influenced by the interface shearing behavior and mainly occurs on the interface and the space nearby. Graphs of interface after shearing are introduced to directly and clearly reflect microbehavior of breakage. The evolution of microstructure including anisotropies and coordination numbers is significantly influenced by normal stress and agglomerate breakage, and composite pattern determines the magnitudes of shear force anisotropy and coordination numbers. The evolution of contact orientation distribution is the forming cause of the “adjustment phase,” during which once the shearing direction changes, the values of contact normal anisotropy and normal force anisotropy will slump to their nadir and then rise back again.
1. Introduction
The soil-structure systems such as cutoff walls, embankment dams, and deep foundations are widely existent in the world. Many soil-structure interface tests have thus been performed to investigate the macrobehaviors such as the stress-displacement relationship (e.g., [1–6]). Basing on these laboratory test results, researchers proposed mathematical formulations and models to describe the macromechanical relationship. The Clough and Duncan model [7] was proposed to describe the nonlinear shear stress-shear displacement relationship for interfaces relying on direct shear tests. Desai et al. [2] proposed that a function of normal stress, tangential displacement, number of loading cycles, and initial density could be adopted to describe the shear stress of a sand-concrete interface. Brandt [8] proposed a rigid plasticity model of an interface, which could be used to describe the cyclic behavior of the interface. The “Disturbed State Concept” was introduced by Desai and Ma [9] to describe the damage of soil during shearing, and such concept presented the constitutive model characterizing the behavior of a soil-structure interface. Zhang et al. [4–6] performed substantial tests under different conditions and proposed the “Elasto-Plasticity Damage Interface model” (EPDI) using the finite element method (FEM) and considering the damage to surrounding soil under the application of different loadings. The model could effectively and accurately anticipate and analyze the macroscopic evolutional process of interfaces between structures and gravels or other coarse-grained soils. Although these methods mentioned above could achieve numerical macromechanical results that have good agreement with experimental results, they fail to investigate the soil-structure interface tests in a microscopic perspective due to the limitation of methodology.
As the discrete element method (DEM) has now developed into a mature and powerful tool for studying granular materials, it can be of great help to study the microscopic characteristics in numerical simulations and has been widely applied in simulating and studying a number of geotechnical fields, such as biaxial tests, triaxial tests, postliquefaction, and rockfill creep (e.g., [10–16]). DEM could also be employed to simulate the breakage behavior of granular materials, and two main methods have been developed. (i) Large agglomerates are generated by bonding a number of disks for 2D cases or spheres (also named subparticles or microparticles in a lot of literature) and for 3D cases, where there are two techniques to arrange the spheres within agglomerates: hexagonal close packing (HCP) and radius expansion. The HCP agglomerate, which can theoretically generate the densest packing for uniformly sized spheres but are created highly anisotropic [17], has been adopted by many researchers to study the particle breakage behavior [10, 18–20]. The radius expansion technique, introduced by Potyondy and Cundall [21], is to iteratively expand the radii of particles which have less than three contacts. Alshibli and Cil [22] combined synchrotron microtomography (SMT), which could help to locate positions and determine the radii of particles, and the radius expansion technique to quickly and substantially generate crushable agglomerates in the numerical sample to investigate the 3D evolution of fracture behavior of sand subjected to 1D compression. There is one special arrangement adopted by Alonso et al. [23], who used 14 spheres with identical radii to form an agglomerate of pyramidal shape to characterize sharp edges and relatively flat faces. (ii) Large particles are replaced with predefined finer grains when a predefined particle strength is reached (e.g., [24–26]). Tapias et al. [27] used 13 subparticles with identical sizes to substitute one spherical particle when one comminution breakage occurred, and Ciantia et al. [28] adopted Apollonian filling to establish five fragmentation multigenerational patterns. Cil and Buscarnera [29] proposed a replacement pattern comprising 20 fragments and conducted oedometric compression tests with different grain sizes to investigate the grain-scale effect of fracture.
Anisotropy is an important characteristic of microstructure, reflecting the fabric composition of granular materials and has attracted many researchers (e.g., [16, 30–35]). Oda [31] derived and generalized the definitions and formulae of fabric tensor. Rothenburg and Bathurst [32] proposed an approach to partition the stress tensor for 2D assemblages and established a relationship between the mobilized angle of friction, the most common measure of shearing resistance and three microstructural anisotropy parameters. Ouadfel and Rothenburg [33] generalized the “stress-force-fabric” relationship in the 3D form and set up a connection between shear strength and anisotropies. Zhao and Guo [34] and Guo and Zhao [35] adopted a three-dimensional DEM to explore the characteristics and evolution of shear-induced anisotropies during the shearing process and successfully built a series of analytical formulae to characterize anisotropies in three characteristic states: liquefaction, phase transformation, and the critical state. Gu et al. [36] and Zhang and Wang [11] used Guo and Zhao’s results [35] to explore the evolution of microstructural anisotropies in soil under shearing and in sand during triaxial creep, respectively. The coordination number (i.e., the average number of contacts per particle) is another important microstructural characteristic and regarded as an assistant role within microstructure in linking macrobehavior and microbehavior. Rothenburg and Kruyt [37] showed that the anisotropy of contact orientation affects the relationship between the void ratio and the coordination number. Gu et al. [36] found that the initial coordination number may be an intrinsic variable reflecting the soil density or liquefaction potential. There are some literature focusing on the effect of breakage on anisotropies or the coordination number in that particle breakage results in lower peak values of anisotropies [20, 38, 39]; Zhou et al. [39] also found that different breakage modes have a minor effect on the tendencies of evolutions of anisotropies but affect the degree of the reductions of anisotropies. Mcdowell and Bolton [40] proposed that the probability of breakage reduces with an increase in the coordination number.
In this paper, we employed a three-dimensional DEM to study the evolution of agglomerate breakage and microstructure of angular gravel under the interface shearing behavior. Two shapes of agglomerates were adopted, and composite pattern was defined as the proportion of these two shapes to constitute a numerical sample. However, such concept, i.e., composite pattern, or similar one has not been found in relevant literature. There are two additional contrast simulation groups, one using breakable agglomerates with different composite patterns under identical constant normal stress and another being unbreakable agglomerates with the same composite pattern under three different normal stress conditions. The effect of composite pattern on agglomerate breakage and microstructural evolution as well as the effect of breakage on the evolution of microstructure is considered and investigated by those two contrast simulation groups. “Adjustment phase,” occurring after the alteration of shearing direction immediately, which has also not been found in other literature, means a phase, during which the values of contact normal anisotropy and normal force anisotropy will slump to their nadir and then rise back again. The paper is organized as follows. Section 2 describes the methodology and procedure of the simulations. Section 3 presents the results on macromechanical responses, and the next section characterizes and quantifies the evolution of agglomerate breakage. Section 5 mainly explores the evolution of microstructure including anisotropies and coordination numbers. This section also investigates the cause of formation of “adjustment phase,” and the last one (i.e., Section 6) provides major conclusions of the study.
2. Numerical Simulation of the Soil-Structure Interface Test
2.1. Modelling of Crushable Angular Gravel
The commercial discrete element code PFC3D 5.0 was used to perform the numerical simulations in this research. Adopting the experimental grain size will lead to the generation of too many soil particles, which significantly undermines the computational efficiency and adds to computational time substantially. Evans and Valdes [41] proposed that simulation with mass scaling (i.e., amplifying particle sizes) showed good agreement with laboratory samples in terms of macroscale response. They further proposed that it could reasonably mimic microscale evolution in granular materials, supported by Evans’s previous research [42]. There are some other literature adopting mass scaling for improving calculation efficiency and achieving desirable results (e.g., [20, 43–46]). The numerical grain sizes in this paper were five times those in laboratory interface tests [4], and they had the same trend in terms of grain size distribution (GSD), as seen in Figure 1.
This paper attempts to reach a balance between accuracy and efficiency through adopting as few spheres as possible to form an agglomerate so that enough number of agglomerates could be generated to constitute a numerical sample without requiring too long computational time and effort. Although an agglomerate with a very large number of spheres was more like rock or soil in terms of morphology [47], it is verified that agglomerates with a limited number of spheres were able to generate reasonable results [20, 22, 27, 48, 49]. As observed from Figure 2, asperity, angularity, and sharp edges are obvious characteristics of angular gravels, and those gravels could be categorized into two kinds: one was slender and another relatively more rounded. For modelling these characteristics, two representative angular gravels, which came from the original place of experimental gravels, Changping, the suburb of Beijing [4], were picked out. Each of these two gravels represented one kind of gravel. They were scanned, and then the two corresponding 3D models were generated, i.e., the upper left one of Figure 3(a) and the upper right one of Figure 3(b). The two 3D models were saved as STL files for DEM simulations and then imported into PFC3D 5.0, where they were set as closed walls, which were filled with spheres. Two agglomerate shapes were thus determined. Shape-A agglomerate consists of 14 spheres, and shape-B agglomerate consists of 12 spheres (Figure 3). The former one was slenderer and more oblate than the latter one while the latter one was relatively more rounded with a higher level of sphericity, both of them successfully simulating characteristics of angular gravels in the experimental sample (Figure 2). Taking the major axis length as the size of the agglomerates, the sizes of the original agglomerates ranged from 25 mm to 50 mm.
(a)
(b)
2.2. Simulation Setup and Model Parameters
According to laboratory interface tests [4], the length and width of the experimental sample were 45 cm and 36 cm, respectively. The range in height of the influence of soil-structure interface on soil was normally no more than five times the average size of the particles (3.5 cm). Therefore, the numerical sample was represented by a cubic space of 45 L · 36 W · 17.5 H cm. The structure material in this test was rough steel with sawtooth (Figure 4(a)). In the numerical simulation, the top wall with sawtooth was used to model the rough steel. Five rigid walls and the top wall made up the box (Figure 4(b)). The density of sphere was 1750 kg/m3, the same as gravel in laboratory tests [4]. The initial porosity of the numerical sample was approximately 0.45. The bonded agglomerates contained intragranular voids, which thus led to a higher porosity than the real value [20, 50].
(a)
(b)
In the laboratory interface tests, the experimental sample generation process was divided into five same and repeated work, and this process was adopted in generating the numerical sample. Each time the specific weight of sample obeying GSD was generated in the designated layer, followed by repeated compaction, until the desired thickness (3.5 cm) was achieved. After the completion of sample generation, each pebble in the “clump” was replaced by a sphere with identical volume. Figure 4(b) exhibits a numerical sample of 4257 agglomerates, comprising 53656 spheres, initially color coded by five layers; shape-A agglomerates constitute 15% of the sample, while shape-B agglomerates account for 85%. Before the numerical sample was cycled to an equilibrium state, parameters were assigned to spheres and contacts in DEM simulations, and these parameters will be discussed later. Then, the numerical sample was brought to a constant normal stress of 200, 400, or 700 kPa by the “servo control” of the top wall and kept by the constant normal stress during the whole shearing process; the “servo control” provides the ability to control the translational velocity of the target wall using a servomechanism in order to apply or maintain a desired force [51]. The equation for the servo control is shown as follows:where is the change of contact stress of the target wall, is the average normal stress of the target wall, is the number of contacts between particles and the target wall, is the velocity of the target wall, A is the area of interface, and is the change of time. The shearing processes in simulation were controlled by displacement, just like those in laboratory tests. The top wall at first moved in the positive X direction, i.e., moving rightwards until reaching 50 mm displacement and then moved in the opposite direction (negative X direction) for 100 mm, namely, reaching −50 mm, followed by a return displacement to its original position, and this was the first shear cycle. Subsequently, the top wall repeated the shear cycle for another three times.
As mentioned above, there were two shapes of agglomerates in the numerical sample. The definition of “composite pattern” is the proportion of these two shapes to constitute a numerical sample. A set of contrast simulations with different composite patterns (i.e., SY_2_0, SY_2_15, SY_2_37.5, and SY_2_60) and with unbreakable agglomerates (i.e., SN_2_15, SN_4_15, and SN_7_15) were performed to study the effect of composite pattern on the evolution of breakage and microstructure and the effect of agglomerate breakage on microstructural evolution, respectively. When modelling unbreakable agglomerates, a very large bond strength was applied so that the bond between spheres within agglomerates could not break. The other parameters and testing procedure were kept the same with those for breakable agglomerates. There are nine simulations with different simulation conditions that are summarized in Table 1. We divide these nine simulations into three simulation groups. (i) The reference group includes those simulations with breakable agglomerates under various constant normal stress conditions, i.e., SY_2_15, SY_4_15, and SY_7_15. (ii) The composite pattern group involves simulations with breakable agglomerates with different composite patterns under one constant normal stress of 200 kPa, i.e., SY_2_0, SY_2_15, SY_2_37.5, and SY_2_60, as seen in Figure 5. (iii) The unbreakable group involves those simulations with unbreakable agglomerates under various constant normal stress conditions, i.e., SN_2_15, SN_4_15, and SN_7_15.
The linear model, where a linear elastic frictional behavior was assumed, was applied to the contacts between agglomerates. The linear contact stiffness model was defined by the normal and shear stiffness, and a no-tension mode was adopted in the linear model. The normal stiffness and shear stiffness of walls were one order of magnitude larger than those of spheres. To simulate flaws in naturally brittle materials, i.e., angular gravel in this paper, there are three main methods. Some researchers (e.g., [18, 19, 52]) adopted the method that a random percentage between 20% and 50% of balls were removed and then, each agglomerate was given a random rotation, to mimic statistical strength variation. The second method is that Zhou and Song [10] set random virtual cracks with their lengths obeying Gaussian distribution. The last approach, used by many researchers (e.g., [20, 22, 53, 54]), was that spheres within agglomerates were connected by parallel bonds or contact bonds that were given the shear and normal strength also complying with Gaussian distribution. In this paper, the third method was adopted, and the linear parallel bond model was applied to contacts between spheres within agglomerates. The magnitudes of the normal and shear contact stresses were limited by the shear strength and the normal strength, respectively, both obeying Gaussian distribution so that a random strength of agglomerates could be obtained to characterize the strength variation. If the shear strength limit or normal strength limit is exceeded, the bond is thus broken.
Referring to previous studies [20, 22, 48, 55], the linear contact stiffness fell into the order of magnitude of ×106 N/m and bond stiffness was set as a higher order of magnitude [17, 22], of which ×109 Pa/m in this paper. Microparameters (i.e., sphere stiffnesses, bond stiffnesses, and frictional coefficient as well as bond strengths in this paper) were iteratively adjusted until macromechanical results of numerical simulations closely agreed with those from laboratory tests, as suggested in other literature (e.g., [20, 22, 45, 46, 55, 56]). Excellent agreement between numerical and experimental results for simulations under five normal stress conditions could be seen in Figure 6, meaning that the material properties in Table 2 are verified.
(a)
(b)
(c)
(d)
(e)
2.3. Strength Variation
The tensile strength of the particles was defined by Jaeger [57] aswhere F is the greatest force and d is the initial separation of the platens in the single-particle crushing test.
Single grain fracture experiments were simulated with 30 shape-A agglomerates and 30 shape-B agglomerates of the same size, while the bond strengths were assigned randomly, obeying Gaussian distribution mentioned above. The survival probability is calculated as the mean rank value [18]:where i is the rank position of grain when sorted in an ascending order of peak stresses and N = 30.
Weibull [58] proposed a function that could be used to describe the survival probability of apparently identical test pieces in tension tests.where is defined as the shear stress at which 1/e ≈ 37% of the grains survive and m is the Weibull modulus.
Figure 7 shows that the survival probability curve fits the Weibull distribution well and the Weibull modulus, i.e., m, of shape-A agglomerate is smaller than that of shape-B agglomerate. m quantifies the scatter of strength values, and high values of m indicate low scatter [10]. As observed from Figure 7, of shape-A agglomerate is up to nearly 5.0, much larger than that of shape-B agglomerate, and shape-A agglomerate has a higher degree of scatter in terms of the Weibull distribution.
(a)
(b)
3. Macromechanical Response
Apart from the simulations in the reference group as mentioned above (i.e., SY_2_15, SY_4_15, and SY_7_15), two additional simulations with identical simulation conditions except under different normal stress conditions were performed for further verification (Figures 6(a) and 6(e)). Figure 6 shows that the macromechanical results of DEM simulations under five constant normal stress conditions, compared with the laboratory test results from Zhang’ s thesis [4], provided the plots of shear stress versus horizontal displacement in the first and tenth shear cycles (i.e., Cycles 1 and 10) and stated these plots within these two cycles were basically uniform. In Cycle 1, there are remarkable strain softening phenomena under these three normal stress conditions in DEM simulations, and from Cycle 2 onwards, it is interesting to note that DEM simulation results show similar trends to experimental ones, although shear stress is relatively much lower in DEM simulations.
Observed from these DEM simulations in four composite patterns (Figure 8), the peak stresses are almost the same, reflecting that composite pattern has a limited effect on macromechanical responses. Figure 9 shows that simulations in the unbreakable group normally have a higher degree of shear stress compared with those in the reference group. This is more evident in high normal stress conditions (SN_4_15 and SN_7_15), where their shear stresses are approximately twice the magnitudes of shear stresses of those simulations with breakable agglomerates, i.e., SY_4_15 and SY_7_15, which demonstrates that intact agglomerates could provide a higher degree of shear resistance, especially in high normal stress conditions.
(a)
(b)
(c)
4. Evolution of Breakage
4.1. Evolution of GSD
It is important to quantify the breakage of particle under different loads, and much work has been done in this study. Lee and Farhoomand [59] proposed to use the control grain size D15 to characterize the index B15 (the ratio of D15i/D15a, with D15i and D15a measured before and after crushing, respectively), to describe the breakage extent of coarse particles. Lade et al. [60] proposed a similar index B10 to characterize this breakage extent. Miura and O-Hara [61] found that there is a good functional relationship between the newborn surface area and the plastic work, and thus they proposed to use the newborn surface area as a characteristic value to judge the breakage extent. And, there are some other characteristic indices, like Sun and Wang [62] used acoustic emission to characterize calibrated sand breakage extent. However, some literature connected breakage indices with the evolution of GSD. Marsal [63] proposed the use of Bg, representing the maximum increase in percentage passing sieve. Despite its convenience of computing, Bg has some disadvantages such that it is easily influenced by particle size and aggregate gradation and hard to present the realistic change in different size fractions. Hardin [64] generalized former studies on breakage mechanics and then proposed breakage potential Bp, total breakage Bt, and relative breakage index Br, which could reflect the total breakage extent. Einav [65] absorbed the fractal theory and developed a continuum theory to modify Harbin’s breakage measurement. He then proposed that there was an ultimate GSD for the ultimate particle breakage state and the fractal dimension of 2.5–2.6 was useful for most cases in simulating granular crushable materials; in this paper, it is set as 2.6. The modified Harbin’s breakage index by Einav [65] is calculated as follows:where and are the minimum and maximum particle sizes, respectively. and are the function of initial GSD and current GSD, respectively, and the function of ultimate GSD is given by
Figure 10 exhibits that the calculated GSD curves at the end of Cycles 2 and 4 not only demonstrate the variation in the sizes of particles in the numerical sample but also reflect the evolution of breakage. The final GSD was derived from the GSD of Cycle 100 of the laboratory test under a constant normal stress of 200 kPa [4], of which grain size (i.e., abscissa value) was amplified by five times according to the similarity law. As seen in Figure 10, the trends of GSD of DEM simulations are very similar to the final GSD. As the normal stress increases, the shear stress rises correspondingly and more agglomerates are crushed into finer agglomerates or independent particles (i.e., floating spheres independent of any agglomerate), indicating that the normal stress has a non-negligible effect on the evolution of breakage (Figure 10(a)). In contrast, composite pattern holds a much smaller impact on this evolution (Figure 10(b)). With the rise in the fraction of shape-A agglomerates in the sample (from 0% to 60%), only a very little bit larger degree of breakage could be observed from the variation in Figure 10(b) and the GSD curves of the four simulations in the composite pattern group are almost uniform.
(a)
(b)
To quantify the breakage evolution, Bg and BrE, as mentioned above, were calculated in this paper. As Figure 11 shows that Bg of all simulations are higher than their BrE and the two breakage indices correlate well with these observations from GSD in Figure 10. The increase rates of Bg and BrE decrease for almost all simulations except for SY_2_37.5 when calculating Bg, its increase rate nearly keeping constant. It indicates that the variation in the early breakage extent of the sample is higher than that in the later one. In comparison of the composite pattern group, Figure 11 is clearer and more direct in indicating that the degree of breakage increases a little with the increase of shape-A agglomerates and shape-A agglomerates are easier to crush.
4.2. Agglomerate Breakage
“Particle” in this paper includes two kinds: agglomerate with at least two spheres and independent sphere (i.e., independent of its original agglomerate). Two types of breakage modes, defined by Tapias et al. [27], are introduced here: splitting breakage and comminution. As for splitting breakage, for example, one original agglomerate with 14 spheres is broken into two particles, one with 8 spheres and one with 6 spheres, which is marked as 14 = 8 + 6. In the “comminution” type, a new independent particle detaches from the original agglomerate and becomes a floating sphere in the numerical sample due to external force, i.e., 14 = 1 + 13 for shape-A agglomerate and 12 = 1 + 11 for shape-B agglomerate. There is one special probability: under a continuous comminution process, taking shape-A agglomerate as an example, the original agglomerate may end in 14 independent spheres.
A close examination in a statistical sense on the evolution of the number of particles is carried out. Figures 12 and 13 present the numbers of independent spheres and agglomerates, respectively, for simulations with breakable agglomerates, and the normalized numbers of independent spheres and agglomerates are presented to help illustrate the variation more clearly and directly. Figure 12 shows that the increase rates of the numbers of independent spheres of these simulations gradually slow down, signifying that the frequency of comminution decreases. It is evident that as the normal stress increases, increases correspondingly. However, in comparison of the composite pattern group, their normalized numbers of independent spheres are almost uniform, demonstrating that composite pattern does not affect the frequency of comminution. The numbers of agglomerates of the four simulations in the composite pattern group increase continuously within the whole four shear cycles (Figure 13). This is because original agglomerates are crushed into finer ones due to splitting breakage. The agglomerate increase process, for example, evolves as follows: 14 = 8 + 6; 8 = 5 + 3; 6 = 3 + 3; it shows that how one original agglomerate becomes four finer agglomerates. In comparison of the composite pattern group, apart from SY_2_37.5, the normalized numbers of agglomerates for SY_2_0, SY_2_15, and SY_2_60 increase with the growth of the fraction of shape-A agglomerates, probably demonstrating that there is a higher probability of splitting breakage with more shape-A agglomerates. As for SY_7_15, there are two stages. In the first stage, its increase rate of is higher than other simulations within Cycle 1, and reaches its peak, from 4257 to 4622, slightly after the end of Cycle 1, when many agglomerates are crushed into finer ones due to splitting breakage. In the second stage, this figure turns to decrease continuously to 4398 finally and many finer agglomerates are crushed into independent spheres, i.e., 3 = 2 + 1; 2 = 1 + 1, which is proved by Figure 12, which indicates that comminution plays an important role in this stage. for SY_4_15 turns to become stable nearly at the end of Cycle 2 after a steady rise, during which process of splitting breakage and comminution jointly work to keep in dynamic equilibrium, but it decreases slightly to 4627 at the end of Cycle 4.
(a)
(b)
(a)
(b)
As observed from Figure 14(a), finer particles mainly appeared on the interface and the area nearby, which indicates that breakage predominantly occurred on and near the interface in this laboratory test; this phenomenon can also be found in numerical simulations (Figure 14(b)). It is thus interesting to visualize the particles of simulations with breakable agglomerates on the interface in the final state, as seen in Figure 15. Notable increase in independent spheres (red) and 2-agglomerates (yellow) is found when comparing Figures 15(a)–15(c), and red color almost occupies the whole of the interface of SY_7_15 finally. This indicates that the breakage extent increases with the rise in normal stress, which could be observed from Figures 10(a) and 11. For Figures 15(d)–15(f), their color composites seem rather similar although several more 14-agglomerates (purple) and fewer 12-agglomerates (royal blue) exist in the final states of SY_2_37.5 and SY_2_60, compared with the final state of SY_2_15 due to a higher fraction of shape-A agglomerates initially, while dozens more 12-agglomerates (royal blue) exit in the final state of SY_2_0. It seems that composite pattern influences the breakage extent vanishingly, which is also supported by the observations from Figures 10(b) and 11.
(a)
(b)
(a)
(b)
(c)
(d)
(e)
(f)
5. Microstructural Insights of Crushable Angular Gravel in Simulations
5.1. Evolution of Anisotropies
The anisotropy concept was introduced here to reflect the fabric composite and microstructural variation of crushable angular gravel in the sample. There are two kinds of anisotropy sources: geometrical anisotropy, which is contact normal anisotropy, and mechanical anisotropy, including normal force anisotropy and tangential force anisotropy. According to the relevant literature mentioned in Section 1, a summary of the relevant anisotropy tensors and parameters used to quantify the degree of anisotropy is given as follows.
The contact normal anisotropy tensor is expressed aswhere is the deviatoric part of , which is given by
Note that is the contact normal (also is the unit vector of the contact branch vector for spheres) and is the total number of contacts. The normal force anisotropy tensor is defined aswhere is the deviatoric part of and is the average normal force derived from . is calculated aswhere is the magnitude of the contact normal force.
Similarly, the tangential force anisotropy tensor iswhere is the deviatoric part of and is expressed as follows:
Note that is the unit vector of the direction of the contact tangential force and is its magnitude. The three deviatoric tensors of anisotropy sources are used to quantify the degree of anisotropy, as can be conveniently expressed by using , , and :where , , and are contact normal anisotropy, normal force anisotropy, and tangential force anisotropy, respectively. Figures 16–18 present the evolution of three anisotropy invariants , , and for all simulations.
(a)
(b)
(c)
(a)
(b)
(c)
(a)
(b)
(c)
5.1.1. Contact Normal Anisotropy ()
Figure 16 shows the following observations. (i) Every time the top wall changes the moving direction, i.e., when N = 0.25, 0.75, 1.25, 1.75, 2.25, 2.75, 3.25, or 3.75, in all simulations will experience a “phase,” during which will slump to a much lower level and then rise back to its peak value again; such phase is called “adjustment phase,” which will be analyzed and discussed later in Section 5.3. (ii) of SY_2_15 and SY_2_37.5 along the negative X direction is lower than that along the positive X direction while this phenomenon does not exist in that of SY_2_0 and SY_2_60. It is also interesting to note that, in the composite pattern group, of SY_2_37.5 has the largest value when shearing is along the positive X direction and the smallest value along the opposite direction. (iii) Simulations in the unbreakable group normally have higher degrees of than those simulations in the reference group, and during the shearing process, simulations in the reference group witness a decrease in peak values. For example, of SY_2_15 reaches its peak at 0.281 initially within Cycle 1, and it finally reaches its peak at 0.221 within Cycle 4. This observation is more noticeable in higher normal stress conditions. It seems that finer agglomerates and independent spheres lead to a lower degree of .
5.1.2. Normal Force Anisotropy ()
As seen in Figure 17, the evolution trend of is fairly similar to , which was also observed in other simulations (e.g., [11, 34, 35]), although there are still some differences between and . (i) In the composite pattern group, unlike of SY_2_37.5 is not always the largest value along the positive X-direction and not all the way the smallest one along the negative X-direction as well. of SY_2_60 is the largest at the beginning of Cycles 1 and 3 and the smallest within Cycles 2, 3, and 4 when shearing is along the negative X direction. (ii) of SY_2_15 along the positive X direction is not bigger than that along the opposite direction.
5.1.3. Tangential Force Anisotropy ()
Figure 18 presents a quite distinct evolution trend in contrast with and and the following characteristics. (i) There is no “adjustment phase” like and in the evolution of , but the alteration of shearing direction causes an increase in most cases. (ii) In the reference group, of both SY_7_15 and SY_4_15 shows a noticeable decreasing trend, but the former is quicker to reach a stable value (when N = 1.5) though within a moderate fluctuation range. In contrast, simulations in the unbreakable group have no such decreasing trend, but their fluctuates around a larger value within a bigger range compared with simulations in the reference group. (iii) Normal stress has a much more evident impact on than and ; increases with the rise in normal stress despite that of SY_4_15 is a little bit larger than that of SY_2_15 in the end. (iv) The degree of decreases with the growth of fraction of shape-A agglomerates, indicating that composite pattern has a big effect on . (v) of simulations in the unbreakable group is close to that of those simulations in the reference group initially while with the shearing process continuing, for the latter simulations is much lower than that for the former ones, and it is more marked in high normal stress conditions. This observation agrees well with and .
5.2. Relationship between Anisotropies and Shear Strength
At first, we follow the definition of stress tensor proposed by Christoffersen et al. [66] aswhere V is the total volume of the assembly, is the total number of contacts, is the contact force vector, and is the branch vector connecting the centers of two contacting particles. The mean effective stress and deviatoric stress q can be expressed as follows:where is the deviatoric part of the stress tensor .
Rothenburg and Bathurst [32] performed biaxial tests and proposed an analytical relationship connecting deviatoric stress and anisotropy invariants:where and are the principal stresses.
Rothenburg and Bathurst [32] further derived from equation (16) using the mobilized angle of friction for cohesionless soil into
Guo and Zhao [35] modified the relationship that was proposed by Chantawarungal [67], and it could be given by
Equation (18) correlated well with simulations results in some researches (e.g., [33, 35, 36]). But, in this paper, there are two observations distinct from DEM simulation results in those researchers’ literature. (i) is much bigger than other anisotropy invariants , , and , which is different from those simulation results [33, 35, 36]. (ii) There is an obvious “adjustment phase” in , , and , but it is not found in relevant literature. This is probably why equation (18) cannot be applied to the research in this paper. However, through the polynomial fitting method, the following relationship is given by
Figure 19 shows good agreement between the fitting relationship, i.e., equation (19), and DEM simulation results. Observed from equation (19), the contribution from the contact normal anisotropy is the least to the deviatoric stress , while normally provides the largest part, but within the “adjustment phase,” makes a larger contribution than , as seen in Figures 17 and 18. almost has no fluctuation because of the constant normal stress (Figure 20), so the fluctuation of is mainly attributed to , which is regarded as the generalized shear stress in soil mechanics; this fluctuation is caused by the interface shearing behavior. The evolution of has a similar trend with shear stress of the interface. and grow together initially, and once the shearing direction changes, both of them decrease to a much lower value. But, will return to its peak with the same sign (always positive), whereas will reach a peak with an opposite sign. It is interesting to find that both and will reach a peak near the 30 mm or −30 mm displacement after the alteration of shearing direction almost simultaneously. Moreover, in comparison of the composite pattern group (Figure 20(b)), the magnitude of increases with the decrease in the fraction of shape-A agglomerates, signifying that shape-B agglomerates are more compressive than shape-A agglomerates.
(a)
(b)
(c)
(a)
(b)
(c)
5.3. Adjustment Phase and Evolution of Contact Orientation Distribution
Figure 21 is straightforward and clear in displaying the “adjustment phase” as mentioned above. It is obvious to see that there is a lagging effect that and are slower than to reach the lowest point. The interface shearing behavior only affects the distribution of contact orientation on the X-Z plane and therefore that on the X-Y and Y-Z plane are not presented here. Figure 22(a) displays the distribution of contact orientation on the X-Z plane in point A that denotes the lowest degrees of , , and , where its major axis (the red dotted line) is nearly parallel to the Z direction and very closely passes through the circular center. However, in point B (Figure 22(b)), where the degrees of , , and are up to 0.248, 1.426, and 0.274, respectively, the major axis nearly passes along the 120°–315° line, which is far away from the circular center. Comparing points A and B, we can find that the included angles between their major axes and the 90°–270° line reflect the change of contact orientations induced by shearing behavior, and the distances between their major axes and the circular center signify the magnitudes of , and . The angle and the distance specifically denote the included angle between the major axis and the 90°–270° line and the length between the major axis and the circular center, respectively, thereinafter.
(a)
(b)
(c)
(a)
(b)
(c)
(d)
(e)
(f)
(g)
When shearing at point C (Figure 22(c)), in which , , and return to a relatively lower degree, the distance narrows, and the distribution of contact orientation is closer to its original state (i.e., the state in point A), compared with point B. But, with the shearing process along the negative X direction continuing, the distance and the angle become wider in point D (Figure 22(d)) and point E (Figure 22(e)), where , , and come back to a higher level. The evolution of contact orientation distribution from point A to point E demonstrates that such evolution itself is the forming reason of “adjustment phase.”
When comparing point C, point F (Figure 22(f)), and point G (Figure 22(g)), the distance in point F is the smallest among these three points and that in point G ranks second. Meanwhile, , , and in point F are the smallest and those in point G are the second smallest correspondingly. It validates that the distance between the major axis and the circular center is able to be used to compare the magnitudes of , , and once again.
5.4. Evolution of Coordination Numbers
is the total number of contacts between particles in the numerical sample. of simulations with unbreakable agglomerates, i.e., simulations in the unbreakable group, demonstrate a consistent decreasing trend (Figure 23). In contrast, simulations with breakable ones, under high normal stress conditions, i.e., SY_4_15 and SY_7_15, show a total opposite trend, increasing from 34202 to 35869 and from 38080 to 42371, respectively. It demonstrates that the generation of independent spheres and finer agglomerates could cause more contacts. However, we can see other simulations with breakable agglomerates, i.e., SY_2_15, SY_2_37.5, and SY_2_60 nearly remain stable, whereas SY_2_0 sees a slight decrease. The coordination number is defined as the number of contacts per particle, and the coordination number for particles (i.e., agglomerates plus independent spheres) is given bywhere is the number of contacts between particles around one particle and is the number of particles in the numerical sample. As Figure 24 presents, of simulations with unbreakable agglomerates decreases much slower and meanwhile, their final values are larger than those simulations with breakable agglomerates. It is also interesting to note that the decrease rates of of simulations with breakable agglomerates become identical after Cycle 1, although they are totally different initially. In comparison of simulations with breakable agglomerates, of SY_7_15 and SY_4_15 is higher, but their is smaller than that of simulations under 200 kPa normal stress, i.e., SY_2_0, SY_2_15, SY_2_37.5, and SY_2_60, due to the existence of too many independent spheres, which increase the contacts but at the same time, reduce the value of . Besides, the magnitude of decreases with the reduction in the fraction of shape-A agglomerates, signifying that its magnitude is determined by composite pattern.
To investigate a coordination number related to breakage, the coordination number for parallel bonds is used here:where stands for the number of parallel bond contacts between spheres within agglomerates around one sphere and stands for the number of all spheres in the numerical sample. Figure 25 shows the evolution of . When comparing simulations in the reference group, of SY_7_15 is always the smallest one while that of SY_2_15 is the biggest, and meanwhile SY_7_15 has the largest degree of breakage while SY_2_15 has the smallest, as seen in Figures 10 and 11. This agrees well with Mcdowell and Bolton’s conclusion [40] that the probability of breakage reduces with an increase in the coordination number. It is also interesting to note that their is very close to each other at the end of Cycle 4, and it seems that there is an ultimate particle breakage state existing that the same sample will reach an identical particle breakage state, whatever their normal stresses are, as proposed by Einav [65]. Comparison of the composite pattern group is still also consistent with their conclusion [40] that SY_2_60 sees the smallest value in and the highest breakage extent correspondingly while and the breakage extent of SY_2_37.5 rank third and second, respectively; of SY_2_0 is the largest among them, and its breakage extent is the smallest one. Figure 25 also shows that normal stress strongly influences the trend of , whereas composite pattern does affect its trend vanishingly. We can conclude that normal stress deeply affects the trend of , while composite pattern determines its magnitude.
6. Conclusions
The microbehavior of crushable angular gravel under the cyclic soil-structure interface test was investigated using DEM. The crushable angular gravels were modelled by filling two scanned angular gravels with spheres connected by parallel bonds. The similar trend of shear stress-shear displacement of DEM simulations to experimental results has been observed. Intact agglomerates could provide a higher degree of shear resistance from the comparison between DEM simulation results with breakable agglomerates and with unbreakable agglomerates. Novel findings from the study in terms of the evolution of breakage and microstructure are summarized as follows:(1)Breakage mainly occurs on the interface and the space nearby, which is deeply influenced by the interface shearing behavior, and there is a strong positive correlation between normal stress and the breakage extent. Higher normal stress leads to the generation of a larger number of finer particles, i.e., finer agglomerates and independent spheres. In contrast, composite pattern has a minor effect on the breakage extent, which increases a little bit with the increase of shape-A agglomerates. This indicates that shape-A agglomerates are easier to crush.(2)Composite pattern determines the magnitudes of and but has no clear and specific effect on , , and . It is worth noting that shape-B agglomerates with a higher level of sphericity are more compressive. In addition, there is a noticeable positive correlation between normal stress and the magnitudes of , , , and . Breakage noticeably reduces the degree of anisotropy invariants, i.e., , , and , because intact agglomerates have an intrinsic higher degree of anisotropies, whereas finer agglomerates and independent spheres have a lower degree of anisotropies.(3)A new stress-force-fabric relationship through polynomial fitting could be able to describe the relationship between anisotropies and shear strength for all numerical simulations in this study. The two mechanical anisotropy invariants play a predominant role in providing shear resistance, in which makes a much larger contribution.(4)The “adjustment phase” could be evidently found in the evolution of , , and , and they bear much more effect from the interface shearing behavior than . The forming cause of “adjustment phase” is the evolution of contact orientation distribution on the X-Z plane induced by the interface shearing behavior.(5)The interface shearing behavior also causes the evolution of the coordination numbers for particles and for parallel bonds . Both the two coordination numbers decrease during the whole shearing process, and they are significantly influenced by the magnitude of normal stress. Breakage increases the number of contacts but reduces the magnitude of because of the generation of a larger number of finer agglomerates and independent spheres. Composite pattern determines the magnitudes of these two coordination numbers, and normal stress deeply influences the trend of . In addition, , also related to breakage, under different normal stress conditions, is very close to each other at the end of simulations, which demonstrates that there is an ultimate particle breakage state exiting that the same numerical sample will reach an identical particle breakage state, no matter what their normal stresses are.
Data Availability
The data used to support the findings of this study are available from the corresponding author upon request.
Conflicts of Interest
The authors declare that they have no conflicts of interest.
Acknowledgments
This work was financially supported by the National Key R&D Program of China (2017YFC0404905) and the National Natural Science Foundation of China (Grant No. 51679264).