Multi-scale investigation of active failure for various modes of wall movement

Abstract

Retained backfill response to wall movement depends on factors that range from boundary conditions to the geometrical characteristic of individual particles. Hence, mechanical understanding of the problem warrants multi-scale analyses that investigate reciprocal relationships between macro and micro effects. Accordingly, this study attempts a multi-scale examination of failure evolution in cohesionless backfills. Therefore, the transition of retained backfills from at-rest condition to the active state is modeled using the discrete element method (DEM). DEM allows conducting virtual experiments, with which the variation of particle and boundary properties is straightforward. Hence, various modes of wall movement (translation and rotation) toward the active state are modeled using two different backfills with distinct particle shapes (spherical and nonspherical) under varying surcharge. For each model, cumulative rotations of single particles are tracked, and the results are used to analyze the evolution of shear bands and their geometric characteristics. Moreover, dependencies of lateral pressure coefficients and coordination numbers, as respective macro and micro behavior indicators, on particle shape, boundary conditions, and surcharge levels are investigated. Additionally, contact force networks are visually determined, and their influences on pressure distribution and deformation mechanisms are discussed with reference to the associated modes of wall movement and particle shapes.

Graphical abstract

Keywords

discrete-element modelling / granular materials / retaining walls / particle shape / arching

1 Introduction

The design of all civil engineering structures requires the consideration of soil-structure interaction. However, the soil-structure interaction problem geotechnical engineers face when designing retaining structures is more complex and interesting since macro, meso, and micro effects are all prominent on response and intertwined. Therefore, understanding the mechanics of this problem necessitates the consideration of behavior at all scales. Macroscale effects include the properties of the structural system, boundary conditions, and loading type, whereas microscale effects involve grain shape, grain size, and grading. On the other hand, micro effects influence macro response through mesoscale mechanisms that control the arching mechanism and define the sizes, shapes, and locations of shear bands and the formations that control bulk behavior.

Owing to the complexity of the problem and the fact that retaining structures are one of the most common civil engineering structures, the behavior and response of retained backfills attract the interest of engineers and researchers. The mechanisms responsible for the loads and stresses that backfills impart on the retaining structures have attracted attention because their results directly affect the design. That is why, even though efforts to understand backfill behavior dates back to the works of Refs. [1,2], the interest in the topic never reduced [3–10]. Recently, attempts to examine retained backfill behavior peaked due to the advances in sensing and monitoring technology and computational power. Notably, the adaptation of particle image velocimetry to geotechnical engineering [11–14], use of photo-elastic grains [15,16], and the introduction of small-sized high-quality sensors [13,14,17] enhanced the amount and quality of information we can glean from physical model studies. This quality of the data allowed the examination of shear band formation mechanisms [18–23] and measuring actual magnitudes and distributions of stresses acting on retaining structures [14,17,24]. Most of these works focused on simpler boundary conditions where the retaining structures’ mode of displacement is rigid lateral translation [17,24–28]. At the same time, several studies used complex mechanisms to investigate the influence of different wall displacements and deformations [13,14,19,29] on backfill response. Nevertheless, the examination of micro properties such as grain shape and meso responses as force chains for different boundary conditions is still a daunting task. At this point, the Discrete Element Method (DEM) [30] provides a powerful medium to investigate backfill behavior at multiple scales. DEM is a particle scale numerical method that considers particles as a system of distinct interacting bodies. DEM makes it possible to observe the interactions among grains under different boundary conditions. It enables simulation of individual grains, which makes it possible to have a better insight into micro-scale. A distinct merit of DEM is the possibility of running virtual experiments without any limitations in changing boundary conditions and particle shape. Additionally, DEM enhances observation of force distribution among discrete particles, and thus it is possible to study complex mechanisms such as arching and deformation localization.

Recently, several numerical methods have been developed for fracture simulation, of which the cracking particle method (CPM) is shown to be one of the reliable approaches [31–33]. Additionally, a dual-horizon peridynamics (DH-PD) method is developed by Ren et al. [34,35], which enhances the study of crack propagation in composite materials. These robust methodologies might be useful for studying the shear failure in cohesive soils mixed with rubbers. However, in this study DEM method is preferred since DEM is capable of simulating individual grains, which enhance the study of different complex mechanism at particle scale. It is already shown that particle interlocking is a direct function of shape factor, which can be easily varied within DEM. Additionally, DEM enables the study of anisotropy of contact force distribution, and this can be utilized to justify the complex failure phenomena such as arching in particulate assemblies. Moreover, the displacement and rotation of individual particles can be tracked within DEM, which provides the precise location of deformation localization in a granular system.

DEM is a valuable method to investigate soil behavior and soil-structure interaction [36–38]. Recently, many researchers utilized DEM to study retained backfill behavior. They examined the influence of the mode of wall displacement [39–41], the effect of boundary conditions [21,42], formation of force chains and arching [23,27,42–46], the impact of grain shape and size [45,47] on backfill response and the distribution of lateral pressures. All these works confirm that DEM is beneficial for observing and understanding the backfill response at different scales.

Accordingly, this study uses DEM to study the response of a backfill retained behind a wall yielding toward active failure at multiple scales. With this work, influences of particle shape, surcharge, and mode of wall movement on the active failure mechanism are considered. For this purpose, backfills with identical initial porosities are prepared with either spherical or sub-rounded particle shapes. These models are subjected to various modes of wall movement: horizontal translation, rotation about the top, and rotation about the bottom. Additionally, to understand the influence of surcharge, effective vertical stress acting on the backfill is varied between models. Results are analyzed to identify shear bands, calculate lateral earth pressure coefficients, and determine lateral force distribution along the wall height. Furthermore, to understand particle level mechanisms and their influences on macroscale response, contact force networks, arch formations, and coordination numbers are extracted for all samples. Finally, based on the available data, interactions between mechanisms at different scales and their effect on the backfill response are discussed.

2 Methodology

2.1 Working principles of discrete element modeling (DEM)

The initial step in DEM is to generate the geometry of the system, where particles can overlap while maintaining inter-particle contact. The fundamental idea of DEM is based on two different physical laws: Newton’s Second Law applied to the particles, and Force-displacement Law at the contacts. DEM performs repeated calculations using basic physical laws. First, the contact forces of the particles are calculated from their displacements using the Force-displacement Law. Considering both the contact forces and body forces acting on the particle, the total normal and shear forces acting on each particle are separately determined at each time step. Following the calculation of the contact forces, Newton’s Second Law of motion is applied to every individual particle. Between two succeeding time-steps, the accelerations of all particles are calculated based on particle masses and the forces acting. Particles’ new velocities and displacements are determined with double integration of acceleration over each time step.

