Quantitatively assessing the pre-grouting effect on the stability of tunnels excavated in fault zones with discontinuity layout optimization: A case study

Abstract

Pre-grouting is a popular ground treatment strategy utilized to enhance the strength and stability of strata during the excavation of a tunnel through a fault zone. Two important questions need to be answered during such an excavation. First, how should the grouting size be determined? Second, when should excavation begin after grouting? These two questions are conventionally addressed through empirical experience and standard criteria because a reliable quantitative approach, which would be preferable, has not yet been developed. To address these questions, we apply a recently proposed numerical approach known as discontinuity layout optimization, an efficient node-based upper bound limit analysis method. A case study is provided utilizing a tunnel located in a stratum characterized by complicated geological conditions, including soft soil and a fault zone. The factor of safety is used to quantitatively assess the stability of the tunnel section. The influences of the grouted zone thickness and the time-dependent material properties of the grouted zone on the stability of the tunnel section are evaluated, thereby assisting designers by quantitatively assessing the effects of pre-grouting.

Keywords

pre-grouting / stability analysis / factor of safety / discontinuity layout optimization

Introduction

With regard to tunnels constructed in fault zones, the instability and potential collapse of the tunnel ceiling have significantly negative impacts on the project schedule and cost [¹,²]. Pre-grouting constitutes a very commonly used strategy for reinforcing the soil and fractured rock mass in advance, and it enhances the stability of the tunnel section during excavation [³–⁸].The procedure is very clear: a series of boreholes are initially drilled into planar excavation surfaces, after which grout is injected into the soil/rock medium [⁹–¹¹] (see Fig. 1). Following a specific period of time, during which the medium is reinforced and the stability of the tunnel section is guaranteed, the grouted media will be excavated. In this case, the stability of the tunnel section depends on the size of the grouted zone and the material properties of the grouted stratum, both of which will change over time. When the geological conditions of the original stratum are known, designers and engineers seek the answers to two important questions. First, how should the size of the grouted region be determined? Second, when should excavation begin after grouting? These questions are conventionally answered through empirical experience and design criteria. Therefore, a robust and efficient numerical procedure for quantitatively evaluating the pre-grouting process would be highly salient for addressing these two questions; however, such a technique has yet to be developed.

The failures of tunnel sections or tunnel surfaces are attributable to discontinuities, which can be analyzed by nonlinear analysis and limit analysis methods. The former include discontinuity deformation analysis and discrete element methods [¹²,¹³] and phase field methods [¹⁴–¹⁹] in addition to the numerical manifold method [²⁰–²²], strong discontinuity embedded approach [²³–²⁸], cracking elements method [²⁹–³¹], cracking particles method [³²–³⁴], and meshless method [³⁵–³⁷] as well as some sophisticated techniques based on peridynamics [³⁸–⁴²], to name a few. Unfortunately, despite their accuracy and sophistication, most nonlinear analysis techniques suffer from numerical instabilities and low efficiencies relative to limit analysis methods, which preliminarily assume the failure mechanisms and can directly obtain the failure loads and patterns without utilizing loading/unloading processes during the relevant calculations [⁴³–⁴⁶].

Regarding the application of limit analysis to the stability of tunnel sections, researchers have presented various numerical approaches. For example, Leca and Dormieux [⁴⁷] proposed a failure criterion for a tunnel face in a cohesive-frictional soil to determine the upper and lower bounds of the stability of the tunnel face. Mollon et al. [⁴⁸] employed a kinematic approach based on limit analysis theory to analyze the collapse and blow-out face pressures of a circular tunnel driven by a pressurized shield. Similarly, Zhang et al. [⁴⁹] described a new 3D failure mechanism using the limit analysis-based kinematic approach to determine the limit support pressure of the tunnel face. Anagnostou and Perazzelli [⁵⁰] proposed a computational method for assessing the reinforcement of a tunnel face situated within a cohesive-frictional soil with bolts based on limit equilibrium considerations. Han et al. [⁵¹] proposed a new 3D failure mechanism to analyze the limit support pressures of tunnel faces in multilayered cohesive-frictional soils based on limit analysis methods. Chen et al. [⁵²] revealed the failure mechanism and limit support pressure of a tunnel face in dry sandy soil by using the discrete element method (DEM). In addition, Dang and Meguid [⁵³] proposed an efficient finite-DEM to analyze quasi-static nonlinear soil-structure interaction problems involving large deformations in 3D space. Nevertheless, most of these techniques were developed for a continuous medium, and thus, they are not very suitable for highly fractured rock masses. Moreover, these numerical methods are mostly element-based and therefore possess a mesh bias.

