Contents lists available at ScienceDirect

Computers and Geotechnics

journal homepage: www.elsevier.com/locate/compgeo

Research Paper

On the use of nonlocal regularisation in slope stability problems

CrossMark

F.C. Summersgill a,1 S. Kontoeb'*, D.M. Pottsc

a Tony Gee & Partners LLP, Hardy House,140 High Street, Esher, Surrey KT10 9QJ, UK b Department of Civil & Environmental Engineering, Imperial College London, SW7 2AZ, UK

c GCG Professor of Geotechnical Engineering, Department of Civil & Environmental Engineering, Imperial College London, SW7 2AZ, UK

ARTICLE INFO

Article history: Received 9 June 2016

Received in revised form 9 September 2016 Accepted 16 October 2016 Available online 24 October 2016

Keywords:

Nonlocal regularisation Slope stability analysis Coupled consolidation

ABSTRACT

This study examines the use of nonlocal regularisation in a coupled consolidation problem of an excavated slope in a strain softening material. The nonlocal model reduces significantly the mesh dependency of cut slope analyses for a range of mesh layouts and element sizes in comparison to the conventional local approach. The nonlocal analyses are not entirely mesh independent, but the predicted response is much more consistent compared to the one predicted by local analyses. Additional Factor of Safety analyses show that for drained conditions the nonlocal regularisation eliminates the mesh dependence shown by the conventional local model.

© 2016 The Authors. Published by Elsevier Ltd. This is an open access article under the CC BY license (http://

creativecommons.org/licenses/by/4.0/).

1. Introduction

The modelling of slope failure in a strain softening material with conventional Finite Element models can be particular challenging. In a typical finite element analysis the strain developed along the shear bands is calculated using the displacement information computed at the nodes of the elements. This affects both the shear band thickness and the direction of its development [1]. The calculated strain is used to assess the degree of softening experienced by the material at that point. There can be a large difference, i.e. a high gradient, between displacement and therefore strain at neighbouring points. The potential strain concentration at one single point can lead to convergence problems and ambiguities in the development of the slip surface. In addition, the sizes of the elements in the mesh restrict the minimum size of the shear band to the distance between two points of known displacement, i.e. two nodes [2]. For example, if 8-noded two-dimensional elements are employed, the minimum shear band thickness is restricted to the width of half an element. These inherent limitations of the finite element method render the solution to be mesh dependent.

Several approaches have been proposed in the literature to try and regularize the numerical solution and model rigorously the shear zone and these can be divided into three main categories; Cosserat theories [3], gradient theories [3-6] and nonlocal

* Corresponding author. E-mail addresses: freya.summersgill@tonygee.com (F.C. Summersgill), stavroula. kontoe@imperial.ac.uk (S. Kontoe), d.potts@imperial.ac.uk (D.M. Potts). 1 Formerly Imperial College London.

approaches [1,3,7,8]. The present study focuses on the latter approach which is particularly attractive because it does not alter the fundamental governing equations, but it does introduce the calculation of a nonlocal strain as a variable by spatially averaging the local strains [9]. This makes the approach of nonlocal strain regularisation more straightforward to implement in an existing finite element code, compared to Cosserrat and gradient theories.

The non-local method [10,11] adopts a distribution function which spreads the strain of the material at a point over a predefined surrounding volume. The local method calculates the extent of strain softening with reference to the strain at that point alone. To define the contribution of nonlocal strains to the yielding of the material requires the additional input of a characteristic length parameter, which controls the contribution of local strains to the nonlocal calculation depending on the distance of the local strains from the calculation point. In the present study the nonlocal model of Galavi and Schweiger [1] (G&S) is employed in a coupled consolidation problem of an excavated slope in a strain softening material aiming to investigate the performance of this nonlocal approach in terms of mesh dependence and computational cost. Furthermore the impact of the two nonlocal parameters, the defined length (DL) and the radius of influence (RI) on the numerical predictions is thoroughly investigated providing guidance for their use in boundary value problems. The defined length is an integral parameter of the G&S method which modifies the rate of softening, while the radius of influence is an optional parameter which makes the computation more efficient by reducing the number of local strains referenced in the nonlocal strain calculation.

http://dx.doi.org/10.1016/j.compgeo.2016.10.016 0266-352X/© 2016 The Authors. Published by Elsevier Ltd.