For particle i, the equation of motion is given by Eq. (1):

(1)m_i\left({\displaystyle \frac{{\mathrm{d}}^{\mathrm{2}}{x}_i}{\mathrm{d}{t}^{\mathrm{2}}}}\right)=f_i+m_{i}g,

where m_i is the mass of the particle, t is time, x_i is its position, g is the acceleration due to gravity and f_i is the force acting on the particle due to particle contacts as defined by Eq. (2):

(2)f_i=\sum f_i^{\mathrm{c}}.

The rotational motion equation for particle i is calculated using Eq. (3):

(3)I_i\left({\displaystyle \frac{\mathrm{d}{\omega }_i}{\mathrm{d}t}}\right)=T_i.

In Eq. (3), I_i is the moment of inertia for particle i, ω_i is its angular velocity and T_i is the total torque acting, which is defined by Eq. (4) where l_i and r_i are the branch vector and radius of particle i (r_j is the radius of particle j), defined by Eq. (5):

(4)T_i=\sum l_i^{\mathrm{c}}\ast f_i^{\mathrm{c}},

(5)I_i=r_i-r_j.

The principles of DEM theory for calculating the positions and forces of all individual particles are explained in detail by O’Sullivan [48]. The summary of the calculation cycle of DEM is given in Fig. 1.

2.2 The numerical model

The numerical models are prepared using a three-dimensional DEM set up using EDEM® 2018 [49] bulk material simulation software. In these models, structural elements (i.e., model retaining wall and model frame) are rigid, and deflection is not allowed. Such restrictions on the modeling properties are required because the primary purpose of this study is to investigate the response of retained granular assemblies to the development of active conditions for different boundary conditions.

Based on the particulate nature of soils, there are various methodologies for multi-scale analyses depending on whether granular or continuum approaches are assumed. To investigate the fracture phenomena in clay nanocomposites, Talebi et al. [50] have utilized semi-concurrent and concurrent multi-scale methods to study fracture mechanisms, which couple two continuum approaches or atomistic to continuum domains, respectively. Besides, Budarapu et al. [51] reveal an adaptive atomistic (fine-scale)-continuum (coarse-scale) numerical method for quasi-static crack growth. The current study prefers a granular approach and adopts DEM to provide a multi-scale overview of the active failure state. This is achieved by varying the particle scale properties (e.g., particle shape) and boundary conditions (e.g., mode of wall movement and surcharge load). This way, the variations in contact force distributions, coordination numbers (micro-scale), -temporally averaged- normalized lateral stress distribution (macro-scale) are considered.

The models utilize periodic boundaries along y axis as seen in Fig. 2, resulting in a state appropriate for predicting plane strain conditions. When periodic boundary condition is defined, the program surrounds the periodic cell with identical copies [48]. The width of the model in this study between periodic boundaries equals 20 mm (Fig. 2). Retained backfill is horizontal and level with the top of the wall. Depending on the model, the wall either horizontally translates or rotates toward active failure. The rotation can be about the top of the wall or the bottom of the wall. Moreover, the adverse effect of the bottom boundary is eliminated by defining a layer of particles with a depth equal to the height of the wall below the bottom of the wall (Fig. 2). The particles used in this layer are identical to the particles used behind the wall. Additionally, a zone referred to as the optimum area of interest (AoI) is defined (Fig. 2). The extent of AoI is selected such that all possible deformations for all models occur within this zone. This will make it possible to focus on the crucial mechanisms when the results of the models are presented in the following sections. These mechanisms include shear bands and contact network forces.

The material parameters for DEM simulations of silica sand are directly obtained from de Bono and McDowell [52] (Table 1), which are in a similar range as reported by Refs. [53–56]. de Bono and McDowell [57] also compared and validated stress-strain curves from DEM simulations with the experimental triaxial test data from Yammuro and Lade [58]. Also, de Bono and McDowell [56] simulated isotropic compression test for silica sand in DEM. They showed a good comparison with the experimental results in the study of Bandini and Coop [59]. Considering all the above studies, one can conclude that the material parameters used in this study have already been calibrated. Additionally, for validation purposes (in Section 2.4), the dimensionless geometry of failure surfaces from DEM simulations are obtained in the current study and compared with the physical model test results provided by Altunbas et al. [26].

The model backfills are prepared using either spherical particles or using overlapping spheres (clumps). Hereafter, single spheres and clumps are referred to as M1 and M2, respectively. M2 consists of four sub-spheres and represents a sub-rounded particle. The shapes of M1 and M2 particles are classified in terms of sphericity and roundness as suggested by Cho et al. [60]; see Table 2 for shape characteristics. In the remainder of this paper, backfills will be referred to according to the type of the particles they are composed of; either as M1 or M2. Both M1 and M2 backfills are poly-disperse granular systems, which are classified as poorly-graded according to USCS (see Table 3 for particle size distribution details).

Backfills are prepared by the air-pluviation method for all models. For this purpose, a dynamic factory is used to ensure a constant dropping distance for all particles during the filling process. All models are prepared at a typical target porosity which is equal to 0.35 (Table 3). This porosity is achieved by adjusting the friction coefficient (μ_s). Accordingly, μ_s = 0 for M1 and μ_s = 0.17 for M2 models during air-pluviation.

A schematic representation for various modes of wall movements is presented in Fig. 3. For simulations with horizontal wall translation (HT), a constant speed of 1 mm/s is considered. For models with wall rotation, either about the top (RT) or the base (RB), a constant rotational speed of 0.0167 rad/s is assigned. This low speed ensures the quasi-static flow of the particulate system [61]. Problems such as direct shear and silo discharge also have similar shearing regimes and are considered quasi-static problems [62]. The data obtained from the analysis is saved with a frequency of 100 Hz.

To observe the influence of surcharge (q) on backfill response, models are analyzed for three different magnitudes of imposed vertical stress. These magnitudes are q = 0 (no surcharge), q = 0.5 kPa, and q = 1 kPa (see Fig. 2). Surcharge magnitudes are selected based on the average magnitude of vertical stress in the retained backfill solely due to self-weight. The average vertical stress acts on the horizontal plane passing through the mid-height of the wall, and its magnitude is 1 kPa. The surcharge is applied on the backfill using ten separate rigid segments to achieve uniformity even during wall displacement and the associated backfill deformation, see Fig. 2.

