Abstract

In this study, the nonlocal, nonordinary state-based peridynamics (NOSBPD) is introduced as a novel regularization technique to study the strain localization and progressive failure of soil slopes due to strain softening behavior. The NOSBPD formulations are introduced first, and then the strain-softening Mohr–Coulomb model is incorporated into the nonlocal framework. The implicit formulations of the nonlinear NOSBPD model under the plane strain condition combined with the energy dissipation-based arc-length method to avoid the limitations of the force and displacement control in the nonlinear analysis are given. Numerical examples, including the plane strain biaxial experiments and slope models, validate the effectiveness of the proposed method in terms of alleviating mesh sensitivity. The results show that this nonlocal method can remove the mesh dependence, and the localized shear band is related to the nonlocal parameter δ. The proposed method involves being able to capture strength mobilization in different parts of the sliding surface during the progressive failure process of the soil slope. It can avoid the dilemma that the peak strength is too dangerous and the residual strength is too conservative due to the rigid body assumption in the limit equilibrium method.

1. Introduction

Strain localization is a common phenomenon in geotechnical engineering, and its occurrence is often accompanied by strain softening and progressive failure of structures. Regarding the study of strain localization in geomaterials, many attempts have been made in experimental, theoretical, and numerical aspects. In terms of experiments [13], there are many studies on soil samples with significant strain-softening behavior, such as dense sand and hard clay, including the plane strain test and the centrifugal model test. From the theoretical point of view [4, 5], the occurrence of shear bands is commonly predicted by the bifurcation theory. In recent years, the introduction of noncoaxial theory [6] has made the prediction more effective, but bifurcation theory cannot simulate the behavior of materials after bifurcation. For the numerical simulation methods for crack initiation and propagation in brittle or quasi-brittle material, there are two types of methods, namely separation and smearing methods. The numerical method for the strain localization analysis can also be roughly divided into separated and continuous models according to the relationship between shear band width (h) and element size (d). Separation models can be further divided into weak discontinuities [7], strong discontinuities [8], shear elements [9, 10], etc. In this method, the relationship between the element size and the width of the shear band is generally expressed as d ≥ h (in the strong discontinuities method, h = 0). However, the continuous model is generally required d < h. The current study focuses on continuous methods.

For the conventional finite element method without internal length scale, when it deals with strain localization problems, the static control equations are transformed from elliptical to hyperbolic, and the dynamic control equations are changed from hyperbolic to elliptical, due to the strain softening or asymmetric stiffness matrix obtained by the nonassociative flow rule. Thus, the solution to the boundary value problem depends heavily on the size of the finite element. In particular, as the size of the finite element decreases, the width of the shear band decreases, which means that with the mesh refinement, the energy consumption would approach zero. It obviously violates the laws of physics. Therefore, considerable regularization mechanisms have been introduced to overcome the mesh sensitivity, such as viscous models [11], the Cosserat continuum model [12], the gradient plasticity model [13], and the nonlocal model [14, 15], among others.

The peridynamic (PD) theory originally proposed by Silling [16] uses the form of differential-integral to describe the equation of motion of the continuum medium, avoiding the displacement derivative term, so it can obtain a continuous-discontinuous unified expression. The PD theory, including PD-FEM coupled models, has great advantages in dealing with discontinuity failure problems, as illustrated in previous studies [1719]. At the same time, PD naturally has the nonlocal parameter, i.e., neighborhood radius (or called the horizon, δ). Thus, it can seem like a consistent nonlocal model in the sense that in some existing nonlocal models, only nonlocal averaging of stress-strain is performed in the postprocessing stage, while constitutive models are still local models. The peridynamic theory is widely used in composite materials [20], concrete [21], and other fields [22], but it is rarely seen in the field of predicting the progressive failure process of soil slopes. Song et al. [23, 24] initially applied peridynamics for the analysis of chemical plastic strain localization in unsaturated soils, but the strain softening properties, one of the main factors causing strain localization [2527], have not yet been discussed thoroughly.

The main objective of the current study is to use the PD theory as a regularization technique for the strain localization analysis in soil induced by strain-softening behavior. The effectiveness of this method for relieving the mesh dependency is emphatically analyzed by a plane strain element test and a slope model, and the influence of a neighborhood radius of δ is discussed. The analysis in this paper is based on a small deformation assumption and does not consider soil-water interactions. The rest of the manuscript is organized as follows: In Section 2, the governing equations of the numerical model and the constitutive model are introduced. The plane strain test is simulated by the proposed method in Section 3. Numerical studies for the progressive failure of soil slopes are conducted in Section 4. A summary and conclusions are drawn in Section 5.