In this paper, we apply the recently proposed discontinuity layout optimization (DLO) numerical approach [⁵⁴–⁵⁷] to quantitatively evaluate the stability of pre-grouted tunnel sections. DLO, which belongs to the family of topology optimization methods [⁵⁸–⁶⁰], introduces over one thousand nodes into the domain and connects every pair of nodes as a potential discontinuity for optimization. The target function is defined based on the minimization of energy, and the final solution is obtained automatically [⁶¹–⁶⁴]; accordingly, this approach exhibits numerical stability and a high efficiency [⁶⁵–⁶⁸]. Compared with the aforementioned methods, DLO does not predefine any blocks, elements or failure mechanisms; furthermore, the final activated discontinuities are automatically calculated as results rather than inputs. Unlike element-based methods, DLO is a node-based approach, meaning that it does not exhibit any mesh bias [⁶⁹]. In addition, by defining a factor of safety [⁷⁰,⁷¹], DLO can be used for stability analysis [⁶⁵], and thus, it represents a suitable tool for this study. Here, a tunnel excavated within a zone characterized by complicated geological conditions, including a fault zone and fractured rock with a rock quality designation of RQD<10, is considered. The size of the grouting region and the effect of the grouting time are taken into account. The factor of safety is used to quantify the stability and the quality of grouting. With the calculated results, our study answers the two previously proposed questions and helps to determine the size of the grouted zone and the time required for excavation. We emphasize here that DLO is a patented method [⁷²,⁷³] and that the case study shown in this work is provided to indicate the capabilities of this method “only for engineering attempts rather than financial purposes” .

The remaining parts of this paper are organized as follows. In Section 2, the numerical DLO method is explained with the governing equations, the target functions, and the optimization constraints. In Section 3, the project background of the tunnel and the zone are described in detail. Numerical studies are conducted in Section 4 to evaluate the influences of the size of the grouted zone, the geological parameters and the grouting material properties on the tunnel stability. Finally, the concluding remarks are presented in Section 5.

Governing equations of DLO

The numerical DLO method was presented in Ref. [⁵⁴]. By introducing the factor of safety, DLO can be further utilized for stability analysis [⁶⁵]. Accordingly, we present the formulation of DLO here. In a 2D domain, n nodes are introduced, and each pair of nodes is connected to represent a potential discontinuity for a total number of discontinuities C(n,2). Every potential discontinuity with a potential displacement jumps along both the normal direction and the shearing direction composing the displacement jump field d, which represents the basic unknown in the governing equations. Meanwhile, the factor of safety l, which is defined as an increasing factor of the gravitational acceleration that triggers failure, is introduced. With regard to the energy balance, l = E/W, where E represents the dissipated energy and W is the work done by gravity, both of which are functions of d. Hence, the general governing equations are given as

(1)\begin{array}{c}\hfill \min \mathrm{\hspace{0.17em}}\lambda ={\displaystyle \frac{E\left(\mathit{d}\right)}{W\left(\mathit{d}\right)}},\hfill \\ \hfill \hfill \\ \hfill \mathrm{s}\mathrm{u}\mathrm{b}\mathrm{j}\mathrm{e}\mathrm{c}\mathrm{t}\mathrm{\hspace{0.17em}}\mathrm{t}\mathrm{o}\mathit{A}\mathrm{\hspace{0.17em}}\mathit{d}=\mathit{u},\hfill \end{array}

where l is the target function, and A is the constraint matrix of the optimization equation. Since Ref. [⁶⁵] introduced the unit work constraint and added W(d) = 1 to the constraints, l = E(d).

As presented in Refs. [⁵⁴,⁵⁷] and represented in matrix form in Refs. [⁶⁵,⁶⁹], regarding the utility of DLO in stability analysis, A includes three types of constraints: compatibility constraints, flow rule constraints, and unit work constraints. Furthermore, to eliminate absolute operation in the flow rule constraint, a plastic multiplier field p is also introduced as an unknown. Finally Eq. (1) can be transformed into

