Analysis of excessive deformations in tunnels for safety evaluation

Tunnelling and Underground Space Technology 45 (2015) 190–202

Contents lists available at ScienceDirect

Tunnelling and Underground Space Technology journal homepage: www.elsevier.com/locate/tust

Analysis of excessive deformations in tunnels for safety evaluation He Manchao a, R. Leal e Sousa b, André Müller c, Eurípedes Vargas Jr. c, L. Ribeiro e Sousa a,d,⇑, Chen Xin a a

State Key Laboratory for Geomechanics and Deep Underground Engineering, Beijing, China MASDAR Institute of Science & Technology, Abu Dhabi, UAE c Catholic University of Rio de Janeiro (PUC Rio), Rio de Janeiro, Brazil d University of Porto, Porto, Portugal b

a r t i c l e

i n f o

Article history: Received 10 March 2014 Received in revised form 18 August 2014 Accepted 20 September 2014

Keywords: Tunnels Excessive deformation Material point method Numerical methods Deep coal mines

a b s t r a c t Tunnel construction is increasing worldwide in mining and civil engineering. There have been several accidents that resulted in delays, cost overruns, some with more severe consequences. To help minimize these accidents, it is necessary to assess and manage the risks associated with tunnel construction and exploration. A particular type of accident, or undesirable event, which can occur during tunnel construction and operation, is associated with the occurrence of excessive deformations occurring inside the tunnel. This can happen due to deficient design, construction defects, and high in situ stresses or due to specific swelling and squeezing grounds. Deep coal mines where large deformations can occur during and after excavations due to the soft properties of the rock and the high in situ stresses, are particularly vulnerable to this type of event. The associated non-linear problems are related with geomechanical behavior of the rock mass, changes in the geometry of cavities and in some cases with developing surface contacts due to large strains. In this paper, the phenomena involved in large material deformations are analyzed in detail and the basic equations for the Chen’s large deformation theory are presented. The application of an FEM based method to simulate large material deformations, the Material Point Method (MPM) to the simulation of large deformation that occurs in tunnels when failure occurs, is also described. An application of MPM to the Jiahe Coal Mine, in China is presented, and the numerical results obtained with MPM compared with solutions using Chen’s large deformation theory. Safety considerations about the excessive deformation scenario in tunnels are drawn and a risk assessment methodology with special use of Bayesian networks is proposed. A simplified schematic example was presented for two case scenarios. Ó 2014 Elsevier Ltd. All rights reserved.

1. Introduction Most accidents and undesirable events are associated with uncertainties, often related to uncertainties in site ground conditions. It is important to develop risk analysis systems to avoid or to mitigate their occurrence. The associated risk has a complex nature resulting from the combination of two sets of factors, the events and their consequences, and the vulnerability factors that determine the probability of an event having certain consequences. Risk analysis illustrates the fact that decision-making must be based on a certain level of uncertainty (Einstein, 2002; Sousa, 2010). Until relatively recently, risk assessment and risk analyses have not assumed particular relevance when evaluating underground projects and other major geotechnical projects. This situation is ⇑ Corresponding author at: China University of Mining and Technology, Beijing, D11 Xueyuan Road, Haidian District, Beijing 100083, China. Tel.: +86 13683553486. E-mail address: [email protected] (L. Ribeiro e Sousa). http://dx.doi.org/10.1016/j.tust.2014.09.006 0886-7798/Ó 2014 Elsevier Ltd. All rights reserved.

changing and risk analysis has been successfully implemented, for major transportation infrastructures projects in UE, USA and Switzerland, using both commercial and research software for risk analysis. Special reference can be made to the system DAT – Decision Aids for Tunneling, developed by Massachusetts Institute of Technology (MIT), USA, in co-operation with Ecole Polytechnique Fédérale de Lausanne (Einstein, 2006). DAT is an interactive program that uses probabilistic modeling to analyze the effect of geotechnical uncertainties and construction uncertainties on construction costs and time through probabilistic modeling. Recent developments were introduced in this software allowing dealing with other structures like embankments, bridges and Engineered Geothermal systems (Moret et al., 2009; Sousa et al., 2012). A type of accident associated to tunnel construction and operation is the occurrence of excessive deformations inside the tunnel (He et al., 2012), often occurring due to particular type of ground, where swelling and squeezing happen. It is particularly important to consider the possibility of occurrence of excessive deformations in deep coal mines where large deformations may develop during

H. Manchao et al. / Tunnelling and Underground Space Technology 45 (2015) 190–202

and after excavations due to the soft properties of the rock and the high in situ stresses. In China, coal resources play a leading role in the energy strategy, and represent more than 70% of the total energy consumption. Deep coal resources below 1000 m represent the majority of resources, about 70% of the total. Fig. 1 confirms this evidence in the development trend of the average mine depth of the major coal mines in China. In 2010, the planned average depth was about 700 m and in 2020 is expected to have an average of 1200 m for the referred mines (He, 2006). As coal mines went deeper and deeper, accidents and their consequences became more severe as they occurred in conditions of large deformations often associated with rockbursts, waterbursts, high temperatures and gas outbursts. In this context, some of the phenomena involved are analyzed in detail and the basic equations are presented for the classic large deformation approaches and for Chen’s large deformation theory (Chen, 1979, 2000; He et al., 2012). The finite element method was applied to the different approaches for large deformation hazards, including non-linear geometrical approaches and commonly used rock mechanics constitutive models. Almost all existing commercial numerical software only include classical large deformation theories for these types of problems, in order to incorporate suitable models for geomaterials such as Mohr–Coulomb, Drucker–Prager, Hoek–Brown amongst others and supports. The paper also describes the use of the Material Point Method (MPM) to the simulation of large deformation that occurs in tunnels of deep coal mine, illustrated with a practical application to roadway tunnel in the Jiahe deep coal mine. The depths at which this roadway tunnel is located vary between 900 m and 1032 m. The application investigates the influence of several coal seams at different depths in the large deformation behavior of the cavities for different in situ state of stresses. Safety considerations regarding the excessive deformation scenario in tunnels are drawn and a risk assessment methodology using Bayesian networks (BN) is proposed (Sousa, 2012). Two different scenarios of probability of occurrence of large deformations are analyzed and conclusions regarding the occurrence of large deformation scenario are drawn. 2. Excessive deformations in tunnels 2.1. General The identification of the geotechnical risks aims to assess all causes that may lead to a hazard-inducing process. The study of