2. Numerical Method and Constitutive Model

2.1. Nonordinary State-Based Peridynamics Theory

The state-based peridynamics theory was proposed by Silling to overcome the limitations of the fixed Poisson’s ratio of the bond-based peridynamic theory [27]. By introducing the state concept, the existing constitutive models expressed in terms of stress and strain quantities in classical continuum mechanics can be easily incorporated into the nonordinary state-based peridynamics (NOSBPD) framework. For self-completeness, the basic formulations of NOSBPD are briefly outlined in this subsection. More details about this theory can be found in [28].

The notations and kinematics of NOSBPD are illustrated in Figure 1. The reference position vector state is defined as

As deformation proceeds, the deformed position vector state is given by

The nonlocal deformation gradient is approximated aswhere the shape tensor is defined aswhere Hx denotes the horizon of a material point x. It is a circle in two dimensions with a cut-off radius of δ.

Based on the principle of virtual work, the force state can be derived as follows in terms of the stress tensor as an intermediate step:

Subsequently, the governing equation of the NOSBPD in the static or quasi-static conditions ignoring the inertial effects reads

2.2. Strain-Softening Constitutive Model for Soil

In the current study, to describe the strain-softening behavior of the soil, the modified Mohr–Coulomb model, considering a variation in the cohesion strength c with the accumulated plastic strain is employed. The celebrated Mohr–Coulomb constitutive model has been widely used for describing the shear failure of soils, and it can reflect the soil failure dependence on the hydrostatic pressure. The dilatancy can be reflected by the nonassociative flow law. It is noted that the model in this paper does not consider hat-shaped yielding under isotropic compression.

The yield criterion of this model (see Figure 2(a)) is given bywhere and are the major and minor principle of effective stress (the compression components of the stress are positive), respectively. c is the cohesion intercept, and φ is the friction angle. The flow rule is a nonassociated type with a dilatancy angle ψ, which satisfies ψ<φ.

Following the previous studies [14], the strain-softening behavior of soil can be described by decreasing the strength parameters (c and/or φ). In this study, the attenuation of the friction angle is not considered during softening (see Figure 2(b)), and the strain softening characteristics are solely reflected by the attenuation of the cohesion force with the development of the accumulated plastic strain. In specific, the cohesion intercept is linearly reduced with the accumulated plastic strain after the occurrence of the initial yield, as shown in Figure 2(c).

2.3. Numerical Implementation

The meshless method proposed by Silling is employed herein for spatial discretization [29]. That is, a particle with a certain volume at the center of the grid is used to represent the movement of the entire cell. The explicit solution algorithm is often used for peridynamics, while the time step of the explicit algorithm is commonly very short and it takes a long time for the structure to reach a stable state. In this study, a nonlinear implicit solution algorithm is used. It differs from the previous study [30] in that it is expressed in a nonlinear framework.

The shape tensor in a discretized form is given bywhere is the source point; is the field point; is the volume of the particle .

Correspondingly, the discretized deformation gradient tensor reads

With the infinitesimal deformation assumption, the strain tensor can be calculated as

Thus, the strain tensor can be approximated as

For the plane strain problem considered herein, the strain is expressed in a matrix form aswhere the matrix B, N, and U are

The stress tensor can be expressed aswhere denotes the state variables, herein it is the accumulated plastic strain.

Consequently, the force state can be rewritten aswhere the matrix is

It is noted that the imposition of boundary conditions for the nonlocal PD model is more complicated than the local models. In this paper, both the boundary displacement constraints and the external force are applied through a particle layer with a thickness of δ (δ is the neighborhood radius). The traction force is applied as a body force aswhere is the thickness of the boundary layer, and is the pressure.

The implementation method of the elastic-plastic model in the NOSBPD framework is almost identical to the local one [20]. The fully implicit backward Euler stress integration algorithm is employed for the strain-softening Mohr–Coulomb constitutive model [31]. It should be noted that in the Mohr–Coulomb model, the flow directions of the plastic potential at the corners of the plane are a combination of multiple yield surface normal vectors [31].