(2)\begin{array}{l}\min \mathrm{i}\mathrm{m}\mathrm{i}\mathrm{z}\mathrm{e}\\ E=\mathit{e}\mathrm{\hspace{0.17em}}\mathit{p},\\ \mathrm{s}\mathrm{u}\mathrm{b}\mathrm{j}\mathrm{e}\mathrm{c}\mathrm{t}\mathrm{\hspace{0.17em}}\mathrm{t}\mathrm{o}\\ \left[\begin{array}{lll}\mathit{B}& & \mathbf{0}\\ & \mathit{N}& \\ \mathit{G}& & \mathbf{0}\end{array}\right]\mathrm{\hspace{0.17em}}\left[\begin{array}{l}\mathit{d}\\ \mathit{p}\end{array}\right]=\left[\begin{array}{l}\mathbf{0}\\ \mathbf{0}\\ \mathbf{1}\end{array}\right],and\{\begin{array}{l}\forall i\in {{\displaystyle \mathit{I}}}_s,{{\displaystyle \mathit{s}}}_i=\mathrm{0}\\ \forall i\in {{\displaystyle \mathit{I}}}_n,{{\displaystyle \mathit{n}}}_i=\mathrm{0}\\ \mathit{p}\ge \mathbf{0}\end{array}\end{array}

where I_s and I_n are sets of the discontinuity index fixed along the shear and normal directions, respectively. In Eq. (2), the dissipated energy is calculated by

(3)E(\mathit{d})=E(\mathit{p})=\mathit{e}\mathrm{\hspace{0.17em}}\mathit{p}=\underset{j}{\mathrm{\Sigma }}[{{\displaystyle \mathit{c}}}_j{{\displaystyle \mathit{l}}}_j\mathrm{\hspace{1em}}{{\displaystyle \mathit{c}}}_j{{\displaystyle \mathit{l}}}_j]\cdot {{\displaystyle \mathit{p}}}_j,

where c_j is the cohesion of the inner discontinuity j, and l_j is the cross section of the discontinuity j with a value equal to the length of j in two dimensions.

B is obtained by assembling the local constraint matrix B_i, where i denotes a potential discontinuity i, as

(4){{\displaystyle \mathit{B}}}_i=\left[\begin{array}{l}\begin{array}{ll}{{\displaystyle \alpha }}_i& -{{\displaystyle \beta }}_i\end{array}\\ \begin{array}{ll}{{\displaystyle \beta }}_i& {{\displaystyle \alpha }}_i\end{array}\\ \begin{array}{ll}-{{\displaystyle \alpha }}_i& {{\displaystyle \beta }}_i\end{array}\\ \begin{array}{ll}-{{\displaystyle \beta }}_i& -{{\displaystyle \alpha }}_i\end{array}\end{array}\right],

where a_i and b_i are the x-axis and y-axis direction cosines, respectively, for a discontinuity i connecting nodes A and B. Figure 2 shows the polar angles of discontinuities 1 through 4, where a_i = cos q_i and b_i = sin q_i. B d = 0 represents the rigid displacements of the sub-domains separated by all potential discontinuities.

N in Eq. (2) is calculated by assembling the local flow rule matrix in consideration of the Mohr-Coulomb flow law as

(5){{\displaystyle \mathit{N}}}_j\left[\begin{array}{l}{{\displaystyle \mathit{d}}}_j\\ {{\displaystyle \mathit{p}}}_j\end{array}\right]=\left[\begin{array}{l}\begin{array}{ll}\begin{array}{ll}\begin{array}{ll}-\mathrm{1}& \mathrm{0}\end{array}& \mathrm{1}\end{array}& -\mathrm{1}\end{array}\\ \begin{array}{ll}\begin{array}{ll}\begin{array}{ll}\mathrm{0}& -\mathrm{1}\end{array}& \tan {{\displaystyle \varphi }}_j\end{array}& \tan {{\displaystyle \varphi }}_j\end{array}\end{array}\right]\left[\begin{array}{l}{{\displaystyle \mathit{d}}}_j\\ {{\displaystyle \mathit{p}}}_j\end{array}\right]=\mathrm{0},

where j is the inner potential discontinuity, and {{\displaystyle \varphi }}_i is the corresponding friction angle (see Refs. [⁵⁴,⁶⁵] for details).

The unit constraint matrix G is constructed by considering the work done by gravity attributable to the potential movements of the sub-domains. G is assembled by a local unit constraint matrix as

(6){{\displaystyle \mathit{G}}}_i\mathrm{\hspace{0.17em}}{{\displaystyle \mathit{d}}}_i={{\displaystyle G}}_i\mathrm{\hspace{0.17em}}\mathrm{sgn}({{\displaystyle \alpha }}_i)\mathrm{\hspace{0.17em}}[-{{\displaystyle \beta }}_i\mathrm{\hspace{1em}}-{{\displaystyle \alpha }}_i]\cdot {{\displaystyle \mathit{d}}}_i,

where G_i is the total weight (G_i≥0) of the strip of material lying vertically above the discontinuity i.

As discussed in Ref. [⁶⁵], the final obtained results are composed of the minimum value of l in addition to the corresponding displacement jump field d and plastic multiplier field p. As a direct criterion, if the obtained l>1, the structure is considered to fail only when the constant of gravitational acceleration is greater than the standard value. In other words, l>1 indicates that the structure is stable and self-supported. For details regarding the numerical platform used in this study, interested readers are referred to Refs. [⁷⁴–⁷⁶].

The case study

Project background

The considered tunnel project is the Zhongjiashan tunnel located in the western part of Jiangxi Province, China, which was excavated by the New Austrian tunneling method (NATM). The project includes two tunnel lines passing through a fault zone named F2, which is considered to have a large width of approximately 15-30 m. F2 dips toward the east with an angle of 84°, and it intersects with the axis of the tunnel at 45°. In addition, the stratum is composed of weak soils and fractured rocks with RQD<10. During the excavation of the tunnel, a collapse occurred near the right line at Y K91+ 371 ~ 389 (as shown in Fig. 3); the volume of the collapse was approximately 1200 m³. The formation lithology of the region is divided into four units (see Table 1).

The grouting material used in this project is a cement paste-silicate grout with a cement paste/silicate ratio of CP/S = 1:1, where the cement paste is composed of a cement paste/water ratio of C/W = 1:1. The cohesion and friction angle of the grouted stratum are obtained by direct shear tests on two groups of in situ specimens: one from the typical grouted stratum and the other from the grouted stratum near F2 (as shown in Fig. 4).

Numerical modeling and discretization

The cross section of the tunnel is shown in Fig. 5. We focus on the crown of the tunnel (section 1) to investigate the stability of the tunnel.

The excavation of the Zhongjiashan tunnel mainly included three conditions: i) the tunnel encountered the F2 fault and other soils; ii) the tunnel was completely located within the F2 fault zone; iii) the tunnel was completely located in soil.