There exist many uncertain parameters that can affect the behavior of the granular backfill, of which we have shortlisted the most important factors as particle shape, boundary conditions, and stress state. All other variables are kept constant while varying the focused factor to quantify the influence of each parameter. However, uncertainty associated with the selected factors must be considered in future studies. Various examples of uncertainty examples can be seen in Refs. [63–65].

2.3 Determination of failure surface geometry

In a sheared particulate assembly, Oda and Kazama [66] observed that the maximum particle rotations occur within shear bands, which is later confirmed by Zheng et al. [67]. Accordingly, this study uses localized deformations as detected from the analyses of cumulative rotations to determine the geometries and positions of failure surfaces. The cumulative rotation θ_i and its magnitude θ_i for each particle i, is calculated using Eq. (6) and then recorded for the full analysis.

(6){\mathit{\theta }}_i\left(t\right)={\displaystyle \frac{\mathrm{1}}{\mathrm{2}\mathrm{\pi }}}{\displaystyle \underset{\mathrm{0}}{\overset{t}{\int }}}{\mathit{\omega }}_i\mathrm{d}\tau ,\mathrm{\hspace{0.17em}}{\theta }_i=\left|\right|{\mathit{\theta }}_i\left|\right|,\mathrm{\hspace{0.17em}}\mathrm{\hspace{0.17em}}i=\overline{\mathrm{1},N_{\mathrm{p}}}\mathrm{,}

where ω_i is the angular velocity of particle i and N_p is the number of particles in the system.

The localized deformation region, obtained from cumulative rotation distribution, can be quantified by establishing a coordinate system. This fixed coordinate system is shown in Fig. 4 and its origin is chosen to coincide with the bottom of the model wall at its initial position. Dimensionless distances are obtained by normalizing the axes with model wall height (H_W). The geometry of the active failure surface is determined based on the measurements of the coordinates of points along the main shear surface, as illustrated in Fig. 4(a). In this respect, coordinates of the points along the failure surfaces are identified as either H (parallel to the model wall) or B (away from the model wall), as shown in Fig. 4(a). The horizontal distance between the retaining wall and the point at which the failure surface emerges at the ground surface is B_f (Fig. 4(b)). Note the measurements to determine the geometry of the failure surface are always conducted at S_average/H_w= 0.003 for all models, at which the backfill reaches a critical state (a more distinct failure wedge is observed at this instance).

2.4 Validation of the numerical model

The purpose of this work is to conduct a comparative modeling study using DEM to reveal the influences of grain shape, wall movements, and uniform surcharge on backfill response at different scales. Validation of the proposed model is needed to understand whether the outputs of models with the assigned parameters lead to realistic responses. Accordingly, results of a previously conducted small-scale physical model are used for comparison. The details of the physical models used in comparison can be found in Altunbas et al. [26]. Figure 5 shows the results of failure surface comparison. This comparison aims not to achieve a perfect fit as this is not a calibration process. General trends and forms of the failure surfaces obtained from DEM models are in good agreement with those obtained in physical model tests.

3 Results and discussions

3.1 Failure surface geometry

This section attempts to highlight the dependency of failure surface characteristics on particle and boundary properties comparing the quantified failure wedges for all backfills. For this purpose, failure wedges are visualized for various modes of wall movement at different surcharge levels for two different shape particles.

The differences in the evolution of failure wedges for different modes of wall movement are shown in Fig. 6 (note that these results are only shown for backfills with M2 particles when q = 0). For the initial stages of wall movement/rotation (Figs. 6(a) and 6(b)), it is noticed that deformation localization is different for different modes of wall movement; for HT: localization initiates at the toe of the wall and spreads toward the backfill surface; for RB: in the early stages of wall rotation, the upper part of the backfill near the backfill surface contains particles with the highest cumulative rotation magnitudes; for RT: deformations are predominantly localized near the toe of the wall from which they spread toward higher positions. As the wall is further translated/rotated (Figs. 6(c) and 6(d)), the failure wedge becomes more apparent (deformations are more localized). After a threshold in the wall translation/rotation is reached (Fig. 6(d), S_average/H_w= 0.03), further continuation of the wall movement does not change the final wedge geometry.

These differences in deformation mechanisms result in distinct failure surfaces. The dissimilarities are further explored by determining active failure surfaces for all models. These are presented in Fig. 7. Initially, the forms of the failure surfaces determined for models with M1 particles (Fig. 7(a)) are evaluated. It is noted in Fig. 7(a) that general forms of active failure surfaces are dependent on the model of wall movement. For HT models, visualized failure surfaces have parabolic forms, similar to the observations of Altunbas et al. [26]. In results of RT models, a more curved lower section is evident (still parabolic), whereas, for RB, the upper section of the failure surface reaches the backfill surface at a farther distance (almost linear). This form is a consequence of the differences in the progression of the failure mechanisms. When RB models are considered, active failure is reached close to the backfill surface at the initial stages of wall rotation. As the wall continues to rotate about the bottom, the active zone increases progressively from the top toward the bottom.

The bigger B_f values observed in RB mode of wall movement originates from the fact that particles, in this case, undergo an uneven displacement (i.e., particles at the top part experience a much larger displacement compared to particles near the bottom of the wall). This non-uniformity of wall displacement along X-axis leads to maximum deformations at the top of the wall (more extensive than those for RT and HT) and narrower distances at the base. This is a consequence of the mechanism in RB models, which creates an active zone that is increasing progressively from the top toward the bottom. The upper portions reach an active state at small magnitudes of rotation, and lower sections follow suit as the wall continues to rotate about the bottom. However, the mechanism is very different in RT models; soil starts to fail at the bottom even at small magnitudes of wall rotation, and failure wedge grows from the bottom toward the top. With this mechanism, later stages of active failure are contained by the initially formed larger wedge, and the resulting magnitudes B_f are similar to HT models. Ultimately, forms of the failure surfaces are more linear in RB models, whereas they are more curved in the case of RT models with greater curvature for lower sections. These results are similar to those experimentally observed by Yang and Tang [29].

