A simple method to reduce mesh dependency in modelling landslides involving brittle soils

Deformation processes of the slopes are usually categorised into four distinct phases (Leroueil, 2001): pre-failure, failure, post-failure and eventual reactivation. Soil brittleness undoubtedly affects the pre-failure phase due to the occurrence of a progressive failure (Bishop, 1967), but it could also play an important role in the kinematics of the unstable soil mass during the post-failure phase (Leroueil, 2001). Deformation and failure mechanisms of the slopes can be effectively analysed using numerical methods capable of simulating both the progressive development of shear zones within the slope and the run-out process of the landslide. In this paper, the material point method (MPM) is used for analysing the deformation processes that occur in slopes consisting of soils with strain-softening behaviour (Yerro et al., 2016; Fern et al., 2019; Conte et al., 2020; Alonso, 2021). In these circumstances, the solution is often affected by a lack of convergence and a strong dependency on the employed mesh (Soga et al., 2016). This problem can be overcome using the so-called regularisation techniques. However, these techniques are generally very expensive from a computational viewpoint. A simple-to-use procedure based on the results of conventional direct shear tests, is proposed herein to reduce the numerical problems due to mesh dependency. This methodology is completely analytical and requires few parameters with a clear physical meaning. A case study is analysed to show the effectiveness of the proposed approach.

An original procedure is presented to reduce the mesh dependency when dealing with brittle soils. This procedure utilises, in the context of a strain-softening Mohr–Coulomb model, a strength reduction law similar to that proposed by Yerro et al. (2016) in which the cohesion intercept c′ and the angle of shearing resistance φ′ of the soil reduce with an increase of the deviatoric plastic strain invariant $εdp$⁠, according to the following equations:

where

$eijp$ is the deviatoric part of the plastic strain tensor, and the model parameter λ controls the strength reduction rate. The proposed procedure is analytical and is based on the results of direct shear tests. Specifically, a complete sequence of shearing cycles necessary to attain the residual strength of the involved soil is required. Figure 1 shows an example of these data. In this figure, δ_p is the horizontal displacement when peak strength occurs, δ_pp is the horizontal displacement (net of δ_p) corresponding to the post-peak condition (in the first shearing cycle) and δ_r is the net horizontal displacement at which the shear strength attains the residual condition. Moreover, τ_p, τ_pp and τ_r are the peak, post-peak and residual shear strength of the soil, respectively (Calabresi & Manfredini, 1973).

In the proposed method, the following equation is assumed to relate the soil shear strength to the plastic shear strain γ^p:

The plastic shear strain is expressed as $γp=δ−δp/l$⁠, where l is the thickness of the shear zone and δ is the total horizontal displacement. As a consequence, equation (4) takes the form:

with $δ>δp$⁠.

Considering the first shearing cycle (Fig. 1), the ratio λ/l can be calculated by replacing $τ=τpp$ and $δ−δp=δpp$ into equation (5). The resulting expression for λ/l is

where the values of τ_p, τ_pp, τ_r and δ_pp are evaluated from the τ–δ curve obtained experimentally (Fig. 1). It is worth noting that equation (6) provides directly the ratio λ/l rather than the parameters λ and l, which are difficult to be evaluated separately. Once λ/l is calculated using equation (6), it needs to check whether equation (5) approaches the residual value τ_r when $δ−δp=δr$⁠. To this end, the ratio β may be considered:

in which $τr∗$ denotes the shear strength calculated using equation (5) with $δ−δp=δr$⁠. The resulting value of β should be close to zero to ensure a good agreement between experimental and analytical curves. Otherwise, new values of τ_pp and δ_pp are assumed and the above-described steps are again performed until β results close to zero. Referring to the experimental data shown in Fig. 1, a ratio λ/l = 300 m⁻¹ is evaluated using this iterative procedure. The τ–δ relationship provided by equation (5) along with this value of λ/l is plotted in blue in the same figure. The resulting relationship approximates the experimental data very well.

It is worth noting that the thickness l of the shear zone in a direct shear test is related to the soil grain size, whereas in the numerical simulations this thickness depends on the mesh size, L (Fig. 2). Likewise, the model parameter λ from the direct shear tests may be different from that adopted in the calculations. This latter is herein denoted λ* to distinguish it from the former. An evaluation of λ* is performed in this study by imposing the same τ–δ relationship (from peak to residual) obtained from the direct shear tests, to the soil of the shear zone in the numerical simulation (Fig. 2). Under this assumption, equation (5) takes the form:

Consequently, equation (6) expresses both λ/l for the direct shear test, and λ*/L for the numerical simulation. To support this assumption, some numerical analyses concerning an ideal shear test under plane strain conditions are performed using the finite-element code, Tochnog Professional (Roddeman, 2022). It is assumed that the tested soil sample consists of a material characterised by the experimental τ–δ relationship shown in Fig. 1 with λ/l = 300 m⁻¹. The shear zone has a thickness equal to L. Meshes with a different size are employed to discretise this shear zone. Specifically, the following values of L are considered: 1, 10, 50 and 100 cm. A comparison between experimental and predicted τ–δ relationships (from peak to residual) is presented in Fig. 3 for different values of λ* and L, but with the same ratio λ*/L = 300 m⁻¹ coincident with λ/l deduced experimentally (Fig. 1). A very close agreement between simulation results and experimental data is obtained under these conditions (λ*/L = λ/l, for a given L). In other words, the model parameter λ* depends on the mesh size L and is related to the ratio λ/l deduced from the direct shear tests.

Summarising, the following simple procedure is proposed to reduce the influence of the mesh dependency on the slope response when soils with strain-softening behaviour are involved. Ratio λ/l is first evaluated from the available data of direct shear tests using equation (6) along with the iterative procedure previously described. Once this ratio is known, the model parameter λ* to be used in the numerical simulation in which a mesh with a given size L is employed, can be evaluated as follows:

This approach was recently used in conjunction with MPM to simulate successfully the deformation processes of a landslide in brittle soils (Troncone et al., 2022).

To assess the capability of the proposed procedure to reduce the mesh-dependency influence on the numerical results, an ideal slope constituted by a brittle soil with different strength reduction rates is considered. This slope is characterised by a height of 15 m and an inclination angle of 45° (Fig. 4), and consists of a homogeneous dry soil whose properties are indicated in Table 1, in which γ is the unit weight, E′ is the Young modulus, ν′ is the Poisson ratio and ψ is the dilatancy angle. Figure 5 shows three τ–δ curves obtained from direct shear tests. These curves have the same peak strength and the same residual strength, but a different strength reduction rate with values of λ/l equal to 100, 300 and 600 m⁻¹, respectively. As shown in Fig. 5, the higher the value of λ/l, the faster the strength reduction from peak to residual. To trigger a deformation process in the slope, an excavation 5 m deep is carried out at the slope toe, as indicated in Fig. 4. For the sake of completeness, the safety factor of the slope is first evaluated using the limit equilibrium method (Morgenstern & Price, 1965). The resulting values are 3·73 and 0·75 when the peak shear strength parameters and the residual ones are used in the calculations, respectively. As a consequence, a progressive failure may occur due to the excavation (Bjerrum, 1967). The deformation processes occurring in the slope (from the pre-failure phase to the post-failure one) are simulated using the MPM code Anura3D (www.anura3d.com). In this context, different meshes with L = 0·5, 1 and 1·5 m are considered. Each mesh is made up of triangular elements and any active element contains initially three material points. Both vertical and horizontal displacements are constrained at the bottom of the domain and the horizontal displacements are prevented on the vertical sides. The initial stress state of the slope is reproduced using the well-known gravity loading procedure. Afterwards, the excavation is simulated by removing the material points located in the area indicated in Fig. 4. The values of the model parameter λ* to be used in the simulation are evaluated as a function of λ/l and L using the procedure described in the previous section. The resulting values of λ* are listed in Table 2.

Some simulation results are documented in Figs. 6–8, where for each value of λ/l a comparison in terms of the final displacement field calculated using a different mesh size is presented. These results show that the slope response depends on the value of λ/l expressing the strength reduction rate in the direct shear test. Specifically, when λ/l = 100 m⁻¹, the soil strength reduction is slow and the slope only undergoes small displacements (Fig. 6). In other words, the excavation performed at the slope toe is not sufficient to reduce the soil strength to the point of causing a slope failure. On the contrary, the slope collapses when the higher values of λ/L are assumed – that is, 300 m⁻¹ (Fig. 7) and 600 m⁻¹ (Fig. 8). In these circumstances, a rapid reduction of the soil strength occurs due to the excavation up to causing a slope failure. This result is found for all employed meshes that are characterised by the same ratio λ*/L assumed equal to the associated λ/l, according to the proposed procedure. In Figs. 7 and 8, the run-out distance (defined as the distance between the tip of the displaced material and the vertical wall of the excavation) is also indicated. As can be observed from the presented results, the final configuration of the displaced material, the displacement magnitude and the run-out distance are slightly affected by the mesh size, provided that the same value of λ*/L (with λ*/L= λ/l) is employed in the simulation. The greatest discrepancy found in terms of soil displacement for the different values of λ*/L considered, is about 10%. In addition, the values of the run-out distance differ from each other of a few decimetres. It is also worthy to note that in Figs. 6(c) and 7(a) the value of λ* = 150 is used with a different mesh size (L = 1·5 m and L = 0·5 m, respectively) and consequently with a different value of λ*/L (100 and 300 m⁻¹, respectively). As can be seen by comparing Figs. 6(c) and 7(a), very different results are obtained when λ* and L are arbitrarily chosen without accounting for the ratio λ/l obtained experimentally. Summarising, shape factor λ* is mesh dependent, but this dependency can be reduced by assuming a constant value of λ*/L that is equal to λ/l obtained from the τ–δ relationship of a direct shear test using the procedure proposed in the present study. This assertion is also corroborated by the results shown in Figs. 9(a)–9(c), which document the time–total displacement relationship at point P, whose initial location is indicated in Fig. 4, for the different mesh sizes and the different values of λ*/L considered.