Fig. 1. Trends of average mining depth of major coal mines in China (adapted from He, 2006).

191

accidents in underground structures is very important in understanding the instability phenomena and mechanisms that are produced by underground construction and, as a result, allowing for the selection of the most appropriate construction methods for future projects. Even though the occurrence of accidents and incidents in underground works is not as unusual as in other geotechnical structures, their dissemination is not always common (Vlasov et al., 2001; Cost C7, 2003; Stallmann, 2005; Seidenfub, 2006; Sousa, 2006, 2010). In fact, there is a tendency to minimize their dissemination and causes. This fact may lead to the repetition of former errors by designers and constructors. Therefore, the number of failure reported in underground works is in comparison to those of other type of construction smaller. An important study on cases of accidents in tunnels during construction was conducted at MIT by Sousa (2010). In this study, data on accidents were collected from the technical literature, newspapers and correspondence with experts in the tunneling domain. The data were stored in a database and analyzed, and the accidents were classified into different categories, their causes and their consequences were evaluated. The main goal of the study was to determine the major undesirable events that may occur during tunnel construction, their causes, chain of events and consequences and ultimately present mitigation measures to avoid accidents on tunnels during construction. To gather information form experts in the field, a survey was created and send to experts worldwide. The survey consisted of two sections: Project Information and Accident Information. The Project Information section contained information related to the project (e.g. tunnel dimensions; geological and geotechnical information, construction method, among others), as well as the location in the tunnel of the accident. The Accident Information section is most important and collects information that is specific to the accident itself (Sousa, 2010). In total 204 cases were registered based on questionnaires responses, literature review and private correspondence with experts from 25 countries in Europe, North and South America and Asia. The study identified 9 types of events which were classified, as: Rock Fall; Collapse; Daylight Collapse; Large Deformations; Waterburst and Flooding; Rockburst; Specific Locations, such as portals and shafts; and Other Types that included slope failures. Fig. 2 shows the distribution undesirable events in the database. The great majority of events reported are Collapses and Daylight Collapses (36% and 28%, respectively). This does not mean that these are the great majority of events that occur during tunneling construction. They are however the most frequently reported in the literature and by the experts because most likely they are the ones with more severe consequences on the construction process, the safety of the workers and people and structures at the surface. The consequences reported ranged from loss of live, equipment damage and damage to the tunnel structure that may lead to a collapse. One of the frequent types of accidents occurring during construction is related with excessive deformations occurring inside the tunnel or at the surface. These events can be caused, for example, by deficient design, construction defects and/or due to particular type of terrains such as swelling and squeezing ground, which had not been predicted. Fig. 3 shows large deformations of the surrounding rock located in China. The serious distortion and rupture of steel occurred at the South water storage tunnel (He et al., 2012). The location of this type of scenario is in Fig. 4. As mentioned previously, one of the main causes of excessive deformation is crossing fault zones composed of squeezing and weak strata. Squeezing ground is characterized by excessive ground pressure that may lead to support failure and sometimes even cause collapse. It generally occurs around the whole cross section frequently involving the invert as well. The development of both rock pressure and rock deformations is time dependent.

192

H. Manchao et al. / Tunnelling and Underground Space Technology 45 (2015) 190–202

Fig. 2. Undesirable event distributions in the database (Sousa, 2010).

Fig. 3. The rupture of steel at the top in the South water storage tunnel (He et al., 2012).

Fig. 4. Location of excessive deformations in tunnels (Sousa, 2010).

This occurrence of large pressure may lead to failure of the lining and/or result in great difficulties for completing underground works, with major delays in construction schedules and cost overruns. Empirical data suggest that low strength, high deformability and the presence of water pressure facilitate squeezing (Kovari and Staus, 1996; Barla, 2001). Depending on the construction method used, the consequences of excessive deformation will be different as well as the mitigation measures that can be used to address this situation. If the tunnel is excavated by conventional means, the excessive deformation will usually cause failure and damage to the primary support, requiring re-excavation to the original tunnel profile (due to the reduction in cross-section) and replacement of the support in the affected section. The support options for tunnel in squeezing ground range from rock bolts (minor squeezing) to shotcrete with longitudinal slots (severe squeezing) that provide stress relieve caused by excessive deformations and allow for deformation to take place in a controlled manner. In the case of mechanized tunneling with a shield or TBM, a possible consequence is for the machine to get trapped during the drive, which in the worst case can lead to abandonment of the TBM. Another possible cause for excessive deformation is swelling. This phenomenon in tunnels is described as a time dependent volume increase of the ground, leading to inward movement of the tunnel perimeter. Some minerals exhibit the property of increasing their volume when absorbing water (swelling). If this volume increase is prevented (for example by a tunnel invert) then a corresponding pressure is exerted against the element preventing the movement. In tunnel construction swelling can cause a longterm heave of the tunnel floor, which can impair the serviceability of the structure, or result in partial failure of the lining (Kolymbas, 2005; Sousa, 2010). Einstein (2000) presents several case studies of tunnels excavated through Opalinus Clayshale and gypsum (Keuper) in the Swiss Jura Mountains, which show how problematic swelling can be during construction and also during operation, if the invert is not strong enough and if water reaches into the shale. Three typical swelling mechanisms have been identified (Einstein, 1996): (i) Mechanical swelling, occurring mostly in clays, silty days, clayey silts and corresponding rocks, caused by the dissipation of negative excess pore pressure. (ii) Osmotic swelling which occurs in argillaceous rocks, and is related to the large difference in concentrations between ions near the clay particles’ surface and those in the clay pore water (double layer effect). (iii) Intra crystalline swelling/hydration which occurs in smectite and mixed layer clays, in anhydrite and in pyrite and marcasite. The mechanisms involved depend on the type of material. Common to all three mechanisms is the important role of pore pressure in the phenomena of swelling. In order to predict the behavior of a tunnel on swelling or squeezing ground, it is necessary to know the natural stress state, stress changes, ground water conditions and material properties. In order to be able to make adequate predictions regarding this type of behavior, the engineer should perform several tests that will allow him to identify and quantify the swelling properties of the ground (Einstein, 1996; Barla, 2008). However, due to interaction of different mechanisms, it is not always very easy to predict the amount of swelling that may occur. Swelling occurs mostly in the tunnel invert, and can develop more or less rapidly depending on the access of water to the excavation. A case of swelling that occurred during tunnel construction is the one of the Chienberg tunnel in Switzerland, where during the time