The only impact of the surcharge (q) application on the geometries of failure surfaces is the change in the magnitude of B_f. Higher magnitudes of applied q yield greater B_f values for the resulting failure surfaces. This phenomenon can be attributed to the fact that the presence of surcharge leads to smaller peak friction angles [68,69], which results in relatively larger failure wedges.

When the failure surfaces of models with M2 backfills are examined, it is noted that the general forms are the same as M1 models for all HT, RT, and RB modes of wall movement. Thus, it is concluded that the influence of the mode of wall movement over the active failure surface is independent of particle shape. But the failure wedges are relatively smaller, as for M2 models particles are nonspherical. Interlocking is enhanced as the particle shape deviates from sphericity. Enhanced interlocking results in greater peak friction angles [70] and leads to narrower failure wedges. The observed results are comparable to experimental studies of Arda and Cinicioglu [71], where it is shown that backfills with more spherical particles yield larger active wedges when all other influences are the same.

3.2 Lateral earth pressure coefficients

The size and shape of the active failure wedge are indicators of active thrust. But the magnitude of lateral thrust on the wall can also be influenced by other factors, such as the mode of wall displacement, level of surcharge, and the backfill characteristics. Accordingly, this study examines the effects of wall movement mode and surcharge modes on normalized resultant lateral earth pressure coefficient. Lateral earth pressure coefficient (K) is the ratio of the horizontal component of the total force acting on the wall (F_h) to the vertical burden force at the bottom of the wall (F_v). Using the ratio of forces for calculating K yields an average K value for the entire depth of the wall. and F_v are calculated using Eqs. (7) and (8), respectively:

(7)F_{\mathrm{h}}={\displaystyle \frac{F_{\mathrm{n}}}{\cos \left(\delta \right)}}\left[\cos \right(\delta \pm \xi )+\sin (\delta )\times \sin (\pm \xi \left)\right],

(8)F_{\mathrm{v}}=\left[\left(q\times H_{\mathrm{w}}\right)+\left({\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}\times {\gamma }_{\mathrm{backfill}}\times H_{\mathrm{w}}^{\mathrm{2}}\right)\right]\times W_{\mathrm{w}}\times \cos \left(\xi \right).

In Eq. (7), F_h is the total normal force of all the contact points between the particles and the wall, δ is the interface friction angle between the backfill and the wall, and ξ is the rotation angle of the wall. In Eq. (7), ξ is added to δ in the RB mode while being subtracted from it in the RT mode (this is shown schematically in Fig. 3). ξ equals zero in the translational movement of the wall (HT). In Eq. (8), γ_backfill is the unit weight of the backfill. In the models without surcharge loads, q in Eq. (8) equals zero. γ_backfill is a function of grain density (ρ_grain), gravitational acceleration (g) and porosity (n) as given in Eq. (9):

(9){\gamma }_{\mathrm{backfill}}={\rho }_{\mathrm{grain}}\times g\times (\mathrm{1}-n).

Before any wall movement, the value of K corresponds to the at-rest coefficient of lateral earth pressure (K₀). Values of K₀ for all models are reported in Table 4. It should be noted here that the magnitude is common for all tests with the same particle shape and surcharge levels. When Table 4 is inspected, it is observed that the magnitude of K₀ for backfills with M1 particles when at q = 0 is higher than 1. This high magnitude of K₀ is a consequence of the material preparation method, wherein the particles are pluviated by setting backfill, which is prepared by setting the inter-particle sliding friction coefficient (μ_s) zero. Following the filling process, the friction coefficient between the particles is set to its final value (0.68) before the commencement of wall movement. At this state, the prepared sample is an overconsolidated backfill, and thus a higher K₀ value is expected. As shown in Table 4, surcharge application before wall movement causes a significant decrease in K₀, as this process reduces the over-consolidation ratio (OCR). On the other hand, backfills with M2 particles are prepared using a higher μ_s during the filling process (μ_s = 0.17). Otherwise, it would not be possible to achieve the same porosity as backfills with M1 particles. This is because deviation from sphericity increases the propensity for forming denser assemblies; hence less energy must be expended during model preparation with nonspherical particles to achieve the same porosity as an assembly composed of spherical particles. The use of less energy in model preparation, in a sense, leads to smaller pre-consolidation pressures. That is the reason why the magnitude of K₀ is lower for models composed of M2 particles. Furthermore, surcharge application also reduces the magnitude of K₀ as it also leads to a reduction in OCR. Table 5 presents values of K_a for all models of this study. Apparently, the magnitude of K_a at residual state is practically independent of surcharge but varies with the mode of wall movement. This is expected since, at the residual state, any interlocking that was present before is eliminated. So, only one average K_a value is presented for all different surcharge levels (0, 0.5, and 1 kPa).

Figure 8 for M1 models and Fig. 9 for M2 models present the variations of lateral earth pressure coefficient concerning wall movement (HT, RB, and RT) under various surcharge loads (0, 0.5, and 1 kPa) as the wall moves toward active failure. The range of K_a values for both backfills (M1 and M2) are compatible with those reported in previous Refs. [40,41]. With the evaluation of the results presented in Figs. 8 and 9, the influences of surcharge, mode of wall displacement, and particle shape can be considered in detail. Both for Figs. 8 and 9, the first three plots (Figs. 8(a), 8(b), 8(c), 9(a), 9(b), and 9(c)) show K-S_average/H_w relationships obtained under different surcharge levels for different modes of wall movements. These results, apart from insignificant differences, are practically the same. This finding is correct for both M1 and M2 models and all modes of wall movement.

On the other hand, when the results in the last three plots of Fig. 8 and Fig. 9 (Figs. 8(d), 8(e), 8(f), 9(d), 9(e), and 9(f)) are examined, it is observed that the mode of wall displacement influences the value of K for both spherical and nonspherical particles. As shown in Figs. 8(d)–8(f) for M1 backfills, both the ultimate values of K (K_a) and the shapes of the K-S_average/H_w relationships are dependent on the mode of wall movement. Regarding K_a, it is interesting to note that the RT mode of wall displacement yields greater values, whereas RB and HT modes yield the same value for all levels of surcharge. This is a consequence of the fact that RT mode of wall movement provides the most stable mechanism of lateral wall motion toward active state; RT mode is in a sense similar to a trapdoor experiment for which the soil flows through the fixed-top position of the wall and the stationary section below the wall. Therefore, the arching mechanism is more substantial in RT mode, as is discussed in the following sections of the paper. Yet, for M2 models, RT and HT modes produce the same K_a, whereas K_a of RB is smaller. RB mode of wall movement provides the smoothest transition to an active state for which the active zone is increasing progressively from the top toward the bottom, as explained in the previous section. However, the magnitude of K_a for HT mode increases to the level in RT mode, owing to the stabilizing effect of the nonspherical particles.