The model of the grouted crown for the DLO analysis is shown in Fig. 6 where the surcharge is set as 3128 kPa. The section of the tunnel ceiling is assumed to be an arch with a height H = 4.4 m; the depth of the tunnel (i.e., the distance from the grouted zone to the upper part of the stratum) is D = 20 m, and the grouted zone is assumed to be perfectly circular with a thickness of T = 8 m. The inner surface of the tunnel and the upper part of the stratum are free boundaries, whereas the lateral and bottom parts of the stratum are fixed along the normal direction. The material parameters are listed in Table 2, in which the cohesion of the interface is assumed to be the same as the smaller cohesion values of the two separate sub-domains.

The grid method [⁵⁴] is used for the DLO pre-processing, i.e., a large number of nodes on regular grids are introduced into the domain, and the discontinuities, which are allowed to intersect one another, are built by connecting each pair of nodes; however, the nodes in different types of soils are prohibited from connecting with each other. The potential discontinuities in both a homogeneous stratum and a heterogeneous stratum are illustrated schematically in Fig. 7. Based on the positions of the potential discontinuities, the cohesion, friction angle and overburden are determined for each discontinuity to construct the optimization equation (Eq. (2)). Zhang et al. [⁶⁵] noted that the grid spacing dp has a substantial influence on the final result. After preliminary experimentation, we take dp = 0.7 min this work, forming approximately 2000 nodes and over one million potential discontinuities. Furthermore, to improve the initial results, refining procedures [⁶⁵] are also used in which new nodes are added around the activated nodes during the last optimization step for re-optimization; this refining procedure is repeated 4 times in every simulation.

Excavation conditions of the tunnel

Considering the 3 possible conditions during the excavation of the tunnel, the soil failure mechanisms are shown in Fig. 8.

It is evident that collapse will occur mainly in the soft layer when the tunnel is excavated within heterogeneous soil. The soil at the junction of the soil layer and the excavated tunnel will be damaged and will tend to fall into the fault zone. Furthermore, the factor of safety of the heterogeneous soil is almost identical to that of the soft homogeneous soil, which explains why the stability of the tunnel mainly relies on the geological properties of the soft soil. When T is equal to 2, 4, 6, 8 and 10 m, the failure mechanisms of the damage zone are almost identical, as shown in Fig. 9. The value of l increases with an increase in the thickness of the grouted zone, which means that the stability of the tunnel increases with the thickness of the grouted zone. When the thickness of the grouted zone is relatively small, the damaged soil tends to collapse outwards, and the damage region becomes wider.

To choose a reasonable grouted zone thickness, the factor of safety is obtained with the thickness varying from 2 to 10 m, as shown in Fig. 10. Evidently, although the geological conditions differ, the relationships between the grouted zone thickness and the safety factor can all be fitted by an exponential function. When the cohesion and friction angle of the soil are identical to those under the project conditions, the most economic grouted zone thickness with the best fit can be determined. If the geological properties are substantially different from the abovementioned soil, the most suitable thickness can be chosen according to these two fitting curves once the influence of the cohesion and friction angle on the safety factor is known. The effects of the geological properties and other influencing factors are discussed in the following section.

Influencing factors on the tunnel stability

The dominant influencing factors on the tunnel stability include the geological properties of the soil above the tunnel and the grout-reinforced strength of the soil. Tunnels excavated in uniform and non-uniform strata with different geological conditions are compared in this section.

Effects of geological properties

In the current simulation, the geological properties are varied by changing the cohesion and the friction angle of the soil. Two dimensionless parameters (\Omega \mathrm{1}={{\displaystyle C}}_{s\mathrm{1}}/{{\displaystyle C}}_{s\mathrm{2}}\mathrm{\hspace{0.17em}}and\mathrm{\hspace{0.17em}}\Omega \mathrm{2}={{\displaystyle \varphi }}_{s\mathrm{1}}/{{\displaystyle \varphi }}_{s\mathrm{2}}) are introduced to evaluate the effect of the degree of non-uniformity on the stability of the tunnel. W1 and W2 vary from 1.4 to 2.5 and then 5.0 by respectively decreasing C_s2 and {{\displaystyle \varphi }}_{s\mathrm{2}}. The material parameters of the soil in the fault zone and the grout are kept constant, as shown in Table 2. Following these parametric studies, the failure mechanisms and associated values of l are shown in Figs. 11 and 12. The results indicate that the geological properties of soil 2, which is a weaker soil, have a greater influence on the tunnel stability.

