Chapter 6

Numerical simulation

6.1 Introduction

To gain a better understanding of the shaft’s performance, the monitoring results (Section 5.3) are compared to a parametric study conducted with the finite difference package FLAC. This chapter provides the details and input parameters for the simulations.

6.2 Overview

FLAC stands for ‘Fast Lagrangian Analysis of Continua’ and is a two-dimensional explicit finite difference program.

The shaft was modelled using axisymmetric mode where the z-axis is the axis of symmetry. The model as displayed in Figure 6.1 is approximately 93 m wide and 120 m deep and contains 12960 elements (Table 6.1). All soil elements are 1 m x 1 m; only the first column of elements is only 0.85 m wide. As FLAC does not allow the use of structural elements in axisymmetric mode, the wall was modelled using solid elements with six 0.2 m wide elements across the width of the wall (Figure 6.1 close-up view). An interface was placed between soil and wall extending 1 m below the toe of the wall. The applied mechanical and fluid boundary conditions are displayed in the figure, as well as the assigned pore water pressure distribution. Further details can be found in Figure 6.1 and Table 6.1.

FLAC was used to run a model similar to the PLAXIS design, which is referred to as ‘FLAC1’. Further analyses incorporating advanced features and the parametric study were run in with the model ‘FLAC2’.

Figure 6.1: Overview of FLAC model

Table 6.1: FLAC model details

6.3 Modelling procedure

Ideally a fully coupled modelling approach would be used. However due to the large permeability contrast the associated time step was very small and it would be very time-consuming to run the full excavation problem. Hence, an uncoupled approach where the fluid flow calculations are uncoupled from the mechanical calculations was applied. This does however not allow real time to be involved in the calculations.

The two model approaches suggested in Itasca (2011b, Section 1.8.2.1) were combined for the FLAC runs carried out in this study: The London Clay and the Lambeth Group Clay were modelled as fully undrained and all other layers as drained. Effective stress based soil models were used and the undrained condition was modelled by adding the stiffness of the pore fluid to the overall material stiffness. Only mechanical calculations were performed, as the fluid flow is negligible and can be switched off. For the drained material the analysis is conducted in two steps: (1) the fluid step where the steady-state pore pressure field is determined and (2) the mechanical step where the model is cycled to equilibrium under the constant pore pressure field.

To model the full excavation problem, the following steps were carried out:

Initial equilibrium: In the first step the mesh is generated, as displayed in Figure 6.1. The mechanical and fluid boundary conditions, soil properties and initial stresses are assigned. Gravity, water density and fluid bulk modulus are also specified (Table 6.1). The model is then cycled to equilibrium.

Wall installation: For the FLAC1 model the wall is assumed to be wished-inplace (WIP). In this case, the procedure detailed below is skipped. For some FLAC2 models the wall installation is modelled. A similar approach as suggested by Lings et al. (1994), Gourvenec and Powrie (1999) and Conti et al. (2012) is applied here. The solid elements of the future wall are removed and the bentonite pressure is applied as total stress, as shown in Figure 6.2 by the orange arrows. After this has reached equilibrium additional pressure which represents the wet concrete is applied, as shown in Figure 6.2 by the purple arrows.

Wall solid elements & interface: In this step the solid elements of the wall are assigned their material model, properties and initial stresses and the interface between soil and wall is activated. The interface behaviour is described with an elastic-plastic Mohr-Coulomb criterion. Its properties are based on the PLAXIS model, which used an interface strength reduction factor Rinter of 0.8. The interface cohesion value cinter is obtained by Rinter · c ′ (Brinkgreve et al., 2009) where c ′ is the ground cohesion; this yields 16 kPa for the interface between the wall and the Chalk and 0 for all other layers. An average friction angle of φ ′ = 35° was used to determine the interface friction angle φinter = arctan(Rinter · tanφ ′ ) (Brinkgreve et al., 2009) to approximately 30° which was assigned to all layers.

Excavation: The excavation process is carried out in 1 m steps.