H. Manchao et al. / Tunnelling and Underground Space Technology 45 (2015) 190–202

193

Fig. 5. Influence diagram for excessive deformation (Sousa, 2010).

that the tunnel construction was stopped due to a previous collapse, the invert was left open. After 4 weeks a heave of 1.5 m was observed in the invert near (behind) the zone of the collapse (Sousa, 2010). Fig. 5 presents an influence diagram summarizing conditions leading to excessive deformations and their consequences, in tunnel construction and exploration. Recommendations for the construction include yielding support, lining with longitudinal slots to allow excessive deformation, and also in some cases rigid supports. For monitoring it is recommended measuring of convergences and use of inclinometers, extensometers, strain gages, load cells, among others. 2.2. Large deformations in deep coal mines

Fig. 6. Characteristics of rock mass geomechanical properties at great depth (He et al., 2012).

Comprehensive investigations have been performed in China on the problem of large deformations in deep mines, with special emphasis in coal mines (He, 2006). According to the investigations, the deformation mechanism of rock mass at great depth varied under different confining pressures, with time effect being of significant to the behavior of rock masses subjected to high confining pressures. In general rock strength is stronger with depth, and the mechanisms of failure, which vary also with depth are, in deep mining controlled by the confining stress. Deep rock mass is a very complex material, which depends on long term rheological rock mass characteristics. Fig. 6 illustrates the different factors that influence the properties of the rock mass. In deep mining activities the major problems are associated with large deformations and overstressing of the rock mass caused by excavations at great depth, and also with rockburst. In the case of large deformations hazards in deep coal mines it is necessary to have a comprehensive investigation of the involved

phenomena due to the soft properties of the rock and the high in situ stresses. The non-linear problems associated are mainly related with geomechanical behavior of the rock mass, the changes in the geometry of the cavities. According to statistics, in China, for large state-owned coal mines, the portion of the mines at a depth between 400 and 800 m are nearly 47%, those between 800 and 1000 m are close to 24%, and those over 1200 m are about 5% (He et al., 2012). Among them, 47 mines reach 1000 m and 30 mines are over 1200 m. Fig. 7 illustrates a severe floor heave of a roadway in China, in which the largest floor heave reached 1 m. A challenging issue, often encountered in deep mining engineering and design, is how to avoid or prevent the occurrence of catastrophic events, many related to large deformations. Consequently, a comprehensive investigation of deep mining mechanics is of interest on in situ monitoring, mechanical performance and deformation characteristics, strength, and post-failure and failure

194

H. Manchao et al. / Tunnelling and Underground Space Technology 45 (2015) 190–202

 where uj i is the covariant derivative of the displacement component uj with respect to co-moving coordinate xi. Then S–R decomposition theorem stated as any invertible linear differential transformation F has unique additive decomposition:

F ¼ S þ R;

F ij ¼ Sij þ Rij ;

ð2Þ

where S is a symmetrical and positive definite sub-transformation representing strain tensor,

Sij ¼

 1 i u jj þ ui jTj  ð1  cos HÞLik Lkj 2

ð3Þ

R is an orthogonal sub-transformation, representing local mean rotation tensor,

Rji ¼ dji þ Lij sin H þ ð1  cos HÞLik Lkj ; Fig. 7. The severe floor heaves of the roadway (He et al., 2012).

mechanism for rock mass which is embedded at deep level (He, 2006). In these deep mines soft rock engineering problems occurred frequently. The geological age of these deep coal strata changes from Paleozoic Carboniferous-Permian to Cenozoic Tertiary. Some examples of large deformation phenomena in coal mines in China are presented and grouped in accordance with their geological age (see He et al., 2012 for more details). 3. Numerical simulation of large deformation problems Numerical simulation of large strain problems in engineering can be rather complex as it leads to intrinsically non-linear problems. These non-linearties are related to outstanding geometry changes including contact development (geometrical non-linearities) and non-linear material (constitutive) properties as occurring in non-linear-elastic and elasto-plastic conditions (material nonlinearities). In the present work all of the above mentioned conditions, which occur in practice, have been taken into account in the analyses carried out using two computer codes, one based on Chen’s large deformation theory and the other on the Material Point Method (MPM), which are described below.

Here, H is mean rotation angle, and Lij is unit vector of rotation axis direction (He et al., 2012). Since the final result of nonlinear finite deformation and finite rotation depends on the stress and strain path, so in the most cases, incremental updated method should be used, in which stress and strain are defined in real time deformed state. In the updated co-moving coordinate txi, i = 1, 2, 3, the strain rate and mean solid rotation speed are given by:

  T 1   V i þ V i ; S_ ij ¼ j j 2

  T i _ ¼ 1 V i R_ ji ¼ Lij H  V  : j j 2

ð5Þ

  where V j  is the covariant derivative of the speed component Vj i with respect to the updated co-moving coordinate. In an updated co-moving coordinates, material derivative of Euler stress tensor is

Drij Dt

  ¼ r_ ij  rik v k j  rij v k k

ð6Þ

which is not an objective derivative respect to time and is related to rigid rotation. The objective derivative of Euler stress is defined by D

rij ¼ r_ ij  rik S_ kj  rij S_ kk

ð7Þ

Equilibrium equation is established in real time deformed state:



rij  þ f j ¼ 0: i

ð8Þ

where fj is the volume force. The rate form constitutive relation is given by:

3.1. Chen’s large deformation theory For large deformation problems, the theories of classical large deformation and Chen’s large deformation have been developed (Truesdell and Noll, 1965; Chen, 1979, 2000; He et al., 2012). More recently, Sulsky and Schreyer (1996) proposed the Material Point Method (MPM), which represents the material contained in a region as a collection of unconnected particles. The fundamentals of the MPM will be presented in detail in Section 3.2. The classical theory of large deformations is based on the use of the Green strain tensor and on the rotation tensor defined by Finger-Truesdell’s polar decomposition theorem and details of the formulation is presented at He et al. (2012). The Chen’s large deformation theory based on S–R decomposition theorem describes the general motion of a deforming body using a co-moving coordinate system method (Chen, 1979, 2000). Suppose a body B0 is identified by a set of co-moving coordinates {x1, x2, x3}, and transforms to B at time t, and the basic vectors   0 0 0 g1 ; g2 ; g3 at a point change to {g1, g2, g3}. The linear differential transformation is denoted by deformation gradient F