To analyze the effects of the geological properties on the uniform soil, the cohesion varies from 10 to 30 and 60 kPa, while the friction angle {{\displaystyle \varphi }}_{s\mathrm{2}} is kept constant as 6.3°. Meanwhile, when the friction angle varies from 5° to 15° and 30°, the cohesion is equal to 35 kPa. The grout-reinforced strength (C_g = 100 kPa, {{\displaystyle \varphi }}_{s\mathrm{2}} = 20°) and the unit weight of the soil (g = 19.8 kPa/m) are kept constant.

The results of the simulations with different cohesion values are illustrated in Fig. 13, which demonstrates that l increases with an increase in the cohesion of the soil. With the same friction angle, the failure mechanism is almost identical for different cohesion values, but the damage region will become wider with an increase in the cohesion of the soil.

The influence of the friction angle on the tunnel stability is shown in Fig. 14, which illustrates that l increases with an increase in the friction angle. Moreover, when the friction angle is smaller, the failure mechanism of the soil tends to be vertical with respect to the ground surface. When the friction of the grout is smaller than that of the soil, the failure in the grouting interface will expand outwards; meanwhile, if the friction angle of the soil is greater than that of the grout, the failure will fold inwards, and the damage region will become narrower.

Effects of the reinforcement strength of grouting

To account for the effects of time, we follow the work presented in Ref. [⁷⁷] and propose an exponential relationship between the maturity of the grout-reinforced strength and the time expressed as

(7){\displaystyle \frac{C(t)}{{{\displaystyle C}}_{\infty }}}=\mathrm{1}-\exp (-{{\displaystyle \mathit{A}}}_{\mathrm{h}}\mathrm{\hspace{0.17em}}t),

where C_∞ is the final cohesion of the reinforced soil and A_h (day ^{- 1} ) is the growth factor. This exponential relationship is also supported by the experimental results provided in Ref. [⁷⁸]. Figure 15 shows an exponential growth curve with a growth factor of A_h = 0.0822 day ^{- 1} in which the maturity reaches 90% after 28 d; this curve will be used in the following analysis.

When C_∞ = 105.8 kPa, the failure mechanisms of the tunnel sections exhibiting similar failure patterns are illustrated in Fig. 16. Furthermore, the evolution curves of lwith different grouted soil cohesion values are shown in Fig. 17. Evidently, the stability of the tunnel increases with an increase in the cohesion of the grouted soil. The results of this analysis can be employed to assist engineers in determining the proper length of time after grouting to begin excavation.

Conclusions

In this paper, the recently proposed DLO numerical method was applied to analyze the effects of pre-grouting on the stability of a section of the Zhongjiashan tunnel. According to the actual formation lithology and the geometrical parameters of the Zhongjiashan tunnel, 3 conditions that could have been encountered during the tunnel excavation were simulated. The results reveal that the factor of safety l is a minimum when the tunnel is constructed at the junction of the fault zone. The effects of the grouted zone thickness were also simulated, and a fitting curve was obtained.

In addition, multiple parameters, namely, the cohesion and friction angle of the soil and the reinforcement strength of the grouting material, were taken into account. The numerical method and the fitting equation presented in this paper maintain the advantages of limit analysis and can therefore assist designers and engineers in the rapid elimination of inappropriate designs and the selection of a reasonable thickness range according to the obtained values of l to reduce the wasting of grouting materials and lower project costs. Following analyses of the factor of safety and the corresponding influencing factors on the application of grouting prior to tunnel excavation, it is advised to employ a larger pre-grouting region or a grouting material with a greater strength when the tunnel crosses a fault zone or is excavated through weak surrounding rock/soil.

References

[1] Hao Y, Azzam R. The plastic zones and displacements around underground openings in rock masses containing a fault. Tunnelling and Underground Space Technology, 2005, 20(1): 49–61

[2] Anagnostopoulos C A. Laboratory study of an injected granular soil with polymer grouts. Tunnelling and Underground Space Technology, 2005, 20(6): 525–533

[3] Fransson Å. Characterisation of a fractured rock mass for a grouting field test. Tunnelling and Underground Space Technology, 2001, 16(4): 331–339

[4] Yesilnacar M. Grouting applications in the Sanliurfa tunnels of GAP, turkey. Tunnelling and Underground Space Technology, 2003, 18(4): 321–330

[5] Funehag J. Sealing of narrow fractures in hard rock – A case study in hallandsås, sweden. Tunnelling and Underground Space Technology, 2004, 19: 1–8