For the excavation above the London Clay the pore pressure at excavation level and 1 m below is set to 0 and fixed, i.e. drainage is allowed (Figure 6.3). FLAC automatically adjusts the total stresses; a flow step is not required as the clay layers act as a impermeable barrier. The top row of elements is then removed and the model is cycled to equilibrium with the mechanical calculations set on and the flow calculations set off.

When excavation reaches the clay layers dewatering is not required anymore. The top row of elements is removed and the fluid boundary condition at the excavation level is set not to allow drainage. The cycling to equilibrium is again conducted with mechanical calculation on and flow calculations off.

For the excavation below the clay two sub-steps are required: (1) the flow step and (2) the mechanical step. For the former the pore pressure is set to 0 at excavation level and 1 m below, similar to Figure 6.3. The model is then stepped to steady-state with mechanical calculations off and flow calculations on. For (2) the top row of elements is removed, the mechanical calculations are turned on, the flow calculations are set off and the model is cycled to equilibrium.

Figure 6.2: Modelling diaphragm wall installation

Figure 6.3: Modelling of excavation

6.4 Constitutive soil models

6.4.1 Overview

The Abbey Mills shaft penetrates several ground layers (Figure 4.1b). The thinner layers, i.e. the Alluvium, the Harwich Formation and the Calcrete are not explicitly modelled. For the simplified FLAC1 model the remaining seven groups, i.e. Made Ground, River Terrace Deposits, London Clay, Lambeth Group Clay, Lambeth Group Sand, Thanet Sand and Chalk, were assigned the Mohr-Coulomb (MC) model. The input parameters can be found in Table 6.4.

In FLAC2 the Made Ground, River Terrace Deposits, London Clay, Lambeth Group Clay are also modelled with the MC model. For the clays a hyperbolic stiffness degradation function was implemented. Lambeth Group Sand and Thanet Sand are modelled with a Strain-Hardening/ Softening (SHS) model; the Hoek-Brown (HB) model was used for the Chalk. This is further detailed below and an overview as well as the input parameters for the soil layers are provided in Table 6.5.

6.4.2 Mohr-Coulomb model

The MC model simulates elastic soil behaviour before yield. The failure envelope that defines the onset of yield comprises two functions: (1) the shear yield function $f^s$ (Equation 6.1) and (2) the tension yield function (Equation 6.2),

where $\sigma'_1$ and $\sigma'_3$ are the principal effective stresses (compressive stresses are negative), $\sigma^t$ is the tension cut-off, c’ is the cohesion intercept, $\phi'$ is the effective stress friction angle and

Yield occurs when either $f^s$ or $f^t$ are equal to 0.

The behaviour after yield is governed by a non-associated flow rule (for shear flow) and an associated flow rule (for tension flow). Further details can be found in Itasca (2008).

All FLAC1 layers are modelled with the MC model and input parameters are tabulated in Table 6.4. In FLAC2 the MC model is assigned to the Made Ground and River Terrace deposits and the input parameters are the same as in FLAC1 (Table 6.5).

6.4.3 Stiffness degradation function

FLAC’s integrated MC model detailed above was assigned to the London Clay and Lambeth Group Clay and small-strain stiffness degradation before yield was incorporated with a so-called ‘FISH-function’. FISH is a programming language in FLAC allowing the user to write functions, new constitutive models, etc.

The FISH-function calculates the shear stress τ according to Equation 6.4 (Pyke, 1979),

where Gmax is the shear modulus at very small strain, a is a constant, n depends on loading/ unloading conditions, γ is the current shear strain, τref and γref are the reference shear stress and strain which are set to the values of the last strain reversal as shown in Figure 6.4a. The equation can be simplified for first time loading (γref = τref = 0) to

The stress-strain relationship used here is displayed in Figure 6.4a; a perfect loop closure was adopted, as no cyclic loading is expected. The value of n is either 1 (for first time loading) or 2 (for unloading and reloading). After the maximum strain is exceeded, n is set to 1 again as shown in Figure 6.4a.