B0 fx1 ; x2 ; x3 ; t 0 g ! Bfx1 ; x2 ; x3 ; tg;  0 gi ¼ F ji gj ; F ji ¼ dji þ uj i ;

ð4Þ

ð1Þ

D

_k rij ¼ Di...l jk Sl

ð9Þ

i...l Djk

where is the tangential stiffness tensor of the material. For the application of the finite element method three types of nonlinearity were considered related to large deformations, with geometrical and contact problems, and physical properties of rock mass and supports and considering a two-dimensional approach (He et al., 2007). A software was developed called LDEAS (Large Deformation Engineering Analysis Software) at China University of Mining and Technology (CUMT), Beijing. The theory and the fundaments of the application of the finite element method are described in detail in He et al. (2012). 3.2. Application of material point method The MPM is an extension of the fluid implicit particle method (FLIP) developed by Brackbill and Ruppel (1986). FLIP, on the other hand, is an extension of the particle in cell method, Harlow (1964). The development of the MPM was introduced by Sulsky et al. (1994, 1995). An increasing number of cases of application of the method in the field of geotechnical engineering can be found in the technical literature. An interesting application of this method

H. Manchao et al. / Tunnelling and Underground Space Technology 45 (2015) 190–202

was made in the study of the failure of Aznalcóllar dam (Zabala and Alonso, 2012) due to earthquake activity. The MPM is a combination between the Lagrangian and Eulerian methods. The method uses Lagrangian material points to carry history-dependent properties and an Eulerian background mesh to calculate derivatives and solve the equations of motion. As the dynamic analysis proceeds, the solution is tracked on to the material points by updating all required properties such as position, velocity, acceleration, strain, stress and others. At each time step, the material point information is extrapolated to a fixed background mesh to solve the equations of motion. Once the equations are solved, the background mesh solution is used to update all material points’ properties. A typical 2D MPM problem setup is shown in Fig. 8. The dotted line represents the boundary of the simulated media X and each closed point in this media represents a material point used to discretize X. The square mesh represents the background mesh. The background mesh nodes are located at the corners of mesh cells. The combination of Lagrangian and Eulerian methods has been proven useful for solving solid mechanics problems including those with large deformations or rotations and involving materials with history dependent, plasticity and viscoelasticity for example. More recently, the MPM has been used to solve fluid mechanical coupling problems (Zhang and Wang, 2009). In the paper a continuum X is considered. This continuum subject to only mechanical conditions is governing by differential equations, firstly by the mass conservation equation

dq þ qdivðv Þ ¼ 0 dt

ð10Þ

integrated over the domain X and subjected to Green´s theorem (see Dow, 1998 or Bathe, 1996). By discretization of the problem, the continuum X is divided into a finite number of material points. Following the discretization, each material point is assigned a mass, velocity, acceleration, stress, strain and other relevant properties defining the state of the material point. In the sequence, the main steps of the MPM are presented. More details can be found in the references cited above. The fixed background mesh is described by the coordinates of the nodes xi, the material points are described by coordinates of the materials points xp, N represents the interpolation functions, P all material points in a determinate cell p is the summation over P of the background mesh, i represents the summation over all nodes of the cell where the material point is located, Nip is the interpolation function to node i of the background mesh evaluated at the material point position, Fp is the deformation gradient at the material point and I is the identity matrix. As previously mentioned, in the initial step of the MPM, each material point is assigned a mass, velocity, acceleration, stress, strain and other relevant properties defining the state of the material point. In this paper, the initial stress conditions are applied in all points of continuum thus as different material properties and the boundary conditions to define to problem. After the initial descriptions, for one time step of the simulation is computed de material points volumes t V p ¼ detðt Fp Þt¼0 V p and are computed the mass, Eq. (12), momentum, Eq. (13), velocity, Eq. (14) internal forces, Eq. (15) and accelerations, Eq. (16), all in the background mesh nodes: t

and also by the momentum conservation equation

195

mi ¼

X Nip t mp

ð12Þ

p

q

dv  divðrÞ  qb ¼ 0 dt

ð11Þ

In Eqs. (10) and (11), q = q(x, t) is the density, v = v(x, t) is the velocity, r = r(x, t) is the Cauchy stress tensor and b = b(x, t) is the specific body forces, all in a x position in time t. Eqs. (10) and (11) are supplemented with a suitable constitutive equation and kinematic relations between strain, stress and displacements. In MPM the conservation of mass is naturally satisfied and the momentum conservation equation is solved by discretization in a weak form. In order to obtain the weak form of the momentum conservation equation, it is multiplied by an arbitrary test function,

X t pi ¼ Nip mp t v p

ð13Þ

p t

v i ¼ t pi =mi

t int fi

ð14Þ

X ¼  rNip t rp t V p

ð15Þ

p t

int

ext

ai ¼ ðt f i þ t f i þ t jt pi Þ=mi

ð16Þ

ext

In Eq. (16), t f i is the set of the external forces acting in continuum at time t, body forces and surface forces. Also, in Eq. (16), t j represents the damping coefficient at time t, in this paper the damping coefficient is update in accord with kinetic energy variation. After this step, the velocities in the background mesh nodes are updated with tþDt v i ¼ t v i þ t ai Dt and the background mesh solution is used to update all material points’ properties, and for that purpose the velocity gradient of the all material points are computed as shown in Eq. (17). Subsequently, the deformation gradients, velocities, positions and densities are updated for all material points through Eqs. (18)–(21) respectively:

rtþDt v p ¼

X

rNip tþDt v i

ð17Þ

i tþDt tþDt

tþDt

DFp ¼ ðI þ rtþDt v p DtÞ X v p ¼ t v p þ Nip t ai Dt xp ¼ t xp þ

tþDt

i X Nip tþDt v i Dt

Fp ¼ tþDt DFp t Fp

ð18Þ ð19Þ ð20Þ

i tþDt

Fig. 8. A typical 2D MPM problem setup.

qp ¼

1 t¼0 qp detðtþDt Fp Þ

ð21Þ

With the updated properties of material points, the new strain and corresponding stress fields can be computed for the same material points. At this point, an appropriate stress–strain or stress–strain rate can be used. Therefore, the method allows with relative ease, the introduction of any constitutive model in the analysis. Updated strain and stress field are described by Eqs.

196

H. Manchao et al. / Tunnelling and Underground Space Technology 45 (2015) 190–202