In the nonlinear analysis of strain softening materials, the force control loading method cannot simulate the descending section, and the displacement control cannot deal with the problem of snap-back, so the arc-length method needs to be used for control. Although the conventional arc-length method is effective for geometric nonlinear control, the arc-length governing equations for material nonlinear problems often have complex roots. Thus, it is difficult to identify the occurrence of local damage and failure. Therefore, the force control and the energy dissipation-based arc-length method are used jointly to control the nonlinear iteration process. That is, a force control constraint is used at first before switching to the energy dissipation-based arc-length method once the plastic energy dissipates. The details of this solution algorithm can be found in [32].

The aforementioned solution scheme has been incorporated into a C++ in-house platform. The flowchart of the computational program is outlined in Figure 3.

3. Evaluation of the Proposed Method

In this section, a plane strain biaxial test is presented to show the capability of the proposed method overcoming the mesh sensitivity. The classical finite element method is also employed for comparison. All simulations are conducted with the plain strain assumption.

The geometric dimensions and boundary conditions of the plane strain specimen are shown in Figure 4. The soil specimen has a width of 40 mm and a height of 80 mm. Herein, the edge friction constraint is applied to the model for stimulating the generation of shear bands; that is, the upper end is a constraint in the x-direction. The calculation parameters are listed in Table 1 by referring to []. The calculation is divided into two steps: the first step is to apply isotropic consolidation pressure σ3, and the second step is to apply deviatoric stress σ1σ3.

Two discretization schemes for FEM are employed for the purpose of comparison, where the fine mesh is composed of 3200 elements with an element size of 1 mm, whereas the element size of the coarse model becomes 2 mm. Figure 5 plots the shear band predicted by FEM with different mesh densities. The stress-strain curves of the point located in the shear band under the two cases are shown in Figure 6. It is readily apparent that the local FEM model yields nonobjective results. Upon mesh refinement, a thinner shear band and a more brittle response at the postpeak stage could be obtained. If the mesh is too fine, then the thickness of the shear band becomes zero, and energy dissipation would also tend to be zero. This nonphysical energy dissipation phenomenon is unacceptable.

With regard to the PD model, comparisons of the simulations with different cases are studied to address the following two questions: (1) can the simulation yield objective results in terms of thickness of the shear band and stress-strain relations? (2) What is the factor to determine the simulation results?

For the first question, we use the PD models, which have different grid densities but the same horizon. That is, the coarse grid model is set Δx = 2 mm and δ = 2Δx = 4 mm, while the fine grid model is set as Δx = 1.333 mm and δ = 3Δx = 4 mm. The shear band width and the stress-strain relation curve under different grid densities were compared in Figures 7 and 8(a), respectively. It can be seen from the two figures that, under the same internal scale parameter (δ), the nonlocal model can obtain a stable shear band width and stress-strain curve independent of mesh density.

For the second question, the PD models have different horizon sizes, but the grid spacing is kept fixed. That is, the small horizon model is set Δx = 2 mm and δ = 2Δx = 4 mm, while the fine grid model Δx = 2 mm and δ = 4Δx = 8 mm. It is employed to investigate the influence of the internal scale parameters (δ) on the simulation results. The width of the shear band and the stress-strain curves obtained from different models are shown in Figures 7(a), 7(c), and 8. It can be seen from that the width of the shear band is basically proportional to the nonlocal parameter δ, and the width of the shear band is roughly equal to δ. In addition, the stress-strain curves also vary at the postpeak stage with different nonlocal parameters. Therefore, in practices, the experimentally measured stress-strain curves can be used to determine the nonlocal parameters, and then one can use the calibrated parameters to investigate the strain localization failure in an objective manner.

From the above analysis, the correctness of the proposed method incorporating the strain-softening Mohr–Coulomb constitutive model is ensured, and it is confirmed that the NOSBPD model can be used as an effective regularization technique for strain localization analysis and that the localized shear band is related to the nonlocal parameter δ.

In addition, the nonlocal PD model is indeed more computationally intensive than the local FEM model under the same condition. For instance, for the FEM model with 40 × 80 elements, the computational time is 18 min, while for the PD model with 40 × 80 particles, the corresponding time is 335 min. Thus, the PD-FEM model [32] with the advantage of enhancing computational efficiency, merits consideration.

4. Progressive Failure Analysis of Soil Slope