The secant shear modulus G^s is then calculated as follows,

To avoid numerical instabilities G^s ≥ Gmax/10.

Gmax and a were determined by matching G^s with triaxial test data. For the London Clay the undrained triaxial test 25gUC carried out by Gasparre (2005) was picked; details on the test are given in Table 6.2. This test was picked, as it was conducted with a mean effective stress in a similar range of what is expected in the London Clay at Abbey Mills. Gasparre (2005) plots shear modulus versus shear strain εs; ε^s had to be converted to match FLAC’s definition of shear strain γ as follows

Figure 6.4b shows Gasparre’s (2005) measured shear modulus G^s versus γ. The hyperbolic function with Gmax = 80 MPa and a = 4500 matched the triaxial data best.

Data for the Lambeth Group Clay was taken from CIRIA C583 (Hight et al., 2004) which provides values of undrained Young’s modulus normalised by the mean effective stress. This was converted into the shear modulus with a mean effective stress of 450 kPa (mid-depth of the 10 m thick layer). The results are plotted in Figure 6.4c. The best fit was achieved with Gmax = 500 MPa and a = 18000.

The complete input parameters for the London Clay and Lambeth Group Clay are listed in Table 6.5.

Table 6.2: Details of triaxial test used to determine stiffness degradation curve for London Clay (Gasparre, 2005)

(a) Stress-strain relation

(b) Triaxial data by Gasparre (2005) and hyperbolic match

(c) Data from CIRIA C583 (Hight et al., 2004) and hyperbolic match

Figure 6.4: Stiffness degradation for London Clay and Lambeth Group Clay

6.4.4 Strain hardening/ softening model

The SHS model is based on the MC model (Section 6.4.2) with elastic material behaviour before yield. After the onset of yield the cohesion, friction, dilation and tensile strength can harden or soften. Piecewise linear functions of these parameters versus axial plastic strain εap can be input into FLAC (Itasca, 2008).

The elastic parameter for Lambeth Group Sand and Thanet Sands were taken from the design report (i.e. same as in Table 5.1). The plastic behaviour is modelled according to one of Ventouras’ (2005) triaxial tests. The test RS17 (details see Table 6.3) is used to curve-fit the piecewise linear functions. This test was picked as the sample’s confining pressure was close to what would be expected and the location of the site was close to Abbey Mills. The stress ratio q/p ′ and the volumetric strain ε^v versus the axial strain ε^a were extracted from Ventouras’ (2005) triaxial test (Figures 6.5a and 6.5c). The onset of yield was picked at an axial strain ε^a = 5 % and a volumetric strain ε^v = 2 %. This yields a Poisson’s ratio ν of 0.3, similar to the design parameters in Table 5.1, which shows that the selected onset of yield was appropriate.

To obtain the input functions for FLAC, the effective stress friction angle φ ′ and the dilation angle ψ were extracted according to the following equations and then plotted versus axial plastic strain, εap = ε^a −5% (Figures 6.5b and 6.5d).

Equation 6.9 was derived from Vermeer and de Borst’s (1984) bilinear idealisation of triaxial test results.

The piecewise linear functions were fitted to the data as shown in Figures 6.5b and 6.5d. These were used for the Thanet Sand and the Lambeth Group Sand; the other input parameters apart from these are given in Table 6.5.

Table 6.3: Details of triaxial test used for SHS model for Thanet Sand and Lambeth Group Sand (Ventouras, 2005)

Figure 6.5: SHS model: fitting piecewise linear functions to the effective stress friction angle and dilation angle from test RS17 taken from Ventouras (2005)

6.4.5 Hoek-Brown model for Chalk

The HB model is a failure criterion for rock material; it is not a full constitutive model. In addition to the failure criterion, which is in FLAC taken as a yield surface, the model requires a stress-strain relation before yield and a flow rule after yield.

Yield onsets when the yield function $f^{hb} = 0$ ,