Similarly, for K-S_average/H_w relationships, the value of K directly drops to its ultimate value, when in fact, the magnitude of K drops to a minimum before slightly increasing to the ultimate value for RB and HT models. This is also an outcome of the progressively expanding active zone in the RT mode of wall movements. Nadukuru and Michalowski [40] presented a similar trend for K variation, even though their backfill is a mixture of spherical and clump particles (30% of non-spheres and 70% of spheres).

3.3 Normalized lateral force distribution along the wall

In DEM, discrete elements exert forces on the wall. Thus, it is necessary to divide the wall into sections and calculate the normalized lateral force (P_L) for each section. For this purpose, the wall is divided into five sections along with its height; the magnitude of total normal force for each section is obtained and normalized by the vertical burden force that acts at the bottom of the wall. In the absence of surcharge, the vertical burden force is calculated using Eq. (8), assuming q = 0. Note that the normalization in the presence of surcharge is carried out using Eq. (8), with q = 1 kPa. This normalization allows a comparison between different models based on pressure distribution. Additionally, it is necessary to note that P_L values are averaged over a certain period of wall movement to reduce temporal fluctuations. For this purpose, results are averaged over various successive saving time-steps. Based on the obtained mean values, an optimum window of 50 saving time-steps is chosen.

Figures 10(a) and 11(a) present the distribution of normalized lateral forces acting on the wall in M1 models for q = 0 and q = 1 kPa conditions, respectively (results for the case with q = 0.5 kPa are not presented in this section, since q = 1 kPa represents well the surcharge effect). Only a single at rest normalized lateral force distribution is used in these figures as all the samples have identical initial packing. Similarly, Figs. 12(a) and 13(a) are prepared for the same purpose, but this time for M2 models. In addition to normalized force distribution plots, Figs. 10, 11, 12, and 13 also present the established force networks at the active state for different modes of wall movement. Visualization of force networks is necessary to understand better and interpret the results. In the following paragraphs, the relationship between mode of wall movement, resulting in force networks, and distribution of lateral earth pressures on the wall.

As expected, for all three modes of wall movement, earth pressures at the active state are significantly lower than the at-rest lateral earth pressures. In the case of HT, P_L distribution has one maximum point at the mid-height of the wall. This can be attributed to strong force chains connecting the midsection of the wall to the stationary section at the toe level, as shown in Fig. 10(b). The force chains form mechanisms similar to arching. However, owing to the constant translation of the wall, force chains constantly disintegrate and reform. The bottom end of the force chains is persistent, as indicated with the darker color in Fig. 10(a). But the top part adjusts to the nonstationary boundary; therefore, it cannot carry as great magnitudes of stress and is shown with red color. Therefore, it is more appropriate to refer to these force chains as pseudo-arches.

The behavior of the RT case (see Fig. 10(c)) is similar to HT, but the pseudo-arches that originate at the bottom of the wall reach higher positions on the wall. This is a result of wall rotation about the top. With RT mode of movement, the top portions displace less, therefore provide relatively more stable positions to transfer loads. As a consequence of this mechanism, pseudo-arches reach higher positions (Fig. 10(c)), and the maximum point of P_L distribution (Fig. 10(a)) is at a higher position along the wall, relative to the HT case. Nadukuru and Michalowski [40] observed similar P_L distribution for RT models at the active state.

On the other hand, for the RB mode of wall translation, the magnitude of P_L at the active state increases by depth. With the RB mode of wall movement, higher positions on the wall displace more, making it more challenging to form force networks that transfer loads on the wall. This condition can be viewed in Fig. 10(d), which clearly shows the absence of arching. The bottom portion of the wall forms a more stable position to transfer backfill forces, hence the greater magnitudes of P_L toward the bottom of the wall.

Once surcharge is applied on the granular assembly, strong contact forces form between the toe of the wall and the backfill surface, see Fig. 11. This occurs because surcharge serves as a reaction to the forces transferred by the force chains networks. Yet arching between the retaining wall and the toe still exists in both RT and HT modes, similar to the cases without surcharge. It is also clear that applying surcharge does not change the general form of lateral pressure distribution.

With nonspherical M2 particles, the magnitude of P_L at the active state reduces, as observed in Figs. 12(a) and 13(a). This is expected as the deviation from sphericity increases friction angle [60], and higher friction angles yield narrower failure wedges that produce smaller active thrusts. The obtained results are in line with those reported by Jiang et al. [41], where it was shown that applying rolling resistance (as an alternative for particle shape factor) reduces the exerted loads on the wall.

The effect of particle shape on the distribution and effectiveness of force chain networks can be examined by comparing Figs. 10 and 11 with Figs. 12 and 13. Backfills with M2 particles produce less prominent pseudo-arches. As seen before in Figs. 10 and 11, backfills with spherical particles have better-defined pseudo-arching formations. This observation is consistent with the results obtained by Chen et al. [45]. They hypothesized that circular particles have fewer contact points and underestimated interlocking behaviors; thus the force redistribution process and particle rearrangement are enhanced at the microscale. In other words, having irregularly shaped particles restrict particular rearrangement, and arching is hindered.

3.4 Coordination number

Coordination number is a micro-scale variable that shows contact density within granular materials. It defines the average number of particles that a single particle is in contact with. The average coordination number of a granular assembly (CN_ave) is not only a micro-scale indicator but also has an influence over macro-scale behavior [72]. For example, an increase in CN_ave can lead to higher shear stiffness in sands [73] and more stable embankment filters [74]. Therefore, it will be insightful to investigate the magnitude of CN_ave and its variation with particle shape, wall movement mode and surcharge application.

Figure 14 shows the results of simulations involving backfills with M1 (spherical) particles. According to Fig. 14, neither wall movement mode nor surcharge acting on the backfill affects CN_ave before and after wall movement. Similarly, Gezgin et al. [75] showed that compaction of an arrangement of spherical particles with no rolling resistance does not change CN_ave.