[6] Funehag J, Gustafson G. Design of grouting with silica sol in hard rock – New design criteria tested in the field, part ii. Tunnelling and Underground Space Technology, 2008, 23(1): 9–17

[7] Faramarzi L, Rasti A, Abtahi S M. An experimental study of the effect of cement and chemical grouting on the improvement of the mechanical and hydraulic properties of alluvial formations. Construction & Building Materials, 2016, 126: 32–43

[8] Sun Z, Yan X, Liu R, Xu Z, Li S, Zhang Y. Transient analysis of grout penetration with time-dependent viscosity inside 3d fractured rock mass by unified pipe-network method. Water, 2018, 10: 1122

[9] Hernqvist L, Fransson Å, Gustafson G, Emmelin A, Eriksson M, Stille H. Analyses of the grouting results for a section of the APSE tunnel at Äspö Hard Rock Laboratory. International Journal of Rock Mechanics and Mining Sciences, 2009, 46(3): 439–449

[10] Lisa H, Christian B, Åsa F, Gunnar G, Johan F. A hard rock tunnel case study: Characterization of the water-bearing fracture system for tunnel grouting. Tunnelling and Underground Space Technology, 2012, 30: 132–144

[11] Lisa H, Gunnar G, Åsa F, Tommy N. A statistical grouting decision method based on water pressure tests for the tunnel construction stage – A case study. Tunnelling and Underground Space Technology, 2013, 33: 54–62

[12] Huang H, Ye G, Qian C, Schlangen E. Self-healingin cementitious materials: materials, methods and service conditions. Materials & Design, 2016, 92: 499–511

[13] Song Y, Huang D, Zeng B. Gpu-based parallel computation for discontinuous deformation analysis (DDA) method and its application to modelling earthquake-induced landslide. Computers and Geotechnics, 2017, 86: 80–94

[14] Wu J Y, Nguyen V P. A length scale insensitive phase-field damage model for brittle fracture. Journal of the Mechanics and Physics of Solids, 2018, 119: 20–42

[15] Wu J Y. Robust numerical implementation of non-standard phase-field damage models for failure in solids article. Computer Methods in Applied Mechanics and Engineering, 2018, 340: 767–797

[16] Wu J, McAuliffe C, Waisman H, Deodatis G. Stochastic analysis of polymer composites rupture at large deformations modeled by a phase field method. Computer Methods in Applied Mechanics and Engineering, 2016, 312: 596–634 doi:10.1016/j.cma.2016.06.010

[17] Wu J Y. A unified phase-field theory for the mechanics of damage and quasi-brittle failure. Journal of the Mechanics and Physics of Solids, 2017, 103: 72–99 doi:10.1016/j.jmps.2017.03.015

[18] Zhou S, Zhuang X, Rabczuk T. A phase-field modeling approach of fracture propagation in poroelastic media. Engineering Geology, 2018, 240: 189–203

[19] Zhou S, Zhuang X, Zhu H, Rabczuk T. Phase field modelling of crack propagation, branching and coalescence in rocks. Theoretical and Applied Fracture Mechanics, 2018, 96: 174–192 doi:10.1016/j.tafmec.2018.04.011

[20] Wu J, Cai Y. A partition of unity formulation referring to the NMM for multiple intersecting crack analysis. Theoretical and Applied Fracture Mechanics, 2014, 72: 28–36

[21] Zheng H, Liu F, Du X. Complementarity problem arising from static growth of multiple cracks and MLS-based numerical manifold method. Computer Methods in Applied Mechanics and Engineering, 2015, 295: 150–171

[22] Zheng H, Xu D. New strategies for some issues of numerical manifold method in simulation of crack propagation. International Journal for Numerical Methods in Engineering, 2014, 97(13): 986–1010

[23] Saloustros S, Cervera M, Pelà L. Challenges, tools and applications of tracking algorithms in the numerical modelling of cracks in concrete and masonry structures. Archives of Computational Methods in Engineering, 2018 (in press)

[24] Saloustros S, Pelà L, Cervera M, Roca P. Finite element modelling of internal and multiple localized cracks. Computational Mechanics, 2017, 59(2): 299–316

[25] Saloustros S, Cervera M, Pelà L. Tracking multi-directional intersecting cracks in numerical modelling of masonry shear walls under cyclic loading. Meccanica, 2017, 53(7): 1757–1776 doi:10.1007/s11012-017-0712-3

[26] Dias-da-Costa D, Alfaiate J, Sluys L, Areias P, Júlio E. An embedded formulation with conforming finite elements to capture strong discontinuities. International Journal for Numerical Methods in Engineering, 2013, 93(2): 224–244

[27] Dias-da-Costa D, Cervenka V, Graça-e-Costa R. Model uncertainty in discrete and smeared crack prediction in RC beams under flexural loads. Engineering Fracture Mechanics, 2018, 199: 532–543

[28] Nikolić M, Do X N, Ibrahimbegovic A, Nikolić Ž. Crack propagation in dynamics by embedded strong discontinuity approach: Enhanced solid versus discrete lattice model. Computer Methods in Applied Mechanics and Engineering, 2018, 340: 480–499