where $\sigma_1'$ and $\sigma_3'$ are the principal effective stresses, $\sigma_{ci}$ is the uniaxial compressive strength and $m_b$ , s and $a_{hb}$ are material constants (Itasca, 2008). These constants can be related to the Geological strength index GSI, the disturbance factor $D_d$ and the material constant $m_i$ (Hoek et al., 2002, 1998). With GSI = 90 (depending on friction angle, here $\phi' = 35^\circ$ , Table 5.1), $D_d = 0.5$ (average disturbance) and $m_i = 7$ (for Chalk), the material parameters were determined to $m_b = 4.3$ , s = 0.26 and $a_{hb} = 0.5002$ . The value for $\sigma_{ci}$ was obtained as follows,

where $\sigma_c$ is the unconfined compressive strength, which was measured to be approximately 4 MPa in the Geotechnical Interpretative Report (GIR) (Morgan Sindall, 2012).

The material behaves elastic before yield. After yield the user can decide whether the material coefficients $\sigma_{ci}$ , $m_b$ , s and $a_{hb}$ remain constant or soften with increasing plastic confining strain. For this study the coefficients remain constant after yield.

6.5 Constitutive wall model

In the PLAXIS design the wall was assumed isotropic where the stiffness in circumferential direction is equal to the vertical stiffness. This may not be adequate as discussed in Section 2.4.2.1.

For the parametric study in the FLAC2 model it is necessary to be able to reduce the circumferential stiffness. FLAC does not provide a constitutive model that allows a lower circumferential stiffness; it however provides the option to write a user-defined constitutive model in its own programming language, called FISH (Itasca, 2011a).

The compliance matrix $[D]^{-1}$ for a material with a stiffness E and a reduced circumferential stiffness $E_{\theta}$ was derived from Hooke’s Law and is shown below,

where the strain vector $\{\dot{\varepsilon}\}$ comprises the incremental radial strain $\Delta\varepsilon_{\rm r}$ , vertical strain $\Delta\varepsilon_{\rm z}$ , circumferential strain $\Delta\varepsilon_{\theta}$ and shear strain $\Delta\varepsilon_{\rm rz}$ and the stress vector $\{\dot{\sigma}\}$ comprises the incremental radial stress $\Delta\sigma_{\rm r}$ , vertical stress $\Delta\sigma_{\rm z}$ , circumferential stress $\Delta\sigma_{\theta}$ and shear stress $\Delta\sigma_{\rm rz}$ . $\nu_{\rm rz}$ is the Poisson’s ratio of the material in the plane of isotropy and $\nu_{\theta \rm r}$ is the value between the anisotropic and isotropic direction.

6.6 FLAC runs and parametric study

6.6.1 FLAC1

As mentioned previously FLAC1 is a simplified model, similar to the PLAXIS design. The input parameters for the analysis are given in Table 6.4. FLAC requires the dry density $\gamma_{\rm dry}$ which was obtained through the bulk density $\gamma_{\rm bulk}$ (Table 5.1) and the porosity from the GIR (Morgan Sindall, 2012). Permeability values could also be estimated from the GIR. All other parameters were the same as in the PLAXIS design with two exceptions: River Terrace Deposits layer was given a dilation angle $\psi = \phi'/2$ and $K_0 = 0.6$ . The stratigraphy was slightly altered to match the borehole record.

The wall is modelled as WIP with reduced $K_0$ values (Table 6.4) to take account for stress relief during wall installation.

The excavation is modelled in steps of 1 m (instead of 10 m as PLAXIS) so that interpolation to the fibre optic (FO) excavation levels is not required.

6.6.2 FLAC2

As mentioned above the FLAC2 model was used to carry out a parametric study. Several analyses with varied parameters were conducted. An overview is provided in Table 6.6; further details on the parameters considered is given below:

• Soil models: The PLAXIS design assumed MC model for all ground layers. The effect of using advanced soil models, as detailed in Section 6.4, is investigated.
• Wall installation/ WIP: It is examined whether modelling the wall installation or assuming a WIP wall has an effect on the results.
• Wall stiffness: The PLAXIS design used a short-term Young’s modulus of 37 GPa for the C50/60 concrete. Due to creep and relaxation over time a Young’s modulus of 26 GPa (70 % of the original value as adivsed by CIRIA C580, Gaba et al., 2003) was used for a long-term analysis. The FLAC2 model is used to investigate how a reduced Young’s modulus influences the results.
• Wall anisotropy: As discussed in Section 2.4.2.1, there is a high uncertainty of the wall stiffness in circumferential direction due to the discontinuous nature of diaphragm walls. Several FLAC2 models were run with decreased circumferential stiffness.
• Wall thickness: The PLAXIS design assumed a 1.2 m thick ring structure as shaft wall. While the walls are likely to be 1.2 m thick, the verticality tolerance may have reduced the overlap between the panels and hence the thickness of the ring structure. To investigate the impact of a reduced thickness of the ring structure, two models were run: (1) FLAC2-tred with a reduced thickness of 1.0 m (constant over depth) and (2) FLAC2-tvert with full thickness of 1.2 m at the top and 0.9 m at the bottom (Figure 6.6). This is a conservative simplification, as each panel is in reality still 1.2 m thick and only the overlap is reduced by the verticality tolerances. It was however not possible to model this in a two-dimensional program and hence the real wall thickness, which would be beneficial for the wall’s bending stiffness, was neglected.

Chalk stiffness: The Chalk stiffness is highly uncertain. The GIR (Morgan Sindall, 2012) provides the results of high pressure dilatometer tests which indicate Young’s moduli between approximately 500 MPa and 4000 MPa depending on the strain level. The effect of varying the Chalk stiffness between the upper and lower limit is investigated.

Lateral earth pressure coefficient K⁰ of the Chalk: The lateral earth pressure in Chalk is also highly uncertain. CIRIA C574 (Lord et al., 2002) states that few measurements of horizontal stresses have been made and that design practice usually assumes a K⁰ between 0.5 and 1.0. However, Watson et al. (1999) reports values between 0.5 and 1.5 based on self-boring pressuremeter tests undertaken for the North Downs Tunnel; Rankin and May (1996) back-analyses K⁰ values from the Channel Tunnel up to 1.7. Hence FLAC is run with varying K⁰ values in the Chalk to investigate the influence of this parameter.

Figure 6.6: Reducing wall thickness due to verticality of panels

6.7 Assumptions and limitations

Uncoupled modelling approach: A fully coupled model involving real time would not have been feasible due to the long computational time in FLAC. Hence an uncoupled approach was used which assumed that, during the excavation stage, the clayey soil layers are ‘undrained’ and the sandy soil layers are ‘fully drained’. In reality, depending on rate of excavation and permeability of the ground, some soil layers are ‘partially’ drained.

Simplified two-dimensional model: FLAC2D in axisymmetric mode was used. This imposes several limitations:

• It is not possible to account for any non-uniform ground conditions or construction activities (e.g. culvert construction, dewatering trial) and hence investigate three-dimensional deformation.
• In reality the wall installation is carried out panel by panel, while these simulations model the wall installation of all panels at the same time.
• Concrete wall linear elastic: The concrete wall was modelled as linear elastic. In reality concrete is non-linear and can crack, which was not included in the model.
• Soil parameters: The soil parameters are based on laboratory experiments and previous experience as detailed in the GIR (Morgan Sindall, 2012). The wide scatter of some experiments, the not always transferability of empirical data and the spatial variability of parameters and ground layers around the shaft cannot be taken into account.

Table 6.4: FLAC1 input parameters

Table 6.5: FLAC2 input parameters

with stiffness degradation parameters change after onset of yield

Table 6.6: Parametric study

^* Linear varying, e.g. $3.7 \triangleright 37(27 \text{ m}) = 3.7 \text{ GPa}$ at top of wall and 37 GPa at 27 mbgl in the wall ** Initially 3.7 GPa, linearly increasing to 37 GPa at 5 mm lateral wall movement