In this section, the proposed NOSBPD method is used to analyze the progressive failure process of the soil slope with strain-softening behavior. The main departure is to illustrate the postfailure behavior of the soil slope in a progressive manner that can be effectively captured by the proposed method.

The geometry and boundary conditions of the slope model are shown in Figure 9. The slope with a height of 10.0 m is subjected to a gradually increasing surcharge load P. The calculation is also divided into two steps. First, the initial stress field is constructed by applying the gravity force, and then a uniform load is applied at the top surface of the slope, covering a length of 3 m. The arc-length method controls the gradual, increasing loading process until the slope fully fails. Material parameters are listed in Table 1. The slope is discretized into uniformly spaced particles with a grid spacing of Δx = 0.3 m. The horizon is set to δ = 2Δx.

Figure 10 shows the displacement contours and typical point (located in the shear band) stress-strain curves. It is observed that the slope with strain-softening behavior under the surcharge load has a tendency to slide outward and downward. Comparisons of the slip surface obtained by the limit equilibrium method (LEM; herein, we use the Morgenstern-Price method) using the peak and residual strengths, respectively, and the NOSBPD model are illustrated in Figure 11. It can be seen from the figure that the position of the slip surface obtained by the NOSBPD model is located between those obtained by the LEM with peak and residual strength. In Figure 11(c), interestingly, different portions of the slip surface exhibit different cohesion strengths. That is, the peak and residual strength are mobilized at the toe and top of the slope, respectively. The middle position of the slip surface near the toe is basically at the residual strength, while the strength parameters of the middle part near the top are intermediate between the peak and residual strength. It means that the slope starts to slide from the top of the slope and gradually expands downward. This nonuniform mobilization of the strength parameter indicates the occurrence of progressive failure of the soil slope.

As carried out in [14], an index referred to as the residual cohesion coefficient is defined to investigate the progressive failure process further:where cf and cr are the peak and residual cohesion strength, respectively. It characterizes the reduction of the material strength parameters of the slope with respect to the peak strength. Figure 12 shows the variation of the residual cohesion coefficient with the progressive failure process of the slope. With the occurrence and development of the plastic zone on the sliding surface of the slope, the material strength parameters gradually weaken to the residual strength.

It can be seen from the abovementioned analysis that there is stress redistribution and strength mobilization near the slip surface during the progressive failure of the slope, which cannot be assumed as a rigid body as made by the limit equilibrium method. Using the proposed method, it can avoid the dilemma that the peak strength is too dangerous and the residual strength is too conservative in the limit equilibrium method.

5. Conclusion

In this paper, the NOSBPD model is used as a regularization technique for the strain localization analysis of soil with strain-softening behavior. The Mohr–Coulomb constitutive model, by reducing the strength parameters linearly with the increase of the accumulative plastic strain, is incorporated into the NOSBPD framework to describe the strain-softening behavior of the soil. The plane strain test and the progressive failure of the soil slope are analyzed. The main conclusions obtained from the analyses are as follows:(1)The nonlocal method can remove the mesh dependence, and the NOSBPD model can be used as an effective regularization technique for strain localization analysis.(2)The localized shear band is related to the nonlocal parameter δ, which can be calibrated by the experimentally measured stress-strain curves.(3)With the occurrence and development of the plastic zone on the sliding surface of the slope, the material strength parameters gradually weaken to the residual strength. Different portions of the slope are mobilized with different strength parameters. The stress redistribution and strength mobilization along the slip surface render the assumption made by the limit equilibrium method incorrect. That is, it is not a rigid body.(4)The proposed method involves being able to capture strength mobilization in different parts of the sliding surface during the progressive failure process of the soil slope. It can avoid the dilemma that the peak strength is too dangerous and the residual strength is too conservative due to the rigid body assumption in the limit equilibrium method.

As a nonlocal theory, NOSBPD is much more computationally intensive than conventional FEM. There is a necessity to apply the FEM-NOSBPD coupled method [32] for soil strain localization analysis. In addition, there is a need to remove the constraint of the small deformation assumption for simulating the progressive failure of slope in the large deformation regime. Also, considering the interaction of soil and water, incorporating a more advanced constitutive model of soil will also require further research.

Data Availability

The datasets generated and analyzed during the current study are available from the corresponding author upon reasonable request.

Conflicts of Interest

The authors declare that they have no conflicts of interest.