On the other hand, the behavior of backfills with M2 (nonspherical) particles is different when CN_ave is concerned. First of all, a comparison between Figs. 14(a) and 15(a) shows that assemblies composed of nonspherical particles (M2) have higher CN_ave before wall movement, even though backfills with spherical particles (M1) (μ_s = 0 during filling) have higher relative densities, owing to the lower μ_s used during pluviation (μ_s = 0 in M1 versus μ_s = 0.17 in M2) under the same boundary conditions. Secondly, Fig. 15(d) indicates that the application of surcharge on a backfill with M2 particles causes an increase in CN_ave value before any wall movement, unlike the case of backfills with M1 particles (see Fig. 14(d)).

Additionally, the difference in CN_ave among backfills assembled with M2 particles for different surcharge levels remains almost constant during wall movement, as shown in Fig. 15(d). The reason for the difference between M1 and M2 models is the presence of interlocking in M2 models. Interlocking hinders the freedom of movement in the case of M2 models. The changes in the values of CN_ave with wall movement may seem small (Fig. 15(d)). However, it must be emphasized that the value of CN_ave is an average for the whole model, whereas the deformation zone (active zone) is just a small portion of the assembly. Thus, the actual change in the value of CN_ave for the actual deformation zone is much more significant.

Overall, it can be commonly deduced from Figs. 14 and 15 that, regardless of the particle shape (M1 or M2), the mode of wall movement does not affect the change in the value of CN_ave.

4 Conclusions

This study considers the influence of the mode of wall movement on the resulting active state for different particle shapes and surcharge levels. For this purpose, DEM models are developed, and the evolution of the failure surface for each backfill is tracked by monitoring the variation of cumulative rotation distribution. The results are analyzed and discussed at different scales to better understand the interactions between macro and micro influences. The main conclusions are as follows.

1) It is confirmed that the geometry of the active failure wedge is a function of the mode of wall movement. RT mode results in curved failure surfaces, whereas RB mode produces more linear ones. This is an outcome of the difference in the progression of failure; with RT mode, the mechanism is bottom toward the top with which later stages of active failure are contained by the initially formed larger wedge; with RB mode, active failure expands progressively from the top toward the bottom; with HT, the mechanism is somewhere in between RT and RB, being more similar to RT mode.

2) The mode of wall movement affects both K–S_average/H_w relationship and the average magnitude of K_a, for both spherical and nonspherical particles. Among different modes of wall movement, RT yields the greatest magnitude of K_a, both for spherical and nonspherical particles. This is because RT mechanism provides more opportunity for the formation of more stable and longer force chains. On the other hand, RB mode of wall movement offers the smoothest transition to active state with a progressively expanding failure zone. Hence, RB models require the greatest S_average/H_w magnitude to reach an active state. It appears that the level of surcharge has a noticeable influence over neither the magnitude of K_a, nor K–S_average/H_w relationship.

3) The mode of wall movement is also influential over the distribution of active pressures. When the variations of active pressures along the wall height are considered, the maximum pressures are observed at about wall midheight for both HT and RT modes. This point is located higher in the case of RT mode. This is a consequence of the resulting force chain networks. RT mode provides more stable conditions toward the portions of the wall, and RB mode is on the opposite end of the spectrum. Therefore, force chain networks that originate at the bottom can grow longer and transfer forces to higher locations in the case of RT mode of movement.

4) Force chain networks generate and disintegrate throughout wall displacement. These networks are never stable but still affect pressure distribution on the wall due to the continuous movement of the wall. Therefore, these force chain networks can be referred to as pseudo-arches. More stable pseudo-arches are produced in the case of models with surcharge, extending to the ground surface, where surcharge provides a more stable reaction surface to transfer forces. It is also noticed that spherical particles, due to ease of particle rearrangement, yield more prominent arches within backfills.

5) The results show that the average coordination number CN_ave reduces with wall movement toward the active state, irrespective of mode. In the case of models with spherical particles, neither surcharge level nor the mode of wall movement has any influence over CN_ave–S_average/H_w relationship. But for models with nonspherical particles, CN_ave is higher under greater magnitudes of the surcharge due to interlocking. Yet, the forms of CN_ave–S_average/H_w relationships are still the same for all modes of wall movement.

As a final one, the current study can lead to optimization in design process of retaining structures against active state of failure by improving engineers’ understanding over: a) the geometry and location of failure surface for different modes of wall movement, which can be used to determine the correct length of supporting anchors in real-world projects; b) the force distribution along the wall and influence of arching phenomenon; c) influence of surcharge on the failure characteristics of the backfill. Additionally, the current study informs the engineers on the importance of considering particle shape effect on granular response.

References

[1] Coulomb C A. An attempt to apply the rules of maxima and minima to several problems of stability related to architecture. Memoirs of the Royal Academy of Sciences, 1776, 7: 343–382

[2] Rankine W J M. On the stability of loose earth. Philosophical Transactions of the Royal Society of London, 1857, 147: 9–27

[3] Terzaghi K. A fundamental fallacy in earth pressure computations. Boston Society of Civil Engineers Journal, 1936, 23(2): 71–88

[4] Tsagareli Z V. Experimental investigation of the pressure of a loose medium on retaining walls with a vertical back face and horizontal backfill surface. Soil Mechanics and Foundation Engineering, 1965, 2(4): 197–200

[5] Narain R J, Saran S, Nandakumaran P. Model study of passive pressure in sand. Journal of the Soil Mechanics and Foundations Division, 1969, 95(4): 969–984

[6] Fang Y S, Ishibashi I. Static earth pressures with various wall movements. Journal of Geotechnical Engineering, 1986, 112(3): 317–333

[7] Tang Z C. A rigid retaining wall centrifuge model test of cohesive soil. Journal of Chongqing Jiaotong University, 1988, 25(2): 48–56

[8] Wang J, Lin X, Chai H, Xu J. Coulomb-type solutions for passive earth pressure with steady seepage. Frontiers of Architecture and Civil Engineering in China, 2008, 2(1): 56–66

[9] Hamrouni A, Dias D, Sbartai B. Soil spatial variability impact on the behavior of a reinforced earth wall. Frontiers of Structural and Civil Engineering, 2020, 14(2): 518–531

[10] Wen H, Cheng Q, Meng F, Chen X. Diaphragm wall-soil-cap interaction in rectangular-closed diaphragm-wall bridge foundations. Frontiers of Architecture and Civil Engineering in China, 2009, 3(1): 93–100