(22) and (23) respectively. In Eq. (22), the increments of the total strain and increments of the plastic strain are presented, where k_ is a positive scaling factor and g is the plastic potential function, see Simo and Hughes (1997) for more details. In this paper an elastoplastic constitutive model is considered with the failure criterion of Mohr–Coulomb. In Eq. (24) the increment of the second Piola Kirchhoff stress tensor is presented, where C is the constitutive tensor, in this paper representing a plane strain behavior. Eq. (24) presents the update of the second Piola Kirchhoff stress tensor. The Cauchy stress tensor in end of the time step is given by Eq. (25). tþDt _

ep ¼

X1 i

tþDt

S_ p ¼ C



2

rNip tþDt v i þ tþDt v Ti rNTip

tþDt _

ep  tþDt e_ pp



e ¼ k_

tþDt _ p p



@g @ rp

Sp ¼ t Sp þ tþDt S_ p Dt 1 tþDt tþDt rp ¼ DFp tþDt Sp tþDt DFTp detðtþDt DFp Þ tþDt

ð22Þ ð23Þ ð24Þ ð25Þ

Once these calculations are completed, a new time step cycle can be started. 4. Numerical analysis in tunnels of Jiahe coal mine 4.1. Large deformation problems in Jiahe mine Xuzhou mining area is an old mining district with 120 years of history (Xi, 2009). Currently, the deepest excavation depth for 8 of 11 mines in Xuzhou mining area are greater than 1000 m, examples are the 1032 m deep Quantai mine, 1032 m deep Qishan mine, 1062 m deep Pangzhuang mine, 1200 m deep Jiahe mine, 1017 m deep Sanhejian mine, 1010 m deep Chacheng mine, 1260 m deep Zhangji mine and 1038 m deep Zhangshuanglou mine (He et al., 2012). Jiahe mine lies 4.5 km away from the northwest of Xuzhou city in Jiangsu Province. The mine field has an area of 24.75 km2, and an annual production capacity of 3,000,000 tons. It has six exploration levels, from 280 until and 1200 m. The gateway at no. 2442 work face from West mining area in Jiahe is at a depth between 850 and 900 m. The coal seam in exploration is fragile blocky rock mass and developed with many fractures. The main roof is a composite stratum with coal and other weak interlayer. The immediate floor is gray-black sandy shale with many fractures, with the thickness of 1.7 m. The main floor is sandstone with the thickness of 10 m. X diffraction experiment found that the contents of clay mineral in roof and floor strata are 50% and 52%, respectively. The high expansion minerals, montmorillonite and illite–smectite are 30%, and microscopic cracks developed largely. The gateway was originally supported by bolts and shotcrete with wire mesh. During 75 days of observation, the average amount of roof subsidence, sidewall convergence and floor heave are 42, 56 and 75 cm, respectively. The length of the gateway that has large sidewall convergence is about 40% of the total length. Large deformation happened that were contained by the wire mesh. However in many sections columns and trapezoidal steel supports had been used. The high strength wire mesh had large deformation which causes the pull-out of many bolts. About 40 bolts were fractured, many cables dropped and the wire mesh on the roof was severely deformed (Figs. 9 and 10), (He et al., 2012).

Fig. 9. Large deformation of the Jiahe coal mine gateway.

model taking a 2D region, 30 m wide and 30 m high was used in both analysis. The cross section of the tunnel has a trapezoidal form with a width of 4.1 m and the heights of the left and right sidewalls of 2.8 and 2.0 m, respectively. A depth of 1500 m was considered for the MPM numerical model where the vertical stress was calculated for an average specific weight of 27 kN/m3 and the horizontal stress was determined by using k0 = 0.5. For the LDEAS model, the depths of 500, 1000 and 1500 m were considered, with values of k = 0.5 and 1.0. Fig. 11 shows the conceptual numerical model, being the geomechanical formations indicated. The excavation process was simulated in one step only. In order to model the excavation process, the material of the tunnel is included in the beginning of the analysis and removed afterwards. The displacements, in the whole of the boundary of the model were fixed. The material parameters of rocks are also listed in Fig. 11. In table of the figure, k is the specific weight, E is Young’s modulus, m is the Poisson radio, c is the cohesion and / the friction angle. The following sequence of analyses was followed. Initially, analysis has been carried out in order to validate the developed code using MPM. In this case, the analyses carried out comprise small strain scenarios. In the sequence, analysis with LDEAS was carried out where large strains are considered under (non-linear) elastic conditions. Finally, analyses were carried out with MPM under large strains and elasto-plastic conditions. The results are compared and discussed. For the initial set of analyses using MPM, horizontal and vertical displacements with linear elastic materials (small strain scenarios) are illustrated in Figs. 12 and 13, respectively. In both figures, a comparison is made with a classical FEM numerical solution (Fig. 12b) showing a good agreement of both solutions. These

4.2. Numerical modeling of deep tunnels Large deformation of a tunnel in Jiahe coal mine was analyzed using both the MPM method at the Catholic University in Rio de Janeiro and the LDEAS software developed at CUMT. The same

Fig. 10. Wire mesh severed deformed at Jiahe coal mine.

197

H. Manchao et al. / Tunnelling and Underground Space Technology 45 (2015) 190–202

Rock

γ (kN/m3)

E ν (GPa)

c (MPa)

φ ( 0)

Finestone

26.5

8.7

0.21

3.8

33.0

Sandy mudstone

20.0

6.6

0.22

2.0

24.0

Sandy shale

26.4

6.2

0.23

1.0

28.0

Coal

16.5

4.9

0.23

0.08

20.0

Coal

16.5

4.0

0.23

0.08

20.0

Sandy shale

26.3

5.7

0.24

1.2

30.0

Sandy mudstone

20.0

6.2

0.23

2.0

24.0

Sandstone

26.2

9.8

0.23

3.0

35.0

Fig. 11. Conceptual model and geomechanical formations.

Table 1 Deformations in points of surrounding rocks for the tunnel at depth 1500 m and k = 0.5. Displacements (m)

Sidewall shrinkage Floor heave Roof subsidence

Large deformation modules Chen’s theory

Classical theory (polar)

0.216 0.423 0.326

0.475 0.699 0.442

results (Figs. 12 and 13) are basically presented here as a validation exercise for the developed computer code using MPM. Using the LDEAS model, different analyses were performed and the corresponding results have already been described in He et al. (2012). Due to instabilities in the calculations, only elastic materials (under large deformations) were considered. Results considering displacements of the tunnel vs. tunnel depth for k = 0.5 are indicated in Fig. 16, for the roof subsidence, floor heave and sidewall shrinkage. Table 1 compares displacements at 1500 m and k = 0.5 using Chen´s theory and classic theory. Finally a comparison between the results obtained by MPM method and the Chen’s theory is illustrated in table.