This is an open access article under the CC BY license (http://creativecommons.org/licenses/by/4.0/).

2. Mesh dependency of local strain softening slope stability analyses

As previously discussed, the modelling of slope failure in a strain softening material with conventional Finite Element models can be particular challenging, as the solution can be very sensitive to the adopted mesh discretization. For an assumed problem geometry, different element layouts can be employed and the mesh is usually refined around the areas of potential strain concentration leading to varying element sizes within the FE model. The sensitivity of a local strain softening model on these two aspects (i.e. the element size and the elements' layout) is first demonstrated, before exploring the use of nonlocal regularisation in strain softening materials. A parametric study employing both biaxial compression analyses,

as well as slope stability analyses was carried out for this purpose. The biaxial compression analyses first examine the role of element size within a uniform mesh discretisation, while in the slope stability analyses both aspects of mesh discretisation are investigated. All the analyses presented in this paper were carried out in plane strain with the Finite Element code ICFEP (Imperial College Finite Element Program) [12], using a strain softening variant of the MohrCoulomb model [13]. This is an elasto-plastic model in which softening behaviour is facilitated through a variation of the angle of shearing resistance u', and the cohesion intercept c' with the devi-atoric plastic strain invariant, as shown in Fig. 1. The limits for peak (up', cp) and residual (ur', cr) strength are specified in the model by a percentage value of the deviatoric plastic strain invariant (Epdp, Epdr, respectively) which is defined in Eq. (1):

<Pp.Cp

(p ftC r

1 -!-w

Eir El

Fig. 1. Variation of the angle of shearing resistance u' and the cohesion intercept c' with the deviatoric plastic invariant Epd adapted from Potts et al. [13].

(a) (b)

Fig. 2. Boundary conditions for biaxial compression.

Table 1

Summary of biaxial compression analysis for local and nonlocal strain softening models.

Mesh identification & element arrangement Element size (m) Local model:strain rate peak to residual (%) Nonlocal model:DL values 0.525 m 1.05 m Ratio of DL over element size 2.1 m

10 x 10 2.1 5-20 - - 1:1

20 x 20 1.05 10-40 - 1:1 1:2

40 x 40 0.525 20-80 1:1 1:2 1:4

Ed (£Pl " £P2)2 + («p2 - ep3)2 + (ep3 - «pi)' (1)

where ep1, sp2, ep3 are the principal plastic strains.

Eight-noded isoparametric elements with reduced integration were used and an accelerated modified Newton-Raphson scheme with a sub-stepping stress point algorithm was employed to solve the nonlinear finite element equations [12].

2.1. Biaxial compression analyses

The biaxial compression analyses, compress a quadrilateral mesh from two opposing sides, whilst leaving the two other opposing sides free to deform. The two lines of symmetry of the problem permit a quarter of it to be analysed, as shown in Fig. 2 together with the adopted boundary conditions. The initial conditions consist of a vertical stress of 50 kPa imposed on the top horizontal boundary and a 100 kPa horizontal stress applied along the right hand side lateral boundary. Restricting the horizontal displacement of the top mesh boundary in Fig. 2(b) and imposing an equal vertical displacement to all nodes along this boundary was sufficient for a slip surface to form.

London Clay was chosen as a representative strain softening material. Kovacevic [14] and Potts et al. [15] employed local strain softening analyses to model London Clay cuttings with the peak

Table 2

Model parameters for biaxial compression and cutting slope analyses [14,15].

Problem analysed Property Assumed value

Biaxial compression & slope Bulk Unit weight, c 18.8 kN/m3

stability Peak strength (bulk) cp = 7 kPa, up = 20°

Residual strength cr = 2 kPa, ur =13°

Poisson's ratio, i 0.2

Angle of dilation, W 0°

Young's modulus, E 25(p' + 100) (min

4000 kPa)

Biaxial compression Plastic deviatoric Varied

strain, Ed

Slope stability Plastic deviatoric Peak 5%, residual 20%

strain, Ed

Coefficient of Ik) = 5 x 10~10 m/s

permeability, k b = 0.003 m2/kN

strength properties retained up to a plastic deviatoric strain of 5% and then linearly reducing to their residual values which are reached at a plastic deviatoric strain of 20%. This produced a realistic failure time and mechanism that agreed with field data for cutting slopes of the same dimensions in London Clay. It should be noted that due to the mesh dependency of the local method, these strain limits are only appropriate for the element sizes employed in the original analyses of Kovacevic [14] and Potts et al. [15]. Their adopted softening rate and the associated strain softening limits were based on comparisons of experimental data with simple shear finite element simulations. The displacements for the limits of softening obtained in the simple shear analyses were related to the size of the elements used in the mesh for the cutting analyses in Potts et al. [15], with 2 m by 1m elements beneath the slope surface.

Based on the work of Kovacevic [14] and Potts et al. [15], in the present study the 5% and 20% strain limits for a local strain softening analysis are employed with a mesh with elements sized 2.1 m by 2.1 m (10 x 10 elements), to give an appropriate softening rate for London Clay. Three square mesh layouts of the form shown in Fig. 2 were considered, with 10 x 10, 20 x 20 and 40 x 40 elements as summarised in Table 1. However, the mesh is 21 m by 21m in size to produce element sizes of 2.1m, 1.05 m and 0.525 m. In order to maintain a consistent softening rate in all the local analyses, the strain limits are scaled so that for each considered element size the softening limits are reached for the same displacement. With smaller element sizes, the softening limits for the mesh dependent local percentage strain softening model must be varied to produce the same softening rate, by maintaining the same displacement at which the softening limits are reached. The local strain softening model can be manipulated in this way, but only for this type of analysis with a mesh that has a regular square format. Therefore, for the 20 x 20 mesh with elements half the size, twice the limits are specified as 10% and 40% and similarly for the 40 x 40 mesh, 20% and 80% strain limits are specified. The biaxial compression analyses are drained and the adopted properties for the London Clay are given in Table 2.

Fig. 3 presents the reaction load on the top of the mesh given the applied vertical displacement for all the biaxial compression analyses. As expected, the local analyses with the altered softening limits for the different element sizes produce coinciding results. This confirms that the softening rate in the local strain softening

- Edp=5% Edr=20%, mesh 10x10

H — - Edp=2t % Edr=80%, m esh 40x40

0 -0.2 -0.4 -0.6 -0.8 -1 -1.2 -1.4 -1.6 -1.8

Vertical displacement at top boundary (m)

Fig. 3. Biaxial compression analyses with local strain softening employing 21 m by 21 m meshes with varying numbers of elements.

Fig. 4. Index of meshes for excavated slope analyses, arranged by mesh layout and element width. Measurements are specified for elements beneath excavated slope area. The area to be excavated is shaded grey.

Fig. 5. Typical finite element mesh and boundary conditions.

model is element size dependant and that in problems with uniform mesh discretisation one can achieve the desired softening rate by adjusting the element size and strain softening limits appropriately.

2.2. Slope stability analyses

In order to investigate the mesh dependency of the local strain softening model predictions in cases with varying element sizes and layouts, a set of eleven meshes were considered, as shown in Fig. 4. These meshes are all designed to simulate an excavated slope 10 m high with a 1 in 3 gradient (vertical to horizontal) and three of the considered layouts can be employed to simulate vertical stabilisation piles, Fig. 4(a)-(h), which is a common stabilisation scenario in slope stability FE analysis. It should be noted that in all the presented cutting analyses a coupled consolidation formulation is adopted unless otherwise stated, which allows the modelling of the generation of excess pore fluid pressures during excavation and their subsequent equilibration with time.

The same soil properties for London Clay as specified in Table 2 are considered for the cutting analyses, but with a variation of Young's modulus, E with confining pressure p' (instead of the constant value of E which was used in biaxial compression analyses). In addition, only one strain rate was specified with the peak strength properties being retained up to a plastic deviatoric strain of 5% and then linearly reducing to their residual values, which are reached at a plastic deviatoric strain of 20%. A typical FE mesh with the slope dimensions and the associated mechanical and hydraulic boundary conditions is shown in Fig. 5. No horizontal displacement was allowed on the vertical boundaries, whereas the bottom boundary was fixed in both horizontal and vertical directions. Before the excavation of the slope, initial stresses are specified in the soil using a bulk unit weight of y = 18.8 kN/m3 and a uniform coefficient of lateral earth pressure, K0 = 2.0. The pore water pressures are hydro-static with 10 kPa suction specified at the ground surface, following the average height expected for the phreatic surface in the UK [16]. The bottom and side boundaries are impermeable. The permeability, k of the soil is modelled as isotropic and is

Fig. 6. Change in horizontal displacement with time for cutting analyses employing the local strain softening model and various meshes.

6 5 4 3 2 1 0 -1 -2 -3 -4 -5 -6

Distance from calculation point (m) Fig. 7. The impact of DL on the weighting function distribution.

linked to the mean effective stress, p' using the non-linear relationship in Eq. (2) [17].

k = koe-bp' (2)

The slope was excavated in horizontal layers over 0.25 years. This unloads the soil surrounding the excavation and the low permeability of the soil creates negative pore water pressures. After excavation 10kPa suction is applied at the free boundary. Time and consolidation allow these excess pore water pressures to slowly dissipate. The changes in pore water pressures and strain softening behaviour of the stiff clay eventually lead to failure of the slope. The point of failure is defined as the last increment of

the analysis that will converge with a time step of 0.01 years. Initially time steps of 1 year are employed and the size of the incremental step is reduced as slope failure is approached.

The computed response employing the local strain softening model is presented in Fig. 6 in terms of horizontal displacement variation with time for the mid-slope surface point. The results are shown to be extremely mesh dependent both in terms of the time to failure and the movement of the slope between excavation and failure. The range of time to failure varies from 14 years for the inclined elements mesh with 1 m wide slope elements, Fig. 4(j), to a slope that did not show signs of failure 250 years after excavation when employing the base pile mesh layout and 0.525 m wide slope elements, Fig. 4(f). In contrast, the same base pile mesh layout with 1.05 m wide slope elements, Fig. 4(e) reached failure in less than 30 years. Two of the analyses employing the mid-slope pile mesh layout produce similar results, with 1.05 m wide elements, Fig. 4 (b) and 0.525 m elements, Fig. 4(c). Slope failure for these meshes occurred after 140 years and 160 years respectively (failure points are not shown in the graph). However, the 2.1 m wide elements for this same layout, took only 40 years to reach failure.

3. Mesh dependency of non-local strain softening slope stability analyses

The biaxial compression and cutting analyses presented so far demonstrate clearly the significant mesh dependency of conventional local strain softening models. On the other hand, non-local models have been implemented in several numerical codes as a regularisation tool for slip surface or fracture formation. Summersgill [18] considered three nonlocal approaches, the original formulation of Eringen [7], the Over-nonlocal of Vermeer and Brinkgreve [8] and the G&S of Galavi and Schweiger [1], in biaxial compression simulations showing that the G&S suffered the least mesh dependency. However, the slip surface developed in biaxial compression analyses are kinematically constrained and therefore the conclusions of such analyses are not directly applicable to any boundary value problem. In the present study, the mesh dependency of the G&S method is examined not only in biaxial compression, but also in the coupled consolidation problem of the excavated slope. The latter is an example of an analysis which does not predetermine the location of the slip surface and therefore provides no kinematic restraint on slip surface development.

3.1. The non-local G&S method

3.2. Comparison of nonlocal and local strain softening analyses

A fully nonlocal model treats both stress and strain as nonlocal components. However, when the model is employed as a regularisation tool, it is common to adopt a partial nonlocal constitutive relation in which only the nonlocal plastic strains control the softening [1]. One of the first formulations for nonlocal plastic strains is presented in Eq. (3). This was proposed by Eringen [7] for strain hardening applications and by Bazant et al. [19] in a strain softening damage model.

ep'(xn) (x(x'n)ep(x'n))dx'idx'2dx3 (3)

where ep is the = accumulated plastic deviatoric strain, * denotes the nonlocal parameter, xn is the global coordinate at which the calculation of the nonlocal plastic strain, ep/ is required, whereas xn refers to all the surrounding locations, i.e. the location of reference strain, with n = 1, 2, 3. Therefore, ep(xn) equals the reference (i.e. local) strain at the reference location. The weighting function, ra (xn) is defined for all the reference locations, but it is centred at the location xn. The weighting function (see Eq. (4)) in the case of the G&S method has been designed to limit the central concentration of strains, assuming that the development of the slip surface is influenced by the directly surrounding areas and not by the concentrated strain at the centre of the slip surface.

\j(X'n - Xnf (X'n - Xn)

x(x'n ) -—-exp

(x'n xn) (xn xn)

The defined length, DL, influences the distribution of the weighting function as shown in Fig. 7. A higher DL is expected to result in a wider slip surface with a lower maximum nonlocal strain. This affects the rate of strain softening, as the same input of local strain will result in a lower nonlocal strain for a larger defined length and therefore a slower rate of softening. The integral of the weighting function in the three dimensions x1, x2 and x3 is referred to as the reference volume, Vx, as shown in Eq. (5). This is used to normalise the calculation of the nonlocal strain and for the G&S distribution is equal to approximately one. The latter attribute of the function ensures that a uniform field of strain would remain unmodified.

^x(xn)dx1dx2dx3

In this section both the biaxial compression and the cutting analyses are repeated with the G&S model to investigate to what extent nonlocal regularisation reduces the mesh dependency.

3.2.1. Biaxial compression analyses

The biaxial compression analyses were repeated with the G&S model considering a range of DL values as summarised in Table 1. The appropriate range of DL values to investigate is influenced by the element sizes in the meshes. A suitable DL value should be at least equal to or larger than the element size, to include sufficient strain reference points. However, it will be demonstrated that subject to the above restriction for a given DL, the softening rate is not influenced by element size and therefore the element sizes of the mesh employed for the biaxial compression analyses can be of any size. The second nonlocal parameter the radius of influence, RI, is taken as three times the DL. RI is an optional parameter that limits the distance from the calculation point of the neighbouring strains that will be included in the nonlocal strain calculation and therefore reduces the computational cost. In Fig. 7, it can be seen that at a distance greater than 3 times DL from the calculation point, the contribution of local strain to the nonlocal strain calculation is extremely small. The selected value of RI is also justified based on the parametric investigation of Summersgill [18], who for biaxial compression showed that an RI equal to three times the defined length, DL, achieves the same accuracy as an analysis with no RI specified, i.e. where all neighbouring strains contribute to the calculation of a nonlocal strain no matter how far away they are from the position for which the nonlocal strain is being evaluated. The role of RI, is investigated in detail in the last part of this study.

Fig. 8 presents the reaction load on the top of the mesh for the applied vertical displacement for all the biaxial compression analyses (local and non-local). As expected, the value of DL influences significantly the softening rate, but for a given value of DL the softening rate is constant irrespective of the number or size of elements. These results demonstrate the low mesh dependency of the nonlocal method, but also the importance of selecting an appropriate value of DL which matches the target softening rate. The DL = 1.05 m analyses produce a close match to local strain softening results and therefore a suitable DL to represent the softening rate of London Clay. It should be clarified that the softening rate

Fig. 8. Biaxial compression analyses with local and nonlocal strain softening employing 21 m by 21 m meshes with varying numbers of elements.

Fig. 9. Change in horizontal displacement with time for cutting analyses employing the nonlocal strain softening method with DL = 2.1 m and various meshes.

Fig. 10. Mesh Index for 6 multiple pile layout meshes with a variation of element sizes as defined in Table 3 for an excavated slope 10 m high with a 1 in 3 (vertical to horizontal) slope and 16 m wide base of excavation.

Table 3

Variation of element sizes and number for the multiple pile layout meshes with different discretisations.

Mesh identification Element size (width to height) Below base of excavation Below slope of excavation Number of elements Reference in Fig. 10

Largest elements (as in Fig. 4g) 1m x 1m 1.2 m x 1 m 1460 (a)

Smallest elements (as in Fig. 4h) 0.5 m x 0.5 m 0.6 m x 0.5 m 5840 (b)

Thin below base 0.5 m x 1 m 1.2 m x 1 m 2180 (c)

Short with thin below base 0.5 m x 0.5 m 1.2 m x 0.5 m 3560 (d)

Thin below slope 1m x 1m 0.6 m x 1 m 2360 (e)

Short with thin below slope 1 m x 0.5 m 0.6 m x 0.5 m 3920 (f)

adopted in the local analyses is considered as the benchmark rate, because is consistent with the simple shear tests simulations of Potts et al. [15]. As previously discussed, Potts et al. [15] showed that this softening rate results in a realistic prediction of the time to failure and failure mechanism that agreed with field data for 10 m high cutting slopes with a 1 in 3 gradient in London Clay. It should be also noted that the selection of DL in granular materials is normally related to the mean grain size. For example, Galavi and Schweiger [1] suggest DL to be 10-20 times the mean grain diameter in dense sands. For boundary value problems in FE and in particular in clay materials such a correlation with the mean grain size becomes impractical, as the element size in any direction should be at most equal to DL for the nonlocal calculation to be effective. Therefore in agreement with Galavi and Schweiger [1] the selection of DL can be based on the softening rate although this might result in unrealistically large shear band thickness.

Another interesting point of the biaxial compression results is that for the same large value of imposed vertical displacement and after material softening has finished for the nonlocal analyses (see for example the case of DL = 1.05), the local results exhibit a higher reaction load. This reaction load for the local results continues to decrease long after the residual plateau in Fig. 8 has been reached by the nonlocal analyses. However, the continued reduction in reaction load occurs at a much slower rate creating a shallower gradient than the initial softening gradient for the local results. These observations are presumably related to the inherent differences in the way the softening parameter (i.e. the deviatoric plastic strain) is calculated in the two approaches.

3.2.2. Slope stability analyses

The cutting analyses presented previously were repeated with the nonlocal G&S method, using the same strain softening limits. In this set of analyses, the nonlocal G&S method employs a defined length, DL of 2.1 m with a radius of influence, RI of 6.3 m (i.e. as previously RI is equal to three times the DL). It should be noted that the minimum value employed for DL should be equal to the

maximum element size found in the meshes. Therefore a DL of 2.1 m, is the smallest value that could be used given the size of the elements in the meshes employed (see Fig. 4). These analyses seek to evaluate the mesh dependence of the nonlocal method. The selection of an appropriate DL and therefore time to failure is addressed subsequently.

The computed response for the nonlocal model is presented in Fig. 9 in terms of horizontal displacement variation with time for the mid-slope surface point. In contrast to the corresponding local results (Fig. 6), the nonlocal analyses demonstrate a more consistent response for the various considered meshes. The range of time to failure is only 40 years, between 61 and 97 years after excavation, as opposed to a range of more than 250 years for the local analyses. Excluding the inclined element meshes and the multiple location pile mesh with 0.5 m elements, the range of time to failure reduces to only 10 years. The increase in horizontal displacement of the mid-slope with time is more gradual and continuous without the sudden changes in the speed of slope movement exhibited by some local analyses in Fig. 6. These sudden changes in displacements are particularly evident for the local analyses employing the mid-slope pile layout with 2.1 m width slope elements, Fig. 4(a), and the inclined elements mesh with 0.5 m height elements, Fig. 4(k). In the first 5 years the nonlocal results are very similar, but the local displacement analyses have already diverged. Overall there is a small variation in the horizontal displacement of mid-slope over time, although the shape produced in Fig. 9 is similar for all of the nonlocal strain softening analyses. These analyses also result in a similar horizontal displacement value at failure, indicating that the failure mechanism and the location of the slip surface should be similar.

Nevertheless, the nonlocal method has shown some mesh dependency for the four distinct mesh layouts that were considered in Fig. 4, in terms of predicting the time to failure and movements. The results for the same mesh layout with a variation in element size provided consistent results, except for the multiple pile layout meshes. The element sizes employed in the two meshes

(e) Thin below slope

Fig. 11. Accumulated local plastic deviatoric strain contours illustrating the effect of element discretisation on the position of the slip surface for six nonlocal analyses for the last converged increment.

of that layout are 1.2 m x 1 m and 0.6 m by 0.5 m for Fig. 4 (g) and (h) respectively. To further investigate the mesh dependency of the local and nonlocal G&S methods in terms of number of elements, the same multiple pile layout mesh is used with a variation in element size and distribution to apply finer discretisation in different areas beneath the excavation. These additional meshes are presented in Fig. 10 and the variation in element sizes is listed in Table 3. For these analyses the nonlocal parameters were changed, adopting a DL = 1.0 m and a corresponding RI = 3.0 m, while all other parameters remained unchanged. As previously discussed, a suitable DL value should be at least equal to or larger than the element size, to include sufficient local devi-atoric plastic strain reference points. Hence, the finer element discretisation employed in this set of analyses, allowed the use of a smaller DL which corresponds to a faster softening rate and hence shorter time to failure. This value is consistent with the investigation of a suitable strain softening rate for London Clay slopes previously discussed and presented in Fig. 8.

Fig. 11 plots contours of accumulated local plastic deviatoric strains for the last converged increment for the 6 nonlocal analyses corresponding to the meshes shown in Fig. 10. Overall the contour plots show the formation of two dominant slip surfaces. Clearly their location is more consistent amongst the analyses with the finest vertical discretisation (i.e. Fig. 11b, d, f), while the coarser meshes show a variation in the development of the slip surfaces.

The predicted horizontal displacement variation with time for the mid-slope surface point is shown in Fig. 12 both for the local and the nonlocal methods. The nonlocal results show that for a DL of 1 m, the analyses employing the coarsest (Fig. 10a or Fig. 4g) and finest (Fig. 10b or Fig. 4h) produce the fastest and slowest times to failure, 38 and 85 years respectively. The remaining four meshes, Fig. 10(c)-(f) produce failure times between these times. The change in horizontal displacement over time varies for each analysis resulting in a displacement just prior to failure of about 0.5 m for the three shortest analyses and of 0.6-0.65 m for the three analyses with the longest times to failure. Clearly, these anal-

Fig. 12. Horizontal displacement of the mid-slope point over time for analyses employing (a) the local and (b) the nonlocal G&S (DL= 1 m and RI = 3m) strain softening models with six meshes of varying discretisation.

Fig. 13. Factor of Safety values evaluated immediately after excavation computed with the local and nonlocal strain-softening methods.

Fig. 14. The effect of varying the radius of influence, RI within a strain softening cutting analysis of DL = 1.0 m.

yses demonstrate some mesh dependency for the nonlocal method, however the nonlocal strain softening method still provides a significant improvement over the local strain softening model. This is demonstrated by the results from 6 local method analyses employing the meshes in Fig. 10. There is a significant variation for local analyses, in time to failure from 15 years up to 270 years and horizontal displacement from 0.35 m to 0.85 m. The relative development of horizontal displacement just prior to failure is also inconsistent for all the local analyses compared to the nonlocal results. Clearly, the local strain softening results are very mesh dependent.

of Safety (FoS) analyses were performed immediately after excavation for the 6 meshes shown in Fig. 10.

Potts and Zdravkovic [20] showed that the most consistent approach in computing the Factor of Safety in FE analysis is to start the analysis with the characteristic strength and at relevant stages of the analysis to gradually increase the safety factor (i.e. the material factor on shear strength) until failure in the soil is fully mobilised. In each increment of the analysis a larger factor of safety (hence lower strength properties) is adopted until failure is reached. The factor of safety is applied by reducing the strength properties of the soil by a factor, F, as shown in Eqs. (6) and (7):

4. Mesh dependency of factor of safety predictions

The analyses presented so far indicate that the prediction of the time to failure is very sensitive to the adopted mesh discretisation when the local strain softening model is employed, while even the time to failure predictions of the non-local G&S model show some mesh sensitivity. In order to further investigate the mesh dependency characteristics of both local and non-local methods, Factor

U = tan

1 itan Uin

where Fs is the current Factor of Safety and Cin and uin are the values of cohesion and angle of shearing resistance respectively in the beginning of the FoS analysis.

In the present study the strength reduction started once the excavation was completed. The FoS was incrementally increased from an initial value of 1.0 until the factored strength resulted in an unstable slope. The time and associated pore pressure changes are stopped on the increment prior to the first factor of safety increment (i.e. at the end of excavation) and the remaining FoS analysis is drained. Excavation is simulated over a time period of 0.25 years and although consolidation is permitted during this period, the time is too short to allow any significant drainage in the low permeability London Clay. In each increment of the analysis, it is the current strength that is factored and so when a strength has been reduced using a strain softening model, the reduced strength is factored. The strain softening behaviour is therefore still captured by the FoS analyses.

Fig. 13 compares the FoS values after excavation as predicted by the local and the nonlocal approaches. The local model results clearly show a significant variation of the FoS values, ranging from 1.370 to 1.507. On the other hand, the nonlocal predictions show negligible mesh dependence with FoS values ranging only from 1.408 to 1.436. The similarity of the nonlocal FoS results for all examined meshes indicates that it is the consolidation process that is responsible for the mesh dependence seen in the predictions of the time to failure. In the consolidation analyses the formation of the slip surface was found to be very sensitive to the pore pressure distribution leading to some mesh dependence even when nonlocal regularisation was employed. However, the nonlocal regularisation eliminated the mesh dependence in the case of the drained FoS analyses where the softening is triggered through the strength reduction.

Accumulated nonlocal plastic deviatoric strain

0 0.05 0.1 0.15 0.2 0.25 0.3 0.35 0.4

Fig. 15. Strain distribution for a vertical cross sections passing through the mid-slope point for the last increment of the cutting analyses varying the RI parameter. DL = 1 m.

5. Radius of influence and relative computational cost

The radius of influence, RI, is an optional parameter that limits the distance from the calculation point of the neighbouring strains that will be included in the nonlocal strain calculation. A small RI reduces the time of the analysis by reducing the number of elements that will be included in the nonlocal calculation. Summersgill [18] investigated the sensitivity of nonlocal biaxial compression analyses on RI, showing that for an RI three times DL or greater, the reaction load response and strain distribution is almost indistinguishable from an analysis performed with no RI specified, i.e. where all neighbouring strains contribute to the calculation of a nonlocal strain no matter how far away they are from the position for which the nonlocal strain is being evaluated.

In the present study the impact of RI both in the accuracy of the solution and the computation cost is examined for the stiff clay cutting analyses. The same soil properties for London Clay as specified in Table 2 are adopted herein for the mesh shown in Fig. 10c (thin below base) with a DL = 1.0 m and various values of RI. The radius of influence was specified as 3.0 m, 4.5 m and 6.0 m and the results are compared against an analysis with no RI specified. The horizontal displacement of the midslope point with time for these four analyses is presented in Fig. 14. The analyses results agree until approximately 20 years after excavation, but show a modest variation from that point onwards. The RI = 3.0 m and no RI analyses give failure at just over 40 years, for RI = 4.5 m and 6.0 m the failure is just over 35 years, within one year of each other. The relative development of slip surfaces within the slope can be used to explain the difference in horizontal displacement with time. The accumulated local plastic deviatoric strain distribution for the last converged increment for these four analyses is shown for a vertical cross section through the mid-slope point in Fig. 15. The corresponding strain contours for the entire model of these analyses are presented in Fig. 16. Figs. 15 and 16 serve to demonstrate the impact of RI on the failure mechanism, indicating

a wider variation of results than was expected based on the biaxial compression results of Summersgill [18].

The sensitivity of the excavated slope results to the RI parameter suggests that the development of strains and therefore slip surfaces within a slope still show a variability when employing the nonlocal strain softening model. It should be noted that these analyses employ the coupled consolidation formulation and as previously discussed, the non-local method shows a mesh dependence in that case. In any case though, the critical slip surface is the same for all four analyses and it is only a secondary deeper slip surface that has developed to a different extent for each analysis and this has caused the difference in the results. The secondary slip surfaces vary in terms of position and development of local and nonlocal strains before slope failure. The analyses with the greatest concurrent development of slip surfaces, RI = 3 m and no RI, take the longest time to reach failure. The strains developed are shown for no RI in Fig. 16(d) and RI = 3 m in Fig. 16(a) and both in Fig. 15. The development of multiple slip surfaces seems to prolong the softening in the slope analyses. The position and strains on the critical slip surface do agree, (see Fig. 15), providing confidence in the repeatability of the results.

The decision of which value of RI to employ for the slope analyses should not be made solely on the similarity of the results to a no RI analysis, a consideration of the time efficiency RI provides is also of interest. Fig. 17 compares the Central Processing Unit (CPU) time in terms of excavation CPU time, total CPU time and total CPU time per increment. For each of these three types of time data, the comparison is made as both a percentage of the analysis with no RI (red2, smaller values) and as percentage of the RI = 6 m analysis time (blue, larger values). The time result for an analysis with the same soil and model parameters, but employing a local strain softening

2 For interpretation of color in Fig. 17, the reader is referred to the web version of this article.

Fig. 16. Local plastic deviatoric strain contours for excavated slopes analyses with various values for the RI parameter. DL = 1 m, excavation 10 m deep for the last converged increment.

model provides further context to assess the time cost of a nonlocal strain softening model and the use of the RI parameter to reduce this time. The analysis employing the local strain softening model is faster than the nonlocal strain softening analyses, although it leads to only a 50% time saving compared to the RI = 3 m analyses. The RI = 3 m analyses are performed in just over 50% of the time for RI = 6 analyses and 70-80% of the RI = 4.5 m analyses. This is a significant time improvement if a large set of slope analyses are to be performed. However, even a large RI of 6 m, 6 times the DL value, provides a significant time saving of around 90% compared to the analysis with no RI. If the final position of the secondary slip surfaces in the slope is important, a RI larger than 6 multiples of DL may provide the required result and would still provide a significant time improvement over not specifying an RI.

The value of RI for this type of analysis should be selected in consideration of the importance of the development of secondary slip surfaces, as well as the comparative time saving provided by the RI value. The use of a radius of influence, RI provides a considerable time advantage, especially for boundary value problems in finite element analyses. An RI of three times the DL value should provide reasonable accuracy of results, as the given form of the

weighting function implies that the contribution of strains at three DL from the point of calculation is very small.

4. Conclusions

In the present study, the use of nonlocal regularisation is examined in a coupled consolidation problem of an excavated slope in a strain softening material. The considered boundary value problem allows for a thorough evaluation of the nonlocal regularisation approach, as it does not entail any kinematic restraint on the slip surface development and therefore the location of the slip surface is not predetermined. Following the recommendations of Summersgill [18], who showed the superiority of the G&S nonlocal strain softening model over two other nonlocal approaches, the G&S method was employed exclusively herein. This model requires the specification of one additional parameter, the defined length DL, which influences the rate of softening. In addition, the optional radius of influence parameter, RI, can be specified to reduce the number of local strains referenced in the nonlocal strain calculation and thus increase the efficiency of the analyses.

Fig. 17. Comparison of time for slope analyses to be performed, comparing the use of local strain softening model to nonlocal G&S strain softening model and selected RI value. (a) CPU time to perform excavation of slope during 25 increments (b) total CPU time for analyses, the number of increments varies for each analysis (c) total CPU time per increment.

The mesh dependency of both the G&S and the conventional local strain softening model was investigated initially for biaxial compression analyses modelling London Clay and then for a range of mesh layouts and element sizes in cutting slope analyses modelling the same material. The biaxial compression results showed that the softening rate in the local strain softening model is element size dependant and that in problems with uniform element discretisation one can achieve the desired softening rate by adjusting accordingly the element size and strain softening limits. In the present study the "target" softening rate for London Clay is based on the findings of Potts et al. [15] resulting from a calibration of their local strain softening analyses on simple shear tests. The nonlocal biaxial compression results showed negligible mesh dependence, confirmed the influence of DL on the softening rate and therefore highlighted the importance of selecting an appropriate value of DL. The DL= 1.05 m analyses produced a close match to the "target" softening rate and therefore this value is considered suitable to represent the softening rate of London Clay. It should be noted that this value is specific to the material softening properties specified herein and with consideration of the application to the subsequent cutting analyses.

In the cutting analyses, the local strain softening results are very dependent on the adopted mesh layout (i.e. element arrangement) and the number of elements. The nonlocal strain softening constitutive model makes a significant improvement on the mesh dependency for a range of mesh layouts and element sizes. The nonlocal analyses are not entirely mesh independent, but the predicted time to failure and horizontal displacement over time are much more consistent compared to analyses that employ the local strain softening constitutive model.

A further parametric study which focused on the computation of the Factor of Safety after excavation for various mesh arrangements, showed that for drained conditions the nonlocal regularisation actually eliminates the mesh dependence shown by the conventional local strain softening model. The similarity of the nonlocal FoS results for all examined meshes indicates that it is the consolidation process that is responsible for the mesh dependence observed in the previous predictions of the time to failure.

Finally a range of RI values was examined in cutting analyses in order to investigate the impact of this parameter on the accuracy of the solution and the computational cost. An RI of 3 multiples of DL was found to provide a suitable compromise between accuracy and time saving. When applied to a slope analysis with DL = 1 m, the RI = 3 m analysis took 5-6% of the time for the no specified RI analysis, 50-60% of the RI = 6 m analysis and approximately twice as long as the conventional local analysis. The RI parameter was found to have an effect on the secondary slip surfaces developed in the slope, although overall the results and the critical slip surface were similar. Overall, the present study shows the superiority

of the nonlocal approach over the conventional local method and the necessity of its use in strain softening problems.

Acknowledgements

This work was part of the PhD research of the first author at Imperial College London through an Industrial Cooperative Award in Science & Technology (CASE) jointly funded by EPSRC, UK and Geotechnical Consulting Group LLP, UK. This support is gratefully acknowledged.

References

[1] Galavi V, Schweiger HF. Nonlocal multilaminate model for strain softening analysis. Int J Geomech 2010;10:30-44.

[2] Vardoulakis IG, Sulem J. Bifurcation analysis in geomaterials. London: Chapman & Hall; 1995.

[3] Bazant ZP, Jirasek M. Nonlocal integral formulations of plasticity and damage: Survey of progress. J Eng Mech 2002;128(11):1119-49.

[4] Zervos A, Papanastasiou P, Vardoulakis I. Finite element displacement formulation for gradient elastoplasticity. Int J Numer Meth Eng 2001;50 (6):1369-88. http://dx.doi.org/10.1002/1097-0207(20010228)50:6.

[5] Di Prisco C, Imposimato S, Aifantis EC. A visco-plastic constitutive model for granular soils modified according to non-local and gradient approaches. Int J Numer Anal Method Geomech 2002;26(2):121-38.

[6] Papanicolopulos SA, Zervos A. Numerical solution of crack problems in gradient elasticity. Eng Comput Mech 2010;163(EM2):73-82. http://dx.doi. org/10.1680/eacm.2010.163.2.73.

[7] Eringen AC. On nonlocal plasticity. Int J Eng Sci 1981;19:1461-74.

[8] Vermeer PA, Brinkgreve RBJ. A new effective non-local strain measure for softening plasticity. In: Chambon R, Desrues J, Vardoulakis I, editors. Localization and bifurcation theory for soils and rocks. Rotterdam, The Netherlands: Balkema; 1994. p. 89-100.

[9] Lu X, Bardet J, Huang M. Spectral analysis of non-local regularization in two-dimensional finite element models. Int J Numer Anal Meth Geomech 2012;36:219-35.

[10] Jirâsek M, Rolshoven S. Comparison of integral-type nonlocal plasticity models for strain-softening materials. Int J Eng Sci 2003;41(13-14):1553-602.

[11] Di Prisco C, Imposimato S. Nonlocal numerical analyses of strain localisation in dense sand. Math Comput Model 2003;37(5-6):497-506.

[12] Potts DM, Zdravkovic L. Finite element analysis in geotechnical engineering: theory. London: Thomas Telford; 1999.

[13] Potts DM, Dounias GT, Vaughan PR. Finite element analysis of progressive failure of Carsington embankment. Geotechnique 1990;40(1):79-101.

[14] Kovacevic N. Numerical analysis of rockfill dams cuts slopes and road embankments Ph.D. thesis. London: Imperial College; 1994.

[15] Potts DM, Kovacevic N, Vaughan PR. Delayed collapse of cut slopes in stiff clay. Géotechnique 1997;47(5):953-82.

[16] Vaughan PR, Walbancke HJ. Pore pressure changes and the delayed failure of cutting slopes in over-consolidated clay. Géotechnique 1973;23(4):531-9.

[17] Vaughan PR. Assumption, prediction and reality in geotechnical engineering. Géotechnique 1994;44(4):573-609.

[18] Summersgill FC. Numerical modelling of stiff clay cut slopes with nonlocal strain regularisation PhD Thesis. London: Imperial College; 2015.

[19] Bazant ZP, Belytschko TB, Chang TP. Continuum model for strain softening. J Eng Mech 1984;110(12):1666-92.

[20] Potts DM, Zdravkovic L. Accounting for partial material factors in numerical analysis. Géotechnique 2012;62(12):1053-65.