[29] Zhang Y, Zhuang X. Cracking elements method for dynamic brittle fracture. Theoretical and Applied Fracture Mechanics, 2019, 102: 1–9 doi:10.1016/j.tafmec.2018.09.015

[30] Zhang Y, Lackner R, Zeiml M, Mang H. Strong discontinuity embedded approach with standard SOS formulation: Element formulation, energy-based crack-tracking strategy, and validations. Computer Methods in Applied Mechanics and Engineering, 2015, 287: 335–366

[31] Zhang Y, Zhuang X. Cracking elements: a self-propagating strong discontinuity embedded approach for quasi-brittle fracture. Finite Elements in Analysis and Design, 2018, 144: 84–100

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

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

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

[35] Rabczuk T, Belytschko T. Adaptivity for structured meshfree particle methods in 2D and 3D. International Journal for Numerical Methods in Engineering, 2005, 63(11): 1559–1582

[36] Zhuang X, Augarde C, Mathisen K. Fracture modeling using meshless methods and level sets in 3D: Framework and modeling. International Journal for Numerical Methods in Engineering, 2012, 92(11): 969–998

[37] Zhuang X, Augarde C, Bordas S. Accurate fracture modelling using meshless methods, the visibility criterion and level sets: Formulation and 2D modelling. International Journal for Numerical Methods in Engineering, 2011, 86(2): 249–268

[38] Silling S. Reformulation of elasticity theory for discontinuities and long-range force. Journal of the Mechanics and Physics of Solids, 2000, 48(1): 175–209

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

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

[41] Han F, Lubineau G, Azdoud Y, Askari A. A morphing approach to couple state-based peridynamics with classical continuum mechanics. Computer Methods in Applied Mechanics and Engineering, 2016, 301: 336–358

[42] Han F, Lubineau G, Azdoud Y. Adaptive coupling between damage mechanics and peridynamics: A route for objective simulation of material degradation up to complete failure. Journal of the Mechanics and Physics of Solids, 2016, 94: 453–472

[43] Sloan S W. Upper bound limit analysis using finite elements and linear programming. International Journal for Numerical and Analytical Methods in Geomechanics, 1989, 13(3): 263–282

[44] Sloan S W. Lower bound limit analysis using finite elements and linear programming. International Journal for Numerical and Analytical Methods in Geomechanics, 1988, 12(1): 61–77

[45] De Borst R, Vermeer P. Possibilities and limitations of finite elements for limit analysis. Geotechnique, 1984, 34(2): 199–210

[46] Le C V, Nguyen-Xuan H, Nguyen-Dang H. Upper and lower bound limit analysis of plates using FEM and second-order cone programming. Computers & Structures, 2010, 88(1–2): 65–73

[47] Leca E, Dormieux L. Upper and lower bound solutions for the face stability of shallow circular tunnels in frictional material. Geotechnique, 1990, 40(4): 581–606

[48] Mollon G, Dias D, Soubra A H. Rotational failure mechanisms for the face stability analysis of tunnels driven by a pressurized shield. International Journal for Numerical and Analytical Methods in Geomechanics, 2011, 35(12): 1363–1388

[49] Zhang C, Han K, Zhang D. Face stability analysis of shallow circular tunnels in cohesive frictional soils. Tunnelling and Underground Space Technology, 2015, 50: 345–357

[50] Anagnostou G, Perazzelli P. Analysis method and design charts for bolt reinforcement of the tunnel face in cohesive-frictional soils. Tunnelling and Underground Space Technology, 2015, 47: 162–181

[51] Han K, Zhang C, Zhang D. Upper-bound solutions for the face stability of a shield tunnel in multilayered cohesive-frictional soils. Computers and Geotechnics, 2016, 79: 1–9

[52] Chen Z, Hofstetter G, Mang H. A Galerkin-type BE-FE formulation for elasto-acoustic coupling. Computer Methods in Applied Mechanics and Engineering, 1998, 152(1–2): 147–155

[53] Dang H, Meguid M A. An efficient finite discrete element method for quasi-static nonlinear soil-structure interaction problems. International Journal for Numerical and Analytical Methods in Geomechanics, 2013, 37(2): 130–149 doi:10.1002/nag.1089

[54] Smith C C, Gilbert M. Application of discontinuity layout optimization to plane plasticity problems. Proceedings: Mathematical, Physical and Engineering Sciences, 2007, 463(2086): 2461–2484

[55] Smith C C, Gilbert M. Identification of rotational failure mechanisms in cohesive media using discontinuity layout optimisation. Geotechnique, 2013, 63(14): 1194–1208

[56] Gilbert M, Casapulla C, Ahmed H. Limit analysis of masonry block structures with non-associative frictional joints using linear programming. Computers & Structures, 2006, 84(13–14): 873–887