Again using MPM and considering at this stage materials with elastoplastic behavior as indicated in Fig. 11, horizontal and vertical displacements are illustrated in Figs. 14 and 15, respectively, for a scenario of large strains. In both figures, results are presented for two time stages of the MPM analysis, intermediate and final instants. The figures clearly show the deformation patterns to be expected in such conditions, adding valuable information for risk analysis assessment as described in the next section. From the numerical point of view it is worth mentioning that the use of MPM can be very convenient as it does not require remeshing no matter the level of deformation attained. Table 2 allows a comparison between linear elastic (small deformation), non-linear elastic (large deformation) and non-linear elasto-plastic formulations. The results make clear that nonlinear effects both geometrical and material should be considered in predicting the evolution of deformation patterns encountered in situations such as in Jiahe coal mine. The differences between displacements in non-linear elastic analysis under large strains (obtained with LDEAS) and displacements obtained in elasto-plastic analysis under large strains (obtained with MPM) are also noteworthy. The use of non-linear elastic materials in the analyzed cases, even under large strain conditions would considerably underestimate the predicted displacements. In this case it would

Fig. 12. Displacements in horizontal direction (a) MPM, (b) FEM, scale in m.

198

H. Manchao et al. / Tunnelling and Underground Space Technology 45 (2015) 190–202

Fig. 13. Displacements in vertical direction (a) MPM, (b) FEM, scale in m.

Fig. 14. Displacements in horizontal direction for two time steps of the MPM analysis, scale in m.

Fig. 15. Displacements in vertical direction for two time steps of the MPM analysis scale in m.

H. Manchao et al. / Tunnelling and Underground Space Technology 45 (2015) 190–202

(a) S-R decomposition module

(b) Polar decomposition module Fig. 16. Displacements of the tunnel vs. tunnel depth for k = 0.5 (He et al., 2012).

Table 2 Comparison of displacements by MPM and Chen’s theory. Displacements (m)

Sidewall shrinkage Floor heave Roof subsidence

MPM

Chen’s theory

Elastic

Elasto-plastic

0.011 0.030 0.029

1.740 1.190 0.960

0.216 0.423 0.326

be important to take into account the possible failure of materials (elasto-plastic conditions) under large strains.

5. Safety considerations about excessive deformation scenario Accidents and other associated problems occur during the construction and exploration of the underground structures and are very often related to uncertainties related to the ground conditions. To help eliminate or at least reduce these accidents, it is necessary to systematically assess and manage the risks associated with construction and operation of underground structures. Traditionally, risks have been managed indirectly through the engineering decisions taken during the project development. However, the guidelines from ITA (Eskesen et al., 2004) consider that present risk management processes can be significantly improved by using systematic risk management techniques throughout the underground project development. By the use of these techniques potential problems can be clearly identified such that appropriate risk mitigation measures can be implemented in a timely manner (Einstein, 2002; Brown, 2012; Sousa, 2012). Risk assessment and management requires an evaluation of the hazard and the likelihood of the harmful effects. It starts with the

199

hazards identification, focusing on the likelihood of damage extent, followed by risk characterization, which involves a detailed assessment of each hazard in order to evaluate the risk associated to each one of them (Sousa, 2010). Risk assessment is developed with the goal of avoiding major problems that can occur. Although there are different techniques to represent knowledge and to perform risk assessment (such as knowledge based systems, events and decision trees), in this paper the proposed risk assessment and decision support framework will be based on Bayesian networks and their extension to influence diagrams. There are several advantages that BN have over other methods. A BN, also known as belief network, is a graphical representation of knowledge for reasoning under uncertainty. Over the last decade, BN have become a popular model for encoding uncertain expert knowledge in expert systems (Heckerman, 1997). BN can be used at any stage of a risk analysis, and may substitute both fault trees and event trees in logical tree analysis. While common cause or more general dependency phenomena pose significant complications in classical fault tree analysis, this is not the case with BN. They are in fact designed to facilitate the modeling of such dependencies. Because of what has been stated, BN provide a good tool for decision analysis, including prior analysis, posterior analysis and pre-posterior analysis. Furthermore, they can be extended to influence diagrams, including decision and utility nodes in order to explicitly model a decision problem. A methodology for risk assessment and decision for tunnels making using BN was proposed by Sousa (2012) for risk assessment and decision making for tunnel projects during design and construction phases. During the design phase, different construction strategies can be defined and evaluated. The tunnel alignment is then divided into different sections of more or less homogeneous ground conditions, and each section is treated independently. For each section the following information is needed: (i) prior probability of geological states; (ii) vulnerability, i.e. the probability of failure mode k, in geology i with construction strategy j; (iii) the consequences of failure mode k, expressed in utilities; and (iv) cost of changing (where relevant) construction strategies. During the design phase, information is available regarding geological, hydrological conditions, as well as regarding construction method costs and times. This information is used to determine, for the different possible alignments the ‘‘optimal’’ construction strategy for each alignment. Most existing tools, determine the ‘‘optimal’’ construction strategy in terms of costs and time. However, in the context of this study the main focus will be to determine the ‘‘optimal’’ construction strategy (or method) in terms of risk of an undesirable event for each given alignment. Once the construction strategy is chosen, the construction phase starts. During construction information becomes available regarding the geological conditions crossed by the tunnel or underground structure, behavior of the excavation as well as information on the construction. This information should and must be used to update the predictions made during the design phase. In the context of the developed methodology emphasis will be placed on updating the geological conditions for the part of tunnel that has not been excavated based on the geological conditions encountered during excavation. This will then be used to update the construction strategy for the remaining unexcavated part of the tunnel. The risk of large deformations can also be analyzed by using a similar approach based on BN probabilistic tools. The referred methodology proposed by Sousa (2012) can be implemented for different hazards including, the occurrence of large deformations. The decision support framework can consist of four models (Fig. 18): the first one, a geological prediction model, would predict the geological conditions ahead of the excavation; the second, a large deformation prediction model would predict the

200

H. Manchao et al. / Tunnelling and Underground Space Technology 45 (2015) 190–202

Fig. 17. Decision support framework for design and construction phases.