Finally, Figs. 10–12 show the calculated deviatoric strain distribution. Consistently with the results in terms of displacements (Figs. 6–8), the induced deformations are low when λ/l = 100 m⁻¹ (Fig. 10), whereas a slope failure occurs when the higher values of λ/l are assumed (Figs. 11 and 12). In these latter circumstances, the shape and location of the failure surface are very similar when the same value of λ*/L is adopted in the analysis. Although the strain magnitude is different, it is important that the product $λ∗εdp$ is mesh-independent because it appears in equations (1) and (2) defining the softening law of the soil. In this connection, Figs. 13(a)–13(c) show the time evolution of $λ∗εdp$ calculated in the material point Q for the considered values of λ* and L. Point Q is located in Fig. 4 and falls within the shear zone, where large deviatoric strains occur when λ*/L = 300  and 600 m⁻¹. As can be seen,  $λ∗εdp$ is essentially unaffected by the mesh size, confirming the effectiveness of the proposed method.

The proposed procedure appears quite attractive for practical purposes because it is simple to use and requires few parameters as input data. However, some aspects concerning its application to real case studies should be clarified. Although the present method is herein used in conjunction with MPM, it could also be employed with different mesh-dependent numerical techniques in which the softening rule considered in the present study is implemented. In addition, the proposed method is based on the results of the direct shear test that is simple and very widespread in engineering applications. However, it is often criticised for the non-uniformity of stress and strain applied to the sample due to the boundary conditions imposed by the box. Potts et al. (1987) analysed in detail this aspect and showed that the above-mentioned non-uniformity slightly affects the load–displacement behaviour of the soil. Finally, considering that the model parameter λ* depends on L, a uniform mesh should be used for homogeneous soils. However, this could be very expensive from a computational viewpoint, especially when the slope is particularly long. To overcome this issue, the layers of homogenous soil may be split into some smaller regions (clusters) with a different size L but with the same value of λ*/L.

A novel methodology has been proposed to reduce the numerical problems due to mesh dependency, which arise when soils with strain-softening behaviour are involved in the deformation processes of slopes. The method is completely analytical and simple to use. It is based on the experimental data from direct shear tests inclusive of the shearing cycles necessary to reduce the soil strength from peak to residual. In addition, a negative exponential law is employed to match the experimental τ–δ curves provided by these tests, and an iterative and quickly converging procedure is developed to evaluate the ratio of the model parameter λ controlling the strength reduction rate, to the thickness of the shear zone in the direct shear tests, l. Once λ/l is known, the model parameter λ* to be used in the numerical simulation can be determined for a given mesh size L, using a simple equation (equation (9)). The numerical simulations performed in the present study using MPM corroborate the reliability of the proposed method.

All data used are available from the corresponding author by request.

$c′$

cohesion intercept

$c′p$

peak value of cohesion intercept

$c′r$

residual value of cohesion intercept

$E′$

Young's modulus

$eijp$

deviatoric part of the plastic strain tensor

$L$

thickness of the shear zone in numerical simulation

$l$

thickness of the shear zone in direct shear test

$β$

ratio defined by equation (7)

$γ$

unit weight

$γp$

plastic shear strain in the shear zone

$δ$

horizontal displacement in direct shear test

$δp$

horizontal displacement corresponding to the peak strength in direct shear test

$δpp$

horizontal displacement (net of $δp$⁠) corresponding to the post-peak condition in direct shear test

$δr$

horizontal displacement (net of $δp$⁠) corresponding to the residual condition in direct shear test

$εdp$

deviatoric plastic strain invariant

$λ$

model parameter controlling the strength reduction rate in direct shear test

$λ∗$

model parameter controlling the strength reduction rate of the shear zone in numerical simulation

$ν′$

Poisson's ratio

$τ$

shear stress

$τp$

peak shear strength

$τpp$

post-peak shear strength

$τr$

residual shear strength

$φ′$

angle of shearing resistance

$φ′p$

peak angle of shearing resistance

$φ′r$

residual angle of shearing resistance

$ψ$

dilatancy angle