[11] White D J, Take W A, Bolton M D. Soil deformation measurement using particle image velocimetry (PIV) and photogrammetry. Geothechnique, 2003, 53(7): 619–631

[12] Stanier S A, Blaber J, Take W A, White D J. Improved image-based deformation measurement for geotechnical applications. Canadian Geotechnical Journal, 2016, 53(5): 727–739

[13] Salehi Alamdari N, Khosravi M H, Katebi H. Distribution of lateral active earth pressure on a rigid retaining wall under various motion modes. International Journal of Mining and Geo-Engineering, 2020, 54(1): 15–25

[14] Patel S, Deb K. Study of active earth pressure behind a vertical retaining wall subjected to rotation about the base. International Journal of Geomechanics, 2020, 20(4): 04020028

[15] Paikowsky S G, Tien H S. Experimental examination of the arching mechanism on the micro level. In: Proceedings of the Third International Conference on Discrete Element Methods. New Mexico: ASCE Press, 2002, 222–228

[16] Fang Y, Guo L, Hou M. Arching effect analysis of granular media based on force chain visualization. Powder Technology, 2020, 363: 621–628

[17] He Z, Liu Z, Liu X, Bian H. Improved method for determining active earth pressure considering arching effect and actual slip surface. Journal of Central South University, 2020, 27(7): 2032–2042

[18] Nuebel K. Experimental and numerical investigation of shear localization in granular material. Dissertation for the Doctoral Degree. Baden Wurttemberg: Karlsruhe University, 2002

[19] Niedostatkiewicz M, Lesniewska D, Tejchman J. Experimental analysis of shear zone patterns in cohesionless for earth pressure problems using particle image velocimetry. Strain, 2011, 47: 218–231

[20] Gutberlet C, Katzenbach R, Hutter K. Experimental investigation into the influence of stratification on the passive earth pressure. Acta Geotechnica, 2013, 8(5): 497–507

[21] Soltanbeigi B, Altunbas A, Cinicioglu O. Influence of dilatancy on shear band characteristics of granular backfills. European Journal of Environmental and Civil Engineering, 2021, 25(7): 1201–1218

[22] Pietrzak M, Leśniewska D. Strains inside shear bands observed in tests on model retaining wall in active state. In: Micro to Macro Mathematical Modelling in Soil Mechanics. Reggio Calabria: Springer, 2018, 257–265

[23] Leśniewska D, Nitka M, Tejchman J, Pietrzak M. Contact force network evolution in active earth pressure state of granular materials: Photo-elastic tests and DEM. Granular Matter, 2020, 22(3): 1–31

[24] Khosravi M H, Pipatpongsa T, Takemura J. Experimental analysis of earth pressure against rigid retaining walls under translation mode. Geotechnique, 2013, 63(12): 1020–1028

[25] Liu F Q. Lateral earth pressures acting on circular retaining walls. International Journal of Geomechanics, 2014, 14(3): 04014002

[26] Altunbas A, Soltanbeigi B, Cinicioglu O. Determination of active failure surface geometry for cohesionless backfills. Geomechanics and Engineering, 2017, 12(6): 983–1001

[27] Li M G, Chen J J, Wang J H. Arching effect on lateral pressure of confined granular material: Numerical and theoretical analysis. Granular Matter, 2017, 19(2): 20

[28] Zhou Q Y, Zhou Y T, Wang X M, Yang P Z. Estimation of active earth pressure on a translating rigid retaining wall considering soil arching effect. Indian Geotechnical Journal, 2018, 48(3): 541–548

[29] Yang M, Tang X. Rigid retaining walls with narrow cohesionless backfills under various wall movement modes. International Journal of Geomechanics, 2017, 17(11): 04017098

[30] Cundall P A, Strack O D L. A discrete numerical model for granular assemblies. Geotechnique, 1979, 29(1): 47–65

[31] Rabczuk T, Zi G, Bordas S, Nguyen-Xuan H. A simple and robust three-dimensional cracking-particle method without enrichment. Computer Methods in Applied Mechanics and Engineering, 2010, 199(37–40): 2437–2455

[32] Rabczuk T, Belytschko T. Cracking particles: A simplified meshfree method for arbitrary evolving cracks. International Journal for Numerical Methods in Engineering, 2004, 61(13): 2316–2343

[33] Rabczuk T, Belytschko T. A three-dimensional large deformation meshfree method for arbitrary evolving cracks. Computer Methods in Applied Mechanics and Engineering, 2007, 196(29–30): 2777–2799

[34] Ren H, Zhuang X, Cai Y, Rabczuk T. Dual-horizon peridynamics. International Journal for Numerical Methods in Engineering, 2016, 108(12): 1451–1476

[35] Ren H, Zhuang X, Rabczuk T. Dual-horizon peridynamics: A stable solution to varying horizons. Computer Methods in Applied Mechanics and Engineering, 2017, 318: 762–782

[36] Jiang M, Zhu H, Li X. Strain localization analyses of idealized sands in biaxial tests by distinct element method. Frontiers of Architecture and Civil Engineering in China, 2010, 4(2): 208–222

[37] Raychowdhury P, Jindal S. Shallow foundation response variability due to soil and model parameter uncertainty. Frontiers of Structural and Civil Engineering, 2014, 8(3): 237–251

[38] Guo Y, Yu X B. Design and analyses of open-ended pipe piles in cohesionless soils. Frontiers of Structural and Civil Engineering, 2016, 10(1): 22–29

[39] Tejchman J, Kozicki J, Leśniewska D. Discrete simulations of shear zone patterning in sand in earth pressure problems of a retaining wall. International Journal of Solids and Structures, 2011, 48(7–8): 1191–1209

[40] Nadukuru S S, Michalowski R L. Arching in distribution of active load on retaining walls. Journal of Geotechnical and Geoenvironmental Engineering, 2012, 138(5): 575–584

[41] Jiang M, He J, Wang J, Liu F, Zhang W. Distinct simulation of earth pressure against a rigid retaining wall considering inter-particle rolling resistance in sandy backfill. Granular Matter, 2014, 16(5): 797–814

[42] Yang M, Deng B. Simplified method for calculating the active earth pressure on retaining walls of narrow backfill width based on dem analysis. Advances in Civil Engineering, 2019, 2019: 1–12

[43] Albaba A, Lambert S, Nicot F, Chareyre B. Relation between microstructure and loading applied by a granular flow to a rigid wall using DEM modeling. Granular Matter, 2015, 17(5): 603–616