strategies and associated utilities; and finally the fourth model, a decision model, one choose the construction methods that minimizes the risk. This methodology was successfully applied to tunnels of Porto Metro where three collapses occurred during construction (Sousa and Einstein, 2012). A simplified example that fits in the decision support framework, presented in Fig. 17, is used in this section to illustrate the use of BN and Monte Carlo Analysis, coupled with numerical calculations, a schematic of that procedure is presented below in Fig. 18. It is a multi-step analysis consisting of a Monte Carlo simulation and BN extended to an influence diagram. The modeling steps are as follows:

= deformaon

P (y>y*)

Step 1: Create a numerical model, using MPM method for example. Step 2: Generate a set of random inputs, xi1, xi2, . . ., xiq, based on probability distribution of the parameters. These can be any geomechanical characteristics, such as module of deformability E, cohesion c, or friction angle /. Step 3: Evaluate the model and store the results as yi. For example roof subsidence, floor heave and sidewall shrinkage. Step 4: Repeat steps 2 and 3 for i = 1 to n, being n the number of Monte Carlo simulations. Step 5: Analyze the results using histograms, summary statistics, confidence intervals, among others. Step 6: Use the results to determine the probability of deformations in a critical point of the tunnel cross section exceeding a certain value y⁄ (i.e. that the probability of excessive deformation). Step 7: Define a BN/influence diagram that represents the decision problem that one is faced with. Step 8: Use the probability determined in Step 6 as an input to the BN defined in the previous step.

Fig. 18. Decision support system workflow schematic (example).

deformations due to construction (to include numerical modeling results); the third one, a construction prediction model would simulate the construction and establish different construction

The initial and conditional probabilities, costs and utilities associated with the BN represented in Fig. 18, are shown in Fig. 19 (the values were chosen for illustrative purposes only). The BN contains three change nodes, one decision node and two utility nodes, representing:

201

H. Manchao et al. / Tunnelling and Underground Space Technology 45 (2015) 190–202

Repair Support Yes No -2000 0

Yes No

Excessive Deformaon Yes

No

0.9 0.1

0.05 0.95

Repair Support Yes Support Damage

Collapse

Support Damage Yes No

No

Yes

No

Yes

No

Yes

0.02

0.02

0.6

0.1

No

0.98

0.98

0.4

0.9

Collapse Yes No -10000 0 Fig. 19. Initial, conditional probability and utility tables.

Fig. 20. Results of the BN represented in Fig. 19 for two case scenarios: (a) P (excessive deformation) = 0.01; (b) P (excessive deformation = 0.70).

5.1. Change nodes 5.3. Utility nodes  Excessive deformation which represents the probability of excessive deformation occurring during construction/operation of the underground structure.  Damage Support which represents the probability of damage occurring in the support given the occurrence of excessive deformation.  Collapse which represents the probability of the structure collapsing given that the support damage, as well as the decision to reinforcing or not the support. These probabilities can be obtained based on numerical or mechanical models, as well as based on expert opinion. 5.2. Decision nodes  Repair/Reinforce Support which represents the decision of repairing or reinforcing structure the support.

 Utilities which represent the utilities associated with the collapse of the structure. It can be expressed in utilities or monetary value.  Cost which represents the cost (or utility) associated with repairing or reinforcing structure the support. Imagine two different case scenarios, in the first one the probability of excessive deformation is 0.02 and the second one this probability equals 0.70 (obtained through Monte Carlo simulations, using the method described previously). The results of these two scenarios are presented in Fig. 20. In the first situation the best decision (i.e. the one with the highest utility associated, 1335) is not to repair or reinforce the support. In the second case, with the high probability of excessive deformations the best decision is to repair or reinforce the support (i.e. the one with the highest utility associated 2200).

202

H. Manchao et al. / Tunnelling and Underground Space Technology 45 (2015) 190–202

Note that the BN can also be used in combination with a deterministic numerical calculation. In that case the obtained deformation value with be compared against a reference value y⁄ and based on that comparison determine whether large deformation occurs or not, which would be fed into the BN as evidence and propagate to determine the ‘‘best’’ decision. Note that even in cases of large deformation one may not have damage or reach collapse. 6. Conclusions The present paper presented a framework for the risk analysis in deep underground excavations associated with large deformations. The framework associates appropriate numerical techniques for analysis of large deformations in scenarios of deep excavations with Bayesian based decision support techniques. The occurrence of large deformations reported in deep underground structures, such as in mining requires the design of the appropriate support systems to minimize the consequences. A study was performed about the importance of the large deformation scenario from a database developed at MIT that shows the distribution of this scenario and their consequences. Comprehensive investigations in China have been performed for deep coal mines. The importance of the numerical simulation of large deformation problems in tunnels was emphasized in order to help in the prediction of this type of accident. For the simulation of large deformation problems, two different formulations were used: Chen’s theory and the MPM method. The simulation exercise related to the conditions encountered in Jiahe coal mine clearly shows the relevance of numerical analysis for (large) deformation predictions. There is an intrinsic risk associated with tunnel construction because of the limited a priori knowledge of the existing subsurface conditions. For this reason, it is therefore important to systematically assess and manage the risks associated with tunnel construction or mining. A decision support framework consisting of different prediction and decision models, based on BN was presented in this paper. This methodology would allow the incorporation of uncertainties (geological and from construction) and numerical modeling into the modeling and decision processes, ultimately make decisions regarding construction strategies based on the minimization of risk. This methodology was successfully applied to tunnels of Porto Metro where three collapses occurred during construction. A simplified example is presented to illustrate the applicability of the methodology to cases of excessive deformation, by combining appropriate numerical models, Bayesian networks and Monte Carlo simulations, allowing the quantification of uncertainty by including probability distributions of ground parameters into numerical calculations, via Monte Carlo simulations, determine hazard probability (i.e. probability of large deformation) and feed that into a decision model, based on Bayesian networks, that will allow one to determine the ‘‘best’’ mitigation measure (in this specific case the ‘‘best’’ decision regarding ground support remedial). References Barla, G., 2001. Lessons learnt from the excavation of a large diameter TBM tunnel in complex hydrogeological conditions. In: GeoEng 2000, International Conference on Geotechnical & Geological Engineering, Melbourne, Australia, 19–24 November. Barla, M., 2008. Numerical simulation of the swelling behaviour around tunnels based on triaxial tests. Tunn. Undergr. 23, 508–521. Bathe, K.J., 1996. Finite Element Procedures in Engineering Analysis. Prentice Hall. Brackbill, J.U., Ruppel, H.M., 1986. FLIP – a method for adaptively zoned, particle in cell calculations in 2 dimensions. J. Comput. Phys. 65, 314–343. Brown, E.T., 2012. Risk assessment and management in underground rock engineering – an overview. J. Rock Mech. Geotech. Eng. 4 (3), 193–204.