[57] Hawksbee S, Smith C C, Gilbert M. Application of discontinuity layout optimization to three-dimensional plasticity problems. Proceedings of The Royal Society A, 2013, 469(2155): 20130009

[58] Lin S, Xie Y, Li Q, Huang X, Zhou S. On the shape transformation of cone scales. Soft Matter, 2016, 12(48): 9797–9802

[59] Nanthakumar S, Lahmer T, Zhuang X, Park H S, Rabczuk T. Topology optimization of piezoelectric nanostructures. Journal of the Mechanics and Physics of Solids, 2016, 94: 316–335

[60] Ghasemi H, Park H S, Rabczuk T. A multi-material level set-based topology optimization of flexoelectric composites. Computer Methods in Applied Mechanics and Engineering, 2018, 332: 47–62

[61] Smith C C, Cubrinovski M. Pseudo-staticlimit analysis by discontinuity layout optimization: Application to seismic analysis of retaining walls. Soil Dynamics and Earthquake Engineering, 2011, 31(10): 1311–1323

[62] Clarke S D, Smith C C, Gilbert M. Modelling discrete soil reinforcement in numerical limit analysis. Canadian Geotechnical Journal, 2013, 50(7): 705–715

[63] Clarke S D, Smith C C, Gilbert M. Analysis of the stability of sheet pile walls using Discontinuity Layout Optimization. In: Benz T, Nordal S., eds. Numerical Methods in Geotechnical Engineering-Proceedings of the 7th European Conference on Numerical Methods in Geotechnical Engineering (NUMGE 2010). CRC Press, 2010, 163–168

[64] Gilbert M, Smith C C, Haslam I, Pritchard T J. Application of discontinuity layout optimization to geotechnical limit analysis problems. In: Benz T, Nordal S., eds. Numerical Methods in Geotechnical Engineering-Proceedings of the 7th European Conference on Numerical Methods in Geotechnical Engineering (NUMGE 2010). CRC Press, 2010, 169–174

[65] Zhang Y, Zhuang X, Lackner R. Stability analysis of shotcrete supported crown of NATM tunnels with discontinuity layout optimization. International Journal for Numerical and Analytical Methods in Geomechanics, 2018, 42(11): 1199–1216

[66] Gilbert M, Smith C C, Pritchard T J. Masonry arch analysis using discontinuity layout optimisation. Proceedings of the Institution of Civil Engineers-Engineering and Computational Mechanics, 2010, 163(3): 155–166

[67] He L, Gilbert M. Automatic rationalization of yield-line patterns identified using discontinuity layout optimization. International Journal of Solids and Structures, 2016, 84: 27–39

[68] He L, Gilbert M, Shepherd M. Automatic yield-line analysis of practical slab configurations via discontinuity layout optimization. Journal of Structural Engineering, 2017, 143(7): 04017036 doi:10.1061/(ASCE)ST.1943-541X.0001700

[69] Zhang Y. Multi-slicing strategy for the three-dimensional discontinuity layout optimization (3D DLO). International Journal for Numerical and Analytical Methods in Geomechanics, 2017, 41(4): 488–507

[70] Zheng H, Tham L, Liu D. On two definitions of the factor of safety commonly used in the finite element slope stability analysis. Computers and Geotechnics, 2006, 33(3): 188–195

[71] Farias M, Naylor D. Safety analysis using finite elements. Computers and Geotechnics, 1998, 22(2): 165–181

[72] Gilbert M, Smith C. Evaluating Displacements at Discontinuities within a Body. UK Patent GB 2442496, 2008

[73] LimitState Ltd. Limitstate: Analysis and Design Software for Engineers. Sheffield: The Innovation Centre, 2019

[74] Zhang Y, Pichler C, Yuan Y, Zeiml M, Lackner R. Micromechanics-based multifield framework for early-age concrete. Engineering Structures, 2013, 47: 16–24

[75] Zhang Y, Zeiml M, Pichler C, Lackner R. Model-based risk assessment of concrete spalling in tunnel linings under fire loading. Engineering Structures, 2014, 77: 207–215

[76] Zhang Y, Zeiml M, Maier M, Yuan Y, Lackner R. Fast assessing spalling risk of tunnel linings under RABT fire: From a coupled thermo-hydro-chemo-mechanical model towards an estimation method. Engineering Structures, 2017, 142: 1–19

[77] Zhang Y, Zhuang X. A softening-healing law for self-healing quasi-brittle materials: Analyzing with strong discontinuity embedded approach. Engineering Fracture Mechanics, 2018, 192: 290–306

[78] Sun Z, Li S, Rentai L, Bo G, Zhang Q, Lewen Z. Quantitative research on grouting reinforcement of soft fluid-plastic stratum. Chinese Journal of Rock Mechanics and Engineering, 2016, 35: 3385–3393 (in Chinese)

RIGHTS & PERMISSIONS

Higher Education Press and Springer-Verlag GmbH Germany, part of Springer Nature