[44] Nitka M, Tejchman J, Kozicki J. Discrete modelling of micro-structural phenomena in granular shear zones. In: Proceedings of the 11th International Workshop on Bifurcation and Degradation in Geomaterials. Limassol: Springer, 2014, 7–12

[45] Chen R P, Liu Q W, Wu H N, Wang H L, Meng F Y. Effect of particle shape on the development of 2D soil arching. Computers and Geotechnics, 2020, 125: 103662

[46] Wu J, Wang L, Cheng Q. Soil arching effect of Lattice-Shaped Diaphragm Wall as bridge foundation. Frontiers of Structural and Civil Engineering, 2017, 11(4): 446–454

[47] Khosravi M H, Azad F H, Bahaaddini M, Pipatpongsa T. DEM analysis of backfilled walls subjected to active translation mode. International Journal of Mining and Geo-Engineering, 2017, 51(2): 191–197

[48] O’Sullivan C. Particulate Discrete Element Modelling: A Geomechanics Perspective. 1st ed. London: Tylor and Francis Group, 2011

[49] EDEM. EDEM 2018 User Guide, 2018

[50] Talebi H, Silani M, Bordas S P A, Kerfriden P, Rabczuk T. A computational library for multiscale modeling of material failure. Computational Mechanics, 2014, 53(5): 1047–1071

[51] Budarapu P R, Gracie R, Bordas S P A, Rabczuk T. An adaptive multiscale method for quasi-static crack growth. Computational Mechanics, 2014, 53(6): 1129–1148

[52] de Bono J P, McDowell G R. DEM of triaxial tests on crushable sand. Granular Matter, 2014, 16(4): 551–562

[53] de Bono J P, McDowell G R. On the micro mechanics of yielding and hardening of crushable granular soils. Computers and Geotechnics, 2018, 97(January): 167–188

[54] de Bono J P, Li H, McDowell G R. A new abrasive wear model for railway ballast. Soil and Foundation, 2020, 60(3): 714–721

[55] Grabowski A, Nitka M, Tejchman J. 3D DEM simulations of monotonic interface behaviour between cohesionless sand and rigid wall of different roughness. Acta Geotechnica, 2021, 16(4):1001–1026

[56] de Bono J P, McDowell G R. Micro mechanics of the critical state line at high stresses. Computers and Geotechnics, 2018(98): 181–188

[57] de Bono J P, McDowell G R. Micro mechanics of drained and undrained shearing of compacted and overconsolidated crushable sand. Geotechnique, 2017, 68(7): 575–589

[58] Yamamuro J A, Lade P V. Drained sand behavior in axisymmetric tests at high pressures. Journal of Geotechnical Engineering, 1996, 122(2): 109–119

[59] Bandini V, Coop M R. The influence of particle breakage on the location of the critical state line of sands. Soils and Foundations, 2011, 51(4): 591–600

[60] Cho G C, Dodds J, Santamarina J C. Particle shape effects on packing density, stiffness, and strength: Natural and crushed sands. Journal of Geotechnical and Geoenvironmental Engineering, 2006, 132(5): 591–602

[61] Soltanbeigi B, Papanicolopulos S, Ooi J Y. Particle shape effect on deformation localization during quasi-static flow at active state. In: 5th International Conference on Geotechnical Engineering and Soil Mechanics. Tehran, 2016

[62] Soltanbeigi B, Podlozhnyuk A, Kloss C, Pirker S, Ooi J Y, Papanicolopulos S A. Influence of various DEM shape representation methods on packing and shearing of granular assemblies. Granular Matter, 2021, 23(2): 26

[63] Hamdia K M, Silani M, Zhuang X, He P, Rabczuk T. Stochastic analysis of the fracture toughness of polymeric nanoparticle composites using polynomial chaos expansions. International Journal of Fracture, 2017, 206(2): 215–227

[64] Hamdia K M, Ghasemi H, Zhuang X, Alajlan N, Rabczuk T. Sensitivity and uncertainty analysis for flexoelectric nanostructures. Computer Methods in Applied Mechanics and Engineering, 2018, 337: 95–109

[65] Vu-Bac N, Lahmer T, Zhuang X, Nguyen-Thoi T, Rabczuk T. A software framework for probabilistic sensitivity analysis for computationally expensive models. Advances in Engineering Software, 2016, 100(October): 19–31

[66] Oda M, Kazama H. Microstructure of shear bands and its relation to the mechanisms of dilatancy and failure of dense granular soils. Geotechnique, 1998, 48(4): 465–481

[67] Zheng H, Wang D, Tong X, Li L, Behringer R P. Granular scale responses in the shear band region. Granular Matter, 2019, 21(4): 1–6

[68] Fannin R J, Eliadorani A, Wilkinson J M T. Shear strength of cohesionless soils at low stress. Geotechnique, 2005, 55(6): 467–478

[69] Cinicioglu O, Abadkon A. Dilatancy and friction angles based on in situ soil conditions. Journal of Geotechnical and Geoenvironmental Engineering, 2015, 141(4): 06014019

[70] Amirpour Harehdasht S, Hussien M N, Karray M, Roubtsova V, Chekired M. Influence of particle size and gradation on shear strength–dilation relation of granular materials. Canadian Geotechnical Journal, 2019, 56(2): 208–227

[71] Arda C, Cinicioglu O. The influences of micro-mechanical properties of cohesionless soils on active failure surface geometries. Teknik Dergi, 2019, 30(5): 9399–9420 (in Turkish)

[72] Fei W, Narsilio G A. Impact of three-dimensional sphericity and roundness on coordination number. Journal of Geotechnical and Geoenvironmental Engineering, 2020, 146(12): 06020025

[73] Fonseca J, O’Sullivan C, Coop M R, Lee P D. Quantifying the evolution of soil fabric during shearing using scalar parameters. Geotechnique, 2013, 63(10): 818–829

[74] Fonseca J, Sim W W, Shire T, O’Sullivan C, Wang Y, Dallo Y A H. Microstructural analysis of sands with varying degrees of internal stability. Geotechnique, 2015, 65(7): 620–623

[75] Gezgin A T, Soltanbeigi B, Cinicioglu O. Discrete-element modelling of pile penetration to reveal influence of soil characteristics. Proceedings of the Institution of Civil Engineers- Geotechnical Engineering, 2020, 1–42

RIGHTS & PERMISSIONS

Higher Education Press