Chen, Z.D., 1979. Geometric theory of finite deformation mechanics for continuum. Acta. Mech. Sin. 2, 107–117 (in Chinese). Chen, Z.D., 2000. Rational Mechanics. Chongqing Publication, Chongqing (in Chinese). Cost C7, 2003. Avoiding Damage Caused by Soil–Structure Interaction: Lessons Learnt from Case Histories. Thomas Telford, Ed. Kastner et al., London, 77p. Dow, J., 1998. Unified Approach to the Finite Element Method and Error Analysis Procedures. Academic Press. Einstein, H., 1996. Tunnelling in difficult ground – swelling behaviour and identification of swelling rocks. Rock Mech. Rock Eng. 29 (3), 113–124. Einstein, H., 2000. A review in Opalinus clayshale. Rev. Case Hist.. Einstein, H., 2002. Risk assessment and management in geotechnical engineering. In: 8th Portuguese Geotechnical Congress, Lisbon, pp. 2237–2262. Einstein, H., 2006. DAT system. Geotechnical Risks in Rock Tunnels, Francis & Taylor, Matos, Sousa, Kleberger & Pinto (Eds.), London. Eskesen, S., Tengborg, P., Kampmann, J., Veicherts, T., 2004. Guidelines for tunnelling risk management: ITA, Working Group No. 2. Tunn. Undergr. Space Technol. 19, 217–237. Harlow, F.H., 1964. The particle-in-cell computing method for fluid dynamics. Method Comput. Phys. 3, 319–343. He, M., 2006. Rock mechanics and hazard control in deep mining engineering in China. In: Leung, Zhou (Eds.), 4th Asian Rock Mechanics Symposium, Singapore, 18p (in CD-Rom). He, M., Chen, X., Liang, G.P., Qian, H.S., Zhou, Y.F., Zhuang, X.Y., 2007. An introduction to version 1.0 of software on large deformation analysis for soft rock engineering at great depth. In: 11th ISRM Congress, Lisbon, Eds. Sousa, Olalla and Grossmann, pp. 169–173. He, M., Chen, X., Zhang, G., Sousa, L.R., 2012. Large deformation analysis in deep coal mines in China. In: Sousa, Vargas, Fernandes, Azevedo (Eds.), Innovative Numerical Modelling in Geomechanics, Chapter 18, CRC Press, pp. 333–353. Heckerman, D., 1997. A tutorial on learning with Bayesian networks. Data Min. Knowl. Disc. 1, 79–119. Kolymbas, D., 2005. Tunneling and Tunnel Mechanics. A Rational Approach to Tunneling. Ed. Springer, 437p. Kovari, K., Staus, J., 1996. Basic considerations on tunnelling in squeezing ground. Rock Mech. Rock Eng. 29 (4), 203–210. Moret, Y., Sousa, R.L., Einstein, H.H., 2009. The Decision Aids for Tunneling (DAT) – Introduction and Recent Application in Portugal. Workshop on Tunnels for High Speed Railways, Porto, p. 17. Seidenfub, T., 2006. Collapses in Tunneling. MSc Thesis, University of Stuttgart, EPFL and ITA. Simo, J.C., Hughes, T.J.R., 1997. Computational Inelasticity. Springer – Verlag, New York. Sousa, L.R., 2006. Learning with accidents and damage associated to underground works. In: Kleberger, Pinto (Eds.), Geotechnical Risks in Rock Tunnels, Francis & Taylor, Matos, Sousa, London, pp. 7–39. Sousa, R.L., 2010. Risk Analysis for Tunneling Projects. MIT, PhD Thesis, Cambridge, p. 589. Sousa, R.L., 2012. Risk Assessment in Tunnels Using Bayesian Networks. In: Sousa, Vargas, Matos Fernandes, Azevedo (Eds.), Innovative Numerical Modeling in Geomechanics, CRC Press/Balkema, Leiden, The Netherlands, pp. 211–244. Sousa, R.L., Einstein, H.H., 2012. Risk analysis during tunnel construction using Bayesian Networks: Porto Metro case study. Tunn. Undergr. Space Technol. 27 (January), 86–100. Sousa, R.L., Ivanova, V., Murrilhy, B., Einstein, H.H., 2012. Decision Aids for Engineered Geothermal Systems: Fracture pattern modeling. Research report, Massachussets Institute of Technology, Cambridge, USA. Stallmann, M., 2005. Verbrüche im Tunnelbau Ursachen und Sanierung. Diplomarbeit.Fachhochschule Hochschule Für Stuttgart Technik Stuttgart University of Applied Sciences, p. 122. Sulsky, D., Schreyer, H.L., 1996. Axisymmetric form of the material point method with applications to upsetting and Taylor impact problems. Comput. Methods Appl. Mech. Eng. 139 (1), 409–429. Sulsky, D., Chen, Z., Schreyer, H.L., 1994. A particle method for history-dependent materials. Comput. Methods Appl. Mech. Eng. 118, 179–196. Sulsky, D., Zhou, S.J., Schreyer, H.L., 1995. Application of a particle-in-cell method to solid mechanics. Comput. Phys. Commun. 87, 236–252. Truesdell, C., Noll, W., 1965. Nonlinear field theories of mechanics. Handbuchder Physic. II/3, Springer. Vlasov, S.N., Makovsky, L.V., Merkin, V.E., 2001. Accidents in Transportation and Subway Tunnels. Construction to Operation. Russian Tunneling Association, Moscow, p. 198. Xi, S.G., 2009. Study on floor heave of soft-rock tunnel mechanism to the large deformation in JiaHe mine by LDEAS. MSc thesis, China University of Mining and Technology (in Chinese). Zabala, F., Alonso, E.E., 2012. The material point method and the analysis of dams and dam failures. In: Sousa, Vargas, Fernandes, Azevedo (Eds.), Innovative Numerical Modelling in Geomechanics, Chapter 8, CRC Press, pp. 171–177. Zhang, H.W., Wang, K.P., 2009. Material point method for dynamic analysis of saturated porous media under external contact/impact of solid bodies. Comput. Methods Appl. Mech. Eng. 198, 1456–1472.