- Email: [email protected]
Investigation of stress wave induced cracking behavior of underground rock mass by the numerical manifold method
Tunnelling and Underground Space Technology 92 (2019) 103032
Contents lists available at ScienceDirect
Tunnelling and Underground Space Technology journal homepage: www.elsevier.com/locate/tust
Investigation of stress wave induced cracking behavior of underground rock mass by the numerical manifold method
T
⁎
L.F. Fana, X.F. Zhoub, Z.J. Wuc, , L.J. Wanga a
College of Architecture and Civil Engineering, Beijing University of Technology, Beijing 100124, China Department of Civil Engineering, Zhejiang University, Yuhangtang Road, Hangzhou 310058, China c School of Civil Engineering, Wuhan University, Wuhan 430072, China b
A R T I C LE I N FO
A B S T R A C T
Keywords: Numerical manifold method Cracking Transient stress field Stress wave propagation Rock mass
The present study investigated the coupling effects of cracking and wave propagation on the underground rock mass by the numerical manifold method (NMM). The traditional NMM was developed for the simulation of dynamic cracking behavior of rock mass. One-dimensional wave propagation was utilized to validate the present code. Subsequently, the transient stress field in the rock mass with an inclined crack was investigated. Finally, two cases of underground caverns with symmetrical or unsymmetrical discontinuities under stress wave were simulated to validate the application potentials. The results show that the transient stress field is not only influenced by the stress concentration around the pre-existing crack, but also influenced by the characteristics of stress wave, e.g. reflection, transmission and diffraction, which are different from the static situation. The results also show that the cracking significantly affects the stress wave duration and amplitude. On the other hand, the stress wave induced cracking depends on the wave propagation distance significantly. Moreover, the NMM can be used to simulate the cracking behavior of underground rock mass under stress wave efficiently.
1. Introduction The discontinuities in the rock mass have significant effects on the wave propagation, which have been widely investigated (Leucci and Giorgi, 2006; Zhao et al., 2008; Fan et al., 2013; Dong et al., 2014). On the other hand, the stress wave may also cause the cracking, which reduces the strength of rock mass. The coupled effects of cracking and stress wave propagation may result in the partial failure or even overall collapse of underground engineering. Therefore, it is of essential significance for the rock engineering to investigate the cracking process of rock mass under stress wave propagation. Cracking was originally observed under the static situation. A large amount of theories and experiments were conducted (Yang and Jing, 2011; Yang et al., 2014a). The studies of the stress field distribution around an inclined crack in the rock mass under static loading showed that the stress field is centrosymmetric with respect to the center of crack under static (or quasi-static) loading such as uniaxial compression or biaxial compression. Moreover, the stress around the crack tips was much larger than that in the middle portion of the crack because of the stress concentration (Hoek and Bieniawski, 1965; Irwin, 1997; Rossmanith, 2014). Besides the stress field, the crack propagation process was further observed based on the different materials, e.g.
⁎
PMMA, gypsum, marble etc. (Wong and Chau, 1998; Sagong and Bobet, 2002; Wong and Einstein, 2009a, 2009b; Lee and Jeon, 2011; Nejati and Ghazvinian, 2014). Seven types of a single crack, including four tensile crack types and three shear crack types, which were characterized by cracking orientation were summarized (Wong and Einstein, 2009a, 2009b). For the cases of more than two pre-existing cracks, the crack spacing considerably affects the crack propagation and connection (Shen et al., 1995). Three main modes including shear mode, the mixed mode and wing tensile mode of two parallel cracks coalescence were classified (Wong and Chau, 1998). Based on the damage mechanics, Wang (1992a) proposed the Wang’s model for the ductile fracture. The model has been applied successfully (Wang, 1991, 1992b). Previous studies showed that the cracking behavior under static loading significantly affected by the length, inclination angle of pre-existing crack and the mechanical properties of the material (Wong and Einstein, 2009a, 2009b). Dynamic loadings were commonly observed in the practical engineering in the forms of seismic wave and explosion wave (Fan et al., 2018; Li et al., 2013; Hu et al., 2018; Xia et al., 2018; Xie et al., 2018). More and more attentions were paid to the cracking behavior under dynamic loadings recently (Biswas et al., 2012; Yang et al., 2014b). The dynamic cracking was experimentally observed by using the Split
Corresponding author. E-mail address: [email protected] (Z.J. Wu).
https://doi.org/10.1016/j.tust.2019.103032 Received 31 December 2018; Received in revised form 17 May 2019; Accepted 6 July 2019 0886-7798/ © 2019 Elsevier Ltd. All rights reserved.
Tunnelling and Underground Space Technology 92 (2019) 103032
L.F. Fan, et al.
where KI and KII are the stress intensity factors (SIFs) of mode I and II, respectively. KIC is the fracture roughness of mode I and θ0 is the crack propagation angle (An, 2010; Wu and Wong, 2012; Wu et al., 2013). The pre-existing crack starts to propagate when
Hopkinson Pressure Bar (SHPB) (Fan et al., 2017) and verified by the numerical simulation (Grégoire et al., 2007). Failure processes of rock specimens under static loading and dynamic loading were compared and it was concluded that more cracks were developed under the dynamic loading compared with the static loading (Zhu and Tang, 2006). However, the dynamic cracking mechanism was much difficultly to be obtained experimentally since the stress field induced by the dynamic loading was transient, so the stress field was influenced not only by the stress concentration but also by the wave transmission and reflection. Moreover, the mechanical properties of rock mass were highly influenced by the loading rate, which further affects the cracking mechanism under dynamic loading. Considering the above difficulties in experiments, numerical methods are under exploring by researchers alternatively (Deng et al., 2014; Li et al., 2018; Zhao and Xia, 2018; Zhao et al., 2018; Yang et al., 2018; Zhao et al., 2019; Zhu et al., 2013). The finite element method (FEM) was usually applied for simulating the crack process(Wu et al., 2019a,b). However, re-mesh was required when using the FEM to simulate the cracking behavior. Extended finite element method (XFEM) and generalized finite element method (GFEM) were also developed based on the FEM to model the discontinuities such as individual crack and intersecting cracks (Strouboulis et al., 2000; Grégoire et al., 2007; Rannou et al., 2009; Zhang and Feng, 2011). Besides the FEM, the numerical manifold method (NMM) proposed by Shi (1992) provided a good way to solve both continuous and discontinuous problems. The method involved a mathematical cover (MC) and a physical cover (PC), which were independent. Therefore, the cracking process could be well investigated using a uniform system. Based on the NMM, Zhang (2010) modeled complex crack propagation involving multiple and branched cracks. Ning (2011) simulated four typical static crack problems and then analyzed the footwall slope in stability. Zhang (2015) further developed the NMM for the simulation of hydraulic fracture in the star-shape. Besides the simulation of mechanical behavior of rock engineering under static loading, the stress wave propagation could also be conveniently modeled, e.g. wave reflection and transmission as well as the attenuation and diffusion when the stress wave propagating through the discontinuities (Fan et al., 2013; Zhou et al., 2017). Wu (2013) simulated the crack initiation and propagation with both the elastic model and the viscoelastic model under dynamic loading. However, the interaction between the stress wave propagation and crack propagation has not been fully explored for the underground engineering. The detailed advantages of NMM was summarized (An, 2010) and a comprehensive review of the NMM can be found (Ma et al., 2010). This paper focuses on both cracking process, stress field variation and wave propagation. The original NMM was developed to simulate the cracking under dynamic loading. Wave propagation through the rock mass was simulated to validate the present code and study the transient stress field, dynamic cracking pattern and the effects of cracking on the wave transmission and reflection. Case studies of underground caverns stability analysis were given to show the coupling effects of cracking and stress wave propagation and further validate the potential application of the NMM on the underground engineering.
(2)
K e = KIC
where Ke represents the equivalent stress intensity factor. The calculation process can be found in the previous research (Wu and Wong, 2012; Wu et al., 2013). The maximum circumferential stress criterion considered that the crack propagates along the angle θ0 which is corresponding to the maximum hoop normal stress and can be gained by
KI sin θ0 + KII (3 cos θ0 − 1) = 0
(3)
Therefore, θ0 is calculated as (An, 2010; Wu and Wong, 2012; Wu et al., 2013)
⎧ 2 arctan (KI − ⎪ θ0 =
⎨ 2 arctan ⎪ ⎪0 ⎩
2) KI2 + 8KII
4KII
2) (KI + KI2 + 8KII
4KII
KII > 0 KII < 0 KII = 0
(4)
2.2. Crack propagation using NMM In this study, each overlapping hexagonal cover having six triangular elements that share the same vertex (called star in the NMM) is termed as a mathematical cover (MC), denoted by Mi (i = 1–5). A MC may be intersected by external boundaries and internal crack into several domains, then each domain is termed as a physical cover (PC), denoted by Pi j (j = 1, 2). The common area of two or more PCs is defined as a manifold element (ME), denoted by e1 and e2 as shown in Fig. 1(a). A simple crack in a physical domain is illustrated in Fig. 1(a). Two different kinds of elements are produced because of the crack named as fully cracked elements (e.g., e1) and partially cracked elements (e.g., e2). The three vertexes of the element (e1) are the stars of three mathematical covers (M1, M2 and M3). M1 and M2 are completely cut by the crack into two independent PCs, P11 and P12 (Fig. 1(b)), P21 and P22 (Fig. 1(c)), respectively, while M3 is only partially cut by the crack and forms one PC of P3 (Fig. 1(d)). Element e1 is divided into two parts: e11 (associated with P11, P21 and P3) and e12 (associated with P12 , P22 and P3). The displacement jump across the crack is expressed as
[[uh (x )]] = ∑i φi (x )·ui1 (x ) − ∑i φi (x )·ui2 (x ) = φ1 (e1 )(u11 (e1) − u12 (e1)) + φ2 (e1 )(u21 (e1) − u22 (e1))
(5)
where [[·]] is the jump function. φi is the weight function corresponding to the physical cover. ui1 and ui2 are the cover functions associated with P11 and P12 or P21 and P22 , respectively (An, 2010; Wu and Wong, 2012; Wu et al., 2013, 2017). When it turns to element e2, the three vertexes of this element are the stars of mathematical covers M3, M4 and M5. All of these mathematical covers are partially cut by a crack. Thus, element e2 is associated with P3 (Fig. 1(d)), P4 (Fig. 1(e)) and P5 (Fig. 1(f)), so Eq. (5) is not applicable anymore. Enrichment methods are proposed to solve this problem as
2. Basic concepts 2.1. Crack propagation criteria
3
3
Extensive studies have been conducted with the cracking initiation criterion (Biswas et al., 2012; Yang et al., 2014b) and crack propagation criterion (Grégoire et al., 2007; Zhang and Feng, 2011). In this study, the maximum circumferential stress criterion is adopted. The criterion of crack instability is
uj (x ) = ϕ·cj
1 θ cos 0 [KI (1 + cos θ0) − 3KII sin θ0] = KIC 2 2
where cj represents the array of additional unknowns.ϕ is the matrix of polynomial bases as
uh (x ) =
∑ φi (x )·ui (x ) + ∑ φi (x )·uj (x ) i=1
i=1
(6)
The uj (x ) is obtained by
(1) 2
(7)
Tunnelling and Underground Space Technology 92 (2019) 103032
L.F. Fan, et al.
Fig. 1. Basic concepts and crack modeling of NMM: a basic concepts and an internal crack, b completely cut of M1, c completely cut of M2, d partially cut of M3, e partially cut of M4, f partially cut of M5.
(An, 2010; Wu and Wong, 2012; Wu et al., 2013).
ϕ1 0 ϕ2 0 ϕ3 0 ϕ4 0 ⎤ ϕ=⎡ ⎢ 0 ϕ1 0 ϕ2 0 ϕ3 0 ϕ4 ⎥ ⎦ ⎣
(8)
where ϕ1, ϕ2 , ϕ3 and ϕ4 can be obtained by
θ [ϕ1, ϕ2 , ϕ3, ϕ4] = ⎡ r sin , 2 ⎣
θ r cos , 2
θ r sin θ sin , 2
2.3. Loop update
θ r sin θ cos ⎤ 2⎦
The loops should be updated once cracks propagated for the contact detection and crack presentation (Wu and Wong, 2012). The modification to the loops includes four types (more details in (Sagong and Bobet, 2002)). Once a new crack exists, the ME loop and PC loop would
(9) where (r, θ) are polar coordinates in the crack tip coordinate system 3
Tunnelling and Underground Space Technology 92 (2019) 103032
L.F. Fan, et al.
both change. Simple examples of ME loop update and PC loop update are given in Fig. 2 for explaining the update procedure. There are two types of loop update in the given examples. New crack would divide the original loop into two parts. Therefore the new loop is added (e.g. e2 to e5 and P1 to P3). On the other hand, the new crack may combine the original two loops into a loop (e.g. P2 to P4) (more details in (Sagong and Bobet, 2002; Wu and Wong, 2012)). 3. Stress wave propagation through intact rock bar The stress wave propagation through an intact rock bar was carried out to calibrate the present code. Fig. 3 shows a rock bar with a length of x = 2.00 m and a width of y = 0.10 m. The incident wave is achieved by dynamic loading at location A of xA = 0.0001 m as shown in Fig. 3. 10 evenly distributed measurement points at each location are arranged to avoid the measurement error. Three measurement locations are placed at xB = 0.50 m, xC = 1.00 m and xD = 1.50 m, respectively. The rock density is ρ = 2640 kg/m3, the Young’s modulus is E = 60 GPa and the Poisson ratio is ν = 0.3. In order to reduce lateral effects of stress wave and to satisfy the one-dimensional wave propagation theory, the length of the numerical model is ten times larger than the width. An incident stress wave of half sinusoidal waveform is applied at the left boundary of the rock bar as
0 ⩽ t ⩽ 1/2f0 σ sin(2πf0 t ) , when ⎧ σ =⎧ 0 ⎨ others ⎨ ⎩0 ⎩
(10)
where σ0 and f0 denote the amplitude and frequency of the incident wave, respectively. In the present study, σ0 = 10.0 MPa and f0 = 2000 Hz. Fig. 4 shows the simulation results. Since the stress wave transmits from the left of the bar to the right in this model, there is a time delay between the incident wave and each transmitted wave. The interval of two waves could be calculated by t = l/C, where l is the distance between the two corresponding measurement locations and C is the wave velocity obtained by C = E / ρ . It is also seen from Fig. 4 that the amplitudes of transmitted wave are 99.82%, 99.02% and 98.83% of the incident wave. The results show that the calculation induced error is within 1.20%, which is acceptable in the rock dynamic engineering. 4. Stress wave propagation through incline cracked rock bar 4.1. Dynamic stress field distribution The existence of crack leads to stress concentration around the tips of crack which results in the cracking of rock mass further influences the stress wave propagation. To investigate the effects of cracking on dynamic stress field, stress wave propagation through incline cracked rock bar is simulated as shown in the Fig. 5. To show the advantages of NMM in the simulation of stress wave propagation, the incline angles of the crack were selected at random. In the present study, angles of 13° and 31° were illustrated. The tensile strength and compressive strength of the rock bar are 30.0 MPa and 200.0 MPa, respectively. The cohesion is c = 15 MPa and the friction angle is φ = 30°. The parameters of the numerical model are selected according to the previous researches (Fan et al., 2013; Wu and Wong, 2012; Wu et al., 2013). Three rows of measurements are arranged with 10 points at location E, 25 points at location F and 25 points at location G, respectively. The measurement points located at locations F and G are parallel to the crack while the measurement points at location E are vertical. All points are numbered from the bottom to top subsequently (1 to 10 at location E and 1 to 25 at locations F and G, respectively). The horizontal locations of the crack tips are x = 0.82 m and x = 1.08 m and the angle of the crack is 13°.
Fig. 2. Update of ME loop and PC loop when cracking: a update of ME loop, b update of PC loop.
4
Tunnelling and Underground Space Technology 92 (2019) 103032
L.F. Fan, et al.
B
0.50 m
D
C
0.50 m
0.50 m
0.50 m
2.00 m
0.10 m
A
Fig. 3. Scheme of an intact rock model.
Fig. 4. Illustration of stress histories of measurement points in rock mass under dynamic non-equilibrium stress field.
The stress histories for the wave propagation through pre-existing oblique cracked rock bar are monitored and shown in Fig. 6. It is observed from Fig. 6(a) that a slight rise was observed at about 0.45 ms in Fig. 6(a). Since the main purpose of the present research is to study the stress field around the crack tips, the stress histories at points E-01 to E10 were studied. It is seen from Fig. 6(a) that the stress histories of ten measurement points E-01 to E-10 of x = 0.1 m show a good agreement with the amplitude of 1.0 MPa, which reveals that the lateral effects can be ignored in the present study. Moreover, it is also shown that the NMM can be utilized to simulate the one-dimensional wave propagation through rock mass with a convergent result. Fig. 6(b) and (c) show the stress distribution histories in front of and behind of the oblique crack, respectively. According to the previous researches under static loading (Li, 2010), the stress distribution around the crack is constant and in the form of centrosymmetric to the center of the crack. While for the stress wave propagation, the stress distribution varies as time increases. It is seen from Fig. 6(b) that the peak stresses around the tip in front of the crack (e.g. F-25) are larger than the amplitude of the incident wave (1.0 MPa). Therefore, stress concentration around crack tips also can be observed under dynamic loading. However, different phenomenon from static loading, it is also seen from Fig. 6(b) that the stress history around the bottom tip in front of the crack (e.g. F-01) are with different waveforms, which is induced by the reflected stress waves. Moreover, obvious tensile stress field (e.g. F-01, 02, 03 and so on) can be observed around the bottom tip in front of the crack.
Fig. 6. Stress distribution inside the rock model without cracking under dynamic non-equilibrium stress field.
The stress history behind the crack shows similar waveforms as the incident wave as shown in Fig. 6(c). It is seen from Fig. 6(c) that distinct stress concentration also can be observed around both the bottom and top tips of the crack (e.g. G-01 and G-25). The concentration rate around the bottom tip behind the crack is larger than that around the top tip behind the crack, which is in accordance with the previous researches under static loading (Li, 2010). It should also be noted that unlike the stress wave propagation through crack normally, the times associated with the peak stresses are different, which are induced by the oblique geometrical characters of the pre-existing crack results in the different wave propagation time. To investigate the stress field variation under dynamic loading, the
Fig. 5. Scheme of a rock model with a long crack and oblique measurement points. 5
Tunnelling and Underground Space Technology 92 (2019) 103032
L.F. Fan, et al.
transmitted waveforms recorded at location H convergent before they were influenced by the cracking. It is seen from Fig. 10 that tensile stress was recorded, which was reflected by the crack. The phenomenon is in accordance with the previous researches. However, it is interesting that the amplitude of the tensile stress decreases for the measurement points from H-01 to H-10, which is caused by the cracking induced dynamic energy loss. Such results cannot be simulated in the previous researches using the traditional manifold methods (Fan et al., 2013). It is also seen from Fig. 10 that the crack propagation has significant effects on the stress field at location I. It is seen that the stress waves at ten measurement points I-01 to I-10 convergent to each other before cracking, which is at about 0.52 ms for the present study. The amplitude of stress waves at measurement points from I to 01 to I-10 decreases after cracking. The phenomenon is similar as the effects of cracking on the wave propagation, which also could not be simulated by the traditional NMM (Fan et al., 2013). 5. Case study
Fig. 7. Inhomogeneous stress field around the crack under dynamic non-equilibrium stress field.
To validate the engineering application potentials of the developed NMM, two cases of underground cavern explosion with symmetrical and unsymmetrical discontinuities under dynamic stress wave propagation were investigated.
stress distribution at four typical times are investigated and shown in Fig. 7. Four typical times are selected as the times of the stresses around top tip and bottom tip in front of the crack, and top tip and bottom tip behind the crack reach their peaks as shown in the Fig. 6(b) and (c), respectively. It is seen from Fig. 7 that the stress increases from the bottom tip to top tip in front of the crack, while decreases from the bottom tip to top tip behind the crack in general. The variation of stress distribution is similar as that under static loading at times at peaks of F25, G-01 and G-25. The inflection point of stress distribution is not located at the tips of the crack, while very close to the tips, such as measurement point 4 at times at peaks of F-01, F-25, G-01 and G-25, and measurement point 22 at times at peaks of F-25, G-01 and G-25, respectively. In addition, it is also seen from Fig. 7 that the stress distribution is not centrosymmetric, which can be observed at time at peak of F-01. The phenomenon is caused by the reflected stress wave from the inclined crack, which is agreed with the results from Fig. 6.
5.1. Symmetrical cracking Fig. 11 shows the cross-section of a test site. The tunnel is circular shaped with a radius of 5.0 m and located at 50 m under the ground. To decrease the boundary effects, the cavern is in the center of a square with a length of 100 m. Displacements at three boundaries are fixed as shown in the Fig. 11. The properties of rock are the same as mentioned in the Section 4. A crack with a length of 27.94 m is located at the ramp above the cavern. The distances of two ends of the crack to the center point of the cavern are equal, namely, l1 = l2. Therefore, the stress wave arrives at both endpoints of the crack simultaneously and the stress value is equal. In the present model, the measurement points of row A and row C are horizontally located while those of row B are vertically located. There are 14 measurement points in each row and the space between two neighboring points is 2.0 m. The measurement points in each row are numbered as 1–14 from the center to the border. A dynamic loading is applied on the internal boundary of the tunnel to simulate a normally incident wave. A stress wave with multiple peaks of different amplitudes is applied with a waveform of
4.2. Cracking mode The failure criteria as shown in Eqs. (1) and (2) is introduced in the present NMM code. Fig. 8 shows the numerical model to investigate the cracking mode and the stress distribution variation under dynamic loading. An initial inclined crack starts at x = 0.90 m and ends at x = 1.00 m with an incline angle of 31° is assumed. The measurement points located at xH = 0.3 m and xI = 1.3 m are set to analyze the effects of cracking on the wave propagation, respectively. Each location with 10 equal-spaced measurement points as in the previous sections. Fig. 9(a) and (b) show a visual result of dynamic crack propagation. The cracking process can be divided into two stages. The first stage is that the crack tip closer to the dynamic loading starts to propagate when the stress reaches its critical value while the farther crack tip maintains the original state as shown in Fig. 9(a). The second stage is that the other crack tip propagates when the stress wave arrives the right tip as shown in Fig. 9(b). Fig. 10 gives the measuring data at locations H and I. The
⎧ 0 ⩽ t ⩽ 1/(2f1 ) ⎧ σ1 sin(2πf1 × t ) σI = σ2 sin(2πf2 × t ) , when 1/(2f1 ) + Δ T⩽ t ⩽ 1/(2f2 ) + 1/(2f1 ) + ΔT ⎨ ⎨ ⎩0 ⎩ others (11) where ΔT denotes the time period between two pulses of the wave. σ1 and σ2 are the amplitudes of the wave with multiple peaks, respectively. f1 and f2 are the corresponding frequencies. The numerical result of two symmetrically wing cracks with equal distances from the crack tips could be seen in Fig. 12. For a symmetrical crack case, the two crack tips start to propagate at the same time and the crack paths are almost symmetrical. The new cracks both through
Fig. 8. Scheme of a rock model with a short crack and vertical measurement points. 6
Tunnelling and Underground Space Technology 92 (2019) 103032
L.F. Fan, et al.
Fig. 9. Crack propagation under dynamic non-equilibrium stress field: a crack propagation at left tip of the crack, b crack propagation at right tip of the crack.
Fig. 11. Scheme of an underground cavern with an equidistant crack.
incident wave, respectively. Different from the one-dimensional wave propagation through rock mass, which is caused by the micro-defect induced attenuation, the attenuation in the present study is mainly caused by the wave front circular expansion, which reduces the dynamic energy density further results in the wave attenuation. Moreover, it is also can be observed from Fig. 13 that tensile stress is introduced between two transmitted pulses, which is caused by the compression stress wave applied on the internal surface of the circle carven. Fig. 10. Stress histories of measurement points during the dynamic crack process.
5.2. Unsymmetrical cracking Just as the symmetrical cracking case, an unsymmetrical crack with a length of 30.04 m is established as shown in Fig. 14. The distances from two crack tips to the center point of the cavern are unequal (e.g. l3 < l4.). Therefore, the stress wave reaches the two tips in sequence. The boundary conditions, material properties and dynamic loading are the same as mentioned in the Section 5.1. Fig. 15 shows the numerical result of crack propagation. It is seen from Fig. 15 that the cracking initiates unsymmetrically, which is
the measurement points of 8th and 9th of row A and row C which influence the wave propagation, respectively. Fig. 13 shows the stress wave histories. The transmitted wave is the stress history recorded at the measurement point A-1 (the 1st measurement point of row A). It can be observed from Fig. 13 that both first pulse and second pulse attenuate. The amplitudes of the first pulse and second pulse are 58.54% and 50.49% comparing to amplitudes of the 7
Tunnelling and Underground Space Technology 92 (2019) 103032
L.F. Fan, et al.
Fig. 14. Scheme of an underground cavern with a non-equidistant crack.
Fig. 12. Simulation result of equidistant crack propagation under dynamic nonequilibrium stress field.
Fig. 13. Waveforms of the incident wave and a typical transmitted wave for underground cavern explosion.
significantly different from the case of cracking under static loading. The crack tip closer to the wave source propagates while the other crack tip keeps constant. The phenomenon can be explained by the dynamic loading attenuation and the interaction between the dynamic loading and cracking. Since the wave front circle enlarges as the stress wave propagates from the wave source, the stress amplitude decreases as the distance from the wave source increases. Therefore, the farther tip of crack may not satisfy the cracking criteria even the closer tip of crack propagates because of the wave attenuation. Moreover, the cracking of the closer tip absorbs the dynamic energy of stress wave, which further increases the wave attenuation. A stress wave as shown in Fig. 13 is applied. Fig. 16 shows the amplitudes of the recorded stress histories at the measurement locations D, E and F. It is seen from Fig. 16 that the new crack passes through the measurement points 8th and 9th of row D according to the simulation result. The second pulse at row D attenuated after propagation of the new crack. Due to the unilateral cracking, the measurement points at row E
Fig. 15. Simulation result of non-equidistant crack propagation under dynamic non-equilibrium stress field.
and row F are not influenced by the pre-existing crack and the new crack. They are agreed well with each other as shown in Fig. 16. In contrast, there is a distinct decrease of stress at row D when the stress wave propagates through the new crack. 6. Conclusions The present study investigates the interaction between dynamic stress wave propagation and the cracking behavior using the NMM. The stress field around the crack tips, cracking modes of pre-existing crack and the stress field after crack propagation are analyzed. The NMM has its advantage in investigating the cracking behavior and dynamic stress field variation inside the rock mass under stress wave. It can be concluded that a transient stress field would be informed around the crack tips during the stress wave propagation. The dynamic 8
Tunnelling and Underground Space Technology 92 (2019) 103032
L.F. Fan, et al.
adjacent tunnel explosion. Tunn. Undergr. Sp. Tech. 35, 227–234. Li, X., Li, X.F., Zhang, Q.B., et al., 2018. A numerical study of spalling and related rockburst under dynamic disturbance using a particle-based numerical manifold method (PNMM). Tunn. Undergr. Sp. Tech. 81, 438–449. Ma, G.W., An, X.M., He, L., 2010. The numerical manifold method: A review. Int. J. Comput. Methods 7 (1), 1–32. Ning, Y., An, X., Ma, G., 2011. Footwall slope stability analysis with the numerical manifold method. Int. J. Rock. Mech. Min. 48 (6), 964–975. Nejati, H.R., Ghazvinian, A., 2014. Brittleness effect on rock fatigue damage evolution. Rock Mech. Rock Eng. 47 (5), 1839–1848. Rannou, J., Gravouil, A., Baïetto-Dubourg, M.C., 2009. A local multigrid x-fem strategy for 3-d crack propagation. Int. J. Numer. Meth. Eng. 77 (4), 581–600. Rossmanith, H.P. (Ed.), 2014. RockOCK FRACT. MECH., vol. 275 Springer. Shi, G., 1992. Modeling rock joints and blocks by manifold method. In: Proceedings of the 33rd US Rock Mechanics Symposium. New Mexico, pp. 639–648. Shen, B., Stephansson, O., Einstein, H.H., et al., 1995. Coalescence of fractures under shear stresses in experiments. J. Geophy. Res. 100 (B4), 5975–5990. Strouboulis, T., Babuška, I., Copps, K., 2000. The design and analysis of the generalized finite element method. Comput. Method. Appl. M. 181 (1), 43–69. Sagong, M., Bobet, A., 2002. Coalescence of multiple flaws in a rock-model material in uniaxial compression. Int. J. Rock. Mech. Min. 39 (2), 229–241. Wang, T.J., 1991. A continuum damage model for ductile fracture of weld heat affected zone. Eng. Fracture Mech. 40 (6), 1075–1082. Wang, T.J., 1992a. Unified CDM model and local criterion for ductile fracture-I. Unified CDM model for ductile fracture. Eng. Fracture Mech. 42 (1), 177–183. Wang, T.J., 1992b. Unified CDM model and local criterion for ductile fracture-II. Ductile fracture local criterion. Eng. Fracture Mech. 42 (1), 185–192. Wong, R.H.C., Chau, K.T., 1998. Crack coalescence in a rock-like material containing two cracks. Int. J. Rock. Mech. Min. 35 (2), 147–164. Wong, L.N.Y., Einstein, H.H., 2009a. Crack coalescence in molded gypsum and Carrara marble: Part 1. Macroscopic observations and interpretation. Rock Mech. Rock Eng. 42 (3), 475–511. Wong, L.N.Y., Einstein, H.H., 2009b. Systematic evaluation of cracking behavior in specimens containing single flaws under uniaxial compression. Int. J. Rock Mech. Min. 46 (2), 239–249. Wu, Z.J., Fan, L.F., Liu, Q.S., Ma, G.W., 2017. Micro-mechanical modeling of the macromechanical response and fracture behavior of rock using the numerical manifold method. Eng. Geol. 225 (20), 49–60. Wu, Z., Wong, L.N.Y., 2012. Frictional crack initiation and propagation analysis using the numerical manifold method. Comput. Geotech. 39 (1), 38–53. Wu, Z., Wong, L.N.Y., Fan, L., 2013. Dynamic study on fracture problems in viscoelastic sedimentary rocks using the numerical manifold method. Rock Mech. Rock Eng. 46 (6), 1415–1427. Yang, S.Q., Jing, H.W., 2011. Strength failure and crack coalescence behavior of brittle sandstone samples containing a single fissure under uniaxial compression. Int. J. Fract. 168 (2), 227–250. Yang, S.Q., Huang, Y.H., Jing, H.W., et al., 2014a. Discrete element modeling on fracture coalescence behavior of red sandstone containing two unparallel fissures under uniaxial compression. Eng. Geol. 178, 28–48. Yang, L., Yang, R., Qu, G., et al., 2014b. Caustic study on blast-induced wing crack behaviors in dynamic–static superimposed stress field. Int. J. Min. Sci. Tech. 24 (4), 417–423. Yang, Y.T., Tang, X.H., Zheng, H., et al., 2018. Hydraulic fracturing modeling using the enriched numerical manifold method. Appl. Math. Model. 53, 462–486. Wu, Z.J., Zhang, P.L., Fan, L.F., 2019a. Numerical study of the effect of confining pressure on the rock breakage efficiency and fragment size distribution of a TBM cutter using a coupled FEM-DEM method. Tunn. Undergr. Sp. Tech. 88, 260–275. Wu, Z.J., Zhou, Y., Fan, L.F., 2019b. A fracture aperture dependent thermal-cohesive coupled model for modelling thermal conduction in fractured rock mass. Comput. Geotech. 114, 103108. Xia, X., Li, H.B., Lu, Y.Q., 2018. A case study on the cavity effect of a water tunnel on the ground vibrations induced by excavating blasts. Tunn. Undergr. Sp. Tech. 71, 292–297. Xie, L.X., Lu, W.B., Zhang, Q.B., et al., 2018. Analysis of damage mechanisms and optimization of cut blasting design under high in-situ stresses. Tunn. Undergr. Sp. Tech. 66, 19–33. Zhu, W., Tang, C., 2006. Numerical simulation of Brazilian disk rock failure under static and dynamic loading. Int. J. Rock Mech. Min. 43 (2), 236–252. Zhao, X., Zhao, J., Cai, J., et al., 2008. UDEC modelling on wave propagation across fractured rock masses. Comput. Geotech. 35 (1), 97–104. Zhang, H., Li, L., An, X., et al., 2010. Numerical analysis of 2-d crack propagation problems using the numerical manifold method. Eng. Anal. Bound. Elem. 34 (1), 41–50. Zhang, Y., Feng, X., 2011. Extended finite element simulation of crack propagation in fractured rock masses. Mater. Res. Innov. 15 (s1), s594–s596. Zhang, G., Xu, L., Haifeng, L., 2015. Simulation of hydraulic fracture utilizing numerical manifold method. Sci China Tech. Sci. 58 (9), 1542–1557. Zhao, G.F., Xia, K., 2018. A study of mode-I self-similar dynamic crack propagation using a lattice spring model. Comput. Geotech. 96, 215–225. Zhao, G.F., Hu, X., Li, Q., et al., 2018. On the four-dimensional lattice spring model for geomechanics. J. Rock Mech. Geotech. Eng. 10 (4), 661–668. Zhao, G.F., Yin, Q., Russell, A.R., et al., 2019. On the linear elastic responses of the 2D bonded discrete element model. Int. J. Number. Anal. Met. 43 (1), 166–182. Zhu, J.B., Deng, X.F., Zhao, X.B., Zhao, J., 2013. A numerical study on wave transmission across multiple intersecting joint sets in rock masses with UDEC. Rock Mech. Rock Eng. 46 (6), 1429–1442. Zhou, X.F., Fan, L.F., Wu, Z.J., 2017. Effects of microfracture on wave propagation through rock mass. Int. J. Geomech. 17 (9), 04017072.
Fig. 16. Peak stress variation before and after crack propagation of non-equidistant crack under dynamic non-equilibrium stress field.
stress field is different from that under static loading. The transient stress field is induced by both stress concentration and the reflection and transmission of stress wave. The characteristics of transient stress field result in special cracking modes. The crack tip closer to the loading points propagates earlier than the farther one during the stress wave propagation. In addition, the cracking process further influences the stress field distribution. Acknowledgements The research is supported by the National Natural Science Foundation of China (NSFC) (Nos. 51778021, 11572282 and 41772309) and the Science Fund for Creative Research Groups of the National Natural Science Foundation of China (No. 51421005). References An, X.M., 2010. Extended numerical manifold method for engineering failure analysis. Ph.D thesis. Technological University. Singapore, Nanyang. Biswas, N., Ding, J., Balla, V.K., et al., 2012. Deformation and fracture behavior of laser processed dense and porous ti6al4v alloy under static and dynamic loading. Mater. Sci. Eng. A. 549, 213–221. Deng, X.F., Zhu, J.B., Chen, S.G., Zhao, Z.Y., Zhou, Y.X., Zhao, J., 2014. Numerical study on tunnel damage subject to blast-induced shock wave in jointed rock masses. Tunn. Undergr. Sp. Tech. 43, 88–100. Dong, Q., Li, X., Zhao, H., 2014. Experimental research on ultrasonic p-wave velocity variation of fractured rock mass under different stress paths. Adv. Mater. Res. 1065–1069, 35–39. Fan, L.F., Wang, L.J., Wu, Z.J., 2018. Wave transmission across linearly jointed complex rock masses. Int. J. Rock Mech. Min. 112, 193–200. Fan, L.F., Wu, Z.J., Wan, Z., Gao, J.W., 2017. Experimental investigation of thermal effects on dynamic behavior of granite. Appl. Therm. Eng. 125, 94–103. Fan, L., Yi, X., Ma, G., 2013. Numerical manifold method (NMM) simulation of stress wave propagation through fractured rock mass. Int. J. Appl. Mech. 05 (02), 1350022. Grégoire, D., Maigre, H., Rethore, J., et al., 2007. Dynamic crack propagation under mixed-mode loading-comparison between experiments and X-FEM simulations. Int. J. Solids. Struct. 44 (20), 6517–6534. Hoek, E., Bieniawski, Z.T., 1965. Brittle fracture propagation in rock under compression. Int. J. Fract. Mech. 1 (3), 137–155. Hu, X.D., Zhao, G.F., Deng, X.F., et al., 2018. Application of the four-dimensional lattice spring model for blasting wave propagation around the underground rock cavern. Tunn. Undergr. Sp. Tech. 80, 135–147. Irwin, G.R., 1997. Analysis of stresses and strains near the end of a crack traversing a plate. Spie Milestone Series MS 137 (167–170), 16. Leucci, G., Giorgi, L.D., 2006. Experimental studies on the effects of fracture on the p and s wave velocity propagation in sedimentary rock (“calcarenite del salento”). Eng. Geol. 84 (3), 130–142. Li, S.Y., 2010. Introduction to Rock Fracture Mechanics. China University of Science and Technology. Lee, H., Jeon, S., 2011. An experimental and numerical study of fracture coalescence in pre-cracked specimens under uniaxial compression. Int. J. Solid. Stut. 48 (6), 979–999. Li, J.C., Li, H.B., Ma, G.W., et al., 2013. Assessment of underground tunnel stability to
9
Contents lists available at ScienceDirect
Tunnelling and Underground Space Technology journal homepage: www.elsevier.com/locate/tust
Investigation of stress wave induced cracking behavior of underground rock mass by the numerical manifold method
T
⁎
L.F. Fana, X.F. Zhoub, Z.J. Wuc, , L.J. Wanga a
College of Architecture and Civil Engineering, Beijing University of Technology, Beijing 100124, China Department of Civil Engineering, Zhejiang University, Yuhangtang Road, Hangzhou 310058, China c School of Civil Engineering, Wuhan University, Wuhan 430072, China b
A R T I C LE I N FO
A B S T R A C T
Keywords: Numerical manifold method Cracking Transient stress field Stress wave propagation Rock mass
The present study investigated the coupling effects of cracking and wave propagation on the underground rock mass by the numerical manifold method (NMM). The traditional NMM was developed for the simulation of dynamic cracking behavior of rock mass. One-dimensional wave propagation was utilized to validate the present code. Subsequently, the transient stress field in the rock mass with an inclined crack was investigated. Finally, two cases of underground caverns with symmetrical or unsymmetrical discontinuities under stress wave were simulated to validate the application potentials. The results show that the transient stress field is not only influenced by the stress concentration around the pre-existing crack, but also influenced by the characteristics of stress wave, e.g. reflection, transmission and diffraction, which are different from the static situation. The results also show that the cracking significantly affects the stress wave duration and amplitude. On the other hand, the stress wave induced cracking depends on the wave propagation distance significantly. Moreover, the NMM can be used to simulate the cracking behavior of underground rock mass under stress wave efficiently.
1. Introduction The discontinuities in the rock mass have significant effects on the wave propagation, which have been widely investigated (Leucci and Giorgi, 2006; Zhao et al., 2008; Fan et al., 2013; Dong et al., 2014). On the other hand, the stress wave may also cause the cracking, which reduces the strength of rock mass. The coupled effects of cracking and stress wave propagation may result in the partial failure or even overall collapse of underground engineering. Therefore, it is of essential significance for the rock engineering to investigate the cracking process of rock mass under stress wave propagation. Cracking was originally observed under the static situation. A large amount of theories and experiments were conducted (Yang and Jing, 2011; Yang et al., 2014a). The studies of the stress field distribution around an inclined crack in the rock mass under static loading showed that the stress field is centrosymmetric with respect to the center of crack under static (or quasi-static) loading such as uniaxial compression or biaxial compression. Moreover, the stress around the crack tips was much larger than that in the middle portion of the crack because of the stress concentration (Hoek and Bieniawski, 1965; Irwin, 1997; Rossmanith, 2014). Besides the stress field, the crack propagation process was further observed based on the different materials, e.g.
⁎
PMMA, gypsum, marble etc. (Wong and Chau, 1998; Sagong and Bobet, 2002; Wong and Einstein, 2009a, 2009b; Lee and Jeon, 2011; Nejati and Ghazvinian, 2014). Seven types of a single crack, including four tensile crack types and three shear crack types, which were characterized by cracking orientation were summarized (Wong and Einstein, 2009a, 2009b). For the cases of more than two pre-existing cracks, the crack spacing considerably affects the crack propagation and connection (Shen et al., 1995). Three main modes including shear mode, the mixed mode and wing tensile mode of two parallel cracks coalescence were classified (Wong and Chau, 1998). Based on the damage mechanics, Wang (1992a) proposed the Wang’s model for the ductile fracture. The model has been applied successfully (Wang, 1991, 1992b). Previous studies showed that the cracking behavior under static loading significantly affected by the length, inclination angle of pre-existing crack and the mechanical properties of the material (Wong and Einstein, 2009a, 2009b). Dynamic loadings were commonly observed in the practical engineering in the forms of seismic wave and explosion wave (Fan et al., 2018; Li et al., 2013; Hu et al., 2018; Xia et al., 2018; Xie et al., 2018). More and more attentions were paid to the cracking behavior under dynamic loadings recently (Biswas et al., 2012; Yang et al., 2014b). The dynamic cracking was experimentally observed by using the Split
Corresponding author. E-mail address: [email protected] (Z.J. Wu).
https://doi.org/10.1016/j.tust.2019.103032 Received 31 December 2018; Received in revised form 17 May 2019; Accepted 6 July 2019 0886-7798/ © 2019 Elsevier Ltd. All rights reserved.
Tunnelling and Underground Space Technology 92 (2019) 103032
L.F. Fan, et al.
where KI and KII are the stress intensity factors (SIFs) of mode I and II, respectively. KIC is the fracture roughness of mode I and θ0 is the crack propagation angle (An, 2010; Wu and Wong, 2012; Wu et al., 2013). The pre-existing crack starts to propagate when
Hopkinson Pressure Bar (SHPB) (Fan et al., 2017) and verified by the numerical simulation (Grégoire et al., 2007). Failure processes of rock specimens under static loading and dynamic loading were compared and it was concluded that more cracks were developed under the dynamic loading compared with the static loading (Zhu and Tang, 2006). However, the dynamic cracking mechanism was much difficultly to be obtained experimentally since the stress field induced by the dynamic loading was transient, so the stress field was influenced not only by the stress concentration but also by the wave transmission and reflection. Moreover, the mechanical properties of rock mass were highly influenced by the loading rate, which further affects the cracking mechanism under dynamic loading. Considering the above difficulties in experiments, numerical methods are under exploring by researchers alternatively (Deng et al., 2014; Li et al., 2018; Zhao and Xia, 2018; Zhao et al., 2018; Yang et al., 2018; Zhao et al., 2019; Zhu et al., 2013). The finite element method (FEM) was usually applied for simulating the crack process(Wu et al., 2019a,b). However, re-mesh was required when using the FEM to simulate the cracking behavior. Extended finite element method (XFEM) and generalized finite element method (GFEM) were also developed based on the FEM to model the discontinuities such as individual crack and intersecting cracks (Strouboulis et al., 2000; Grégoire et al., 2007; Rannou et al., 2009; Zhang and Feng, 2011). Besides the FEM, the numerical manifold method (NMM) proposed by Shi (1992) provided a good way to solve both continuous and discontinuous problems. The method involved a mathematical cover (MC) and a physical cover (PC), which were independent. Therefore, the cracking process could be well investigated using a uniform system. Based on the NMM, Zhang (2010) modeled complex crack propagation involving multiple and branched cracks. Ning (2011) simulated four typical static crack problems and then analyzed the footwall slope in stability. Zhang (2015) further developed the NMM for the simulation of hydraulic fracture in the star-shape. Besides the simulation of mechanical behavior of rock engineering under static loading, the stress wave propagation could also be conveniently modeled, e.g. wave reflection and transmission as well as the attenuation and diffusion when the stress wave propagating through the discontinuities (Fan et al., 2013; Zhou et al., 2017). Wu (2013) simulated the crack initiation and propagation with both the elastic model and the viscoelastic model under dynamic loading. However, the interaction between the stress wave propagation and crack propagation has not been fully explored for the underground engineering. The detailed advantages of NMM was summarized (An, 2010) and a comprehensive review of the NMM can be found (Ma et al., 2010). This paper focuses on both cracking process, stress field variation and wave propagation. The original NMM was developed to simulate the cracking under dynamic loading. Wave propagation through the rock mass was simulated to validate the present code and study the transient stress field, dynamic cracking pattern and the effects of cracking on the wave transmission and reflection. Case studies of underground caverns stability analysis were given to show the coupling effects of cracking and stress wave propagation and further validate the potential application of the NMM on the underground engineering.
(2)
K e = KIC
where Ke represents the equivalent stress intensity factor. The calculation process can be found in the previous research (Wu and Wong, 2012; Wu et al., 2013). The maximum circumferential stress criterion considered that the crack propagates along the angle θ0 which is corresponding to the maximum hoop normal stress and can be gained by
KI sin θ0 + KII (3 cos θ0 − 1) = 0
(3)
Therefore, θ0 is calculated as (An, 2010; Wu and Wong, 2012; Wu et al., 2013)
⎧ 2 arctan (KI − ⎪ θ0 =
⎨ 2 arctan ⎪ ⎪0 ⎩
2) KI2 + 8KII
4KII
2) (KI + KI2 + 8KII
4KII
KII > 0 KII < 0 KII = 0
(4)
2.2. Crack propagation using NMM In this study, each overlapping hexagonal cover having six triangular elements that share the same vertex (called star in the NMM) is termed as a mathematical cover (MC), denoted by Mi (i = 1–5). A MC may be intersected by external boundaries and internal crack into several domains, then each domain is termed as a physical cover (PC), denoted by Pi j (j = 1, 2). The common area of two or more PCs is defined as a manifold element (ME), denoted by e1 and e2 as shown in Fig. 1(a). A simple crack in a physical domain is illustrated in Fig. 1(a). Two different kinds of elements are produced because of the crack named as fully cracked elements (e.g., e1) and partially cracked elements (e.g., e2). The three vertexes of the element (e1) are the stars of three mathematical covers (M1, M2 and M3). M1 and M2 are completely cut by the crack into two independent PCs, P11 and P12 (Fig. 1(b)), P21 and P22 (Fig. 1(c)), respectively, while M3 is only partially cut by the crack and forms one PC of P3 (Fig. 1(d)). Element e1 is divided into two parts: e11 (associated with P11, P21 and P3) and e12 (associated with P12 , P22 and P3). The displacement jump across the crack is expressed as
[[uh (x )]] = ∑i φi (x )·ui1 (x ) − ∑i φi (x )·ui2 (x ) = φ1 (e1 )(u11 (e1) − u12 (e1)) + φ2 (e1 )(u21 (e1) − u22 (e1))
(5)
where [[·]] is the jump function. φi is the weight function corresponding to the physical cover. ui1 and ui2 are the cover functions associated with P11 and P12 or P21 and P22 , respectively (An, 2010; Wu and Wong, 2012; Wu et al., 2013, 2017). When it turns to element e2, the three vertexes of this element are the stars of mathematical covers M3, M4 and M5. All of these mathematical covers are partially cut by a crack. Thus, element e2 is associated with P3 (Fig. 1(d)), P4 (Fig. 1(e)) and P5 (Fig. 1(f)), so Eq. (5) is not applicable anymore. Enrichment methods are proposed to solve this problem as
2. Basic concepts 2.1. Crack propagation criteria
3
3
Extensive studies have been conducted with the cracking initiation criterion (Biswas et al., 2012; Yang et al., 2014b) and crack propagation criterion (Grégoire et al., 2007; Zhang and Feng, 2011). In this study, the maximum circumferential stress criterion is adopted. The criterion of crack instability is
uj (x ) = ϕ·cj
1 θ cos 0 [KI (1 + cos θ0) − 3KII sin θ0] = KIC 2 2
where cj represents the array of additional unknowns.ϕ is the matrix of polynomial bases as
uh (x ) =
∑ φi (x )·ui (x ) + ∑ φi (x )·uj (x ) i=1
i=1
(6)
The uj (x ) is obtained by
(1) 2
(7)
Tunnelling and Underground Space Technology 92 (2019) 103032
L.F. Fan, et al.
Fig. 1. Basic concepts and crack modeling of NMM: a basic concepts and an internal crack, b completely cut of M1, c completely cut of M2, d partially cut of M3, e partially cut of M4, f partially cut of M5.
(An, 2010; Wu and Wong, 2012; Wu et al., 2013).
ϕ1 0 ϕ2 0 ϕ3 0 ϕ4 0 ⎤ ϕ=⎡ ⎢ 0 ϕ1 0 ϕ2 0 ϕ3 0 ϕ4 ⎥ ⎦ ⎣
(8)
where ϕ1, ϕ2 , ϕ3 and ϕ4 can be obtained by
θ [ϕ1, ϕ2 , ϕ3, ϕ4] = ⎡ r sin , 2 ⎣
θ r cos , 2
θ r sin θ sin , 2
2.3. Loop update
θ r sin θ cos ⎤ 2⎦
The loops should be updated once cracks propagated for the contact detection and crack presentation (Wu and Wong, 2012). The modification to the loops includes four types (more details in (Sagong and Bobet, 2002)). Once a new crack exists, the ME loop and PC loop would
(9) where (r, θ) are polar coordinates in the crack tip coordinate system 3
Tunnelling and Underground Space Technology 92 (2019) 103032
L.F. Fan, et al.
both change. Simple examples of ME loop update and PC loop update are given in Fig. 2 for explaining the update procedure. There are two types of loop update in the given examples. New crack would divide the original loop into two parts. Therefore the new loop is added (e.g. e2 to e5 and P1 to P3). On the other hand, the new crack may combine the original two loops into a loop (e.g. P2 to P4) (more details in (Sagong and Bobet, 2002; Wu and Wong, 2012)). 3. Stress wave propagation through intact rock bar The stress wave propagation through an intact rock bar was carried out to calibrate the present code. Fig. 3 shows a rock bar with a length of x = 2.00 m and a width of y = 0.10 m. The incident wave is achieved by dynamic loading at location A of xA = 0.0001 m as shown in Fig. 3. 10 evenly distributed measurement points at each location are arranged to avoid the measurement error. Three measurement locations are placed at xB = 0.50 m, xC = 1.00 m and xD = 1.50 m, respectively. The rock density is ρ = 2640 kg/m3, the Young’s modulus is E = 60 GPa and the Poisson ratio is ν = 0.3. In order to reduce lateral effects of stress wave and to satisfy the one-dimensional wave propagation theory, the length of the numerical model is ten times larger than the width. An incident stress wave of half sinusoidal waveform is applied at the left boundary of the rock bar as
0 ⩽ t ⩽ 1/2f0 σ sin(2πf0 t ) , when ⎧ σ =⎧ 0 ⎨ others ⎨ ⎩0 ⎩
(10)
where σ0 and f0 denote the amplitude and frequency of the incident wave, respectively. In the present study, σ0 = 10.0 MPa and f0 = 2000 Hz. Fig. 4 shows the simulation results. Since the stress wave transmits from the left of the bar to the right in this model, there is a time delay between the incident wave and each transmitted wave. The interval of two waves could be calculated by t = l/C, where l is the distance between the two corresponding measurement locations and C is the wave velocity obtained by C = E / ρ . It is also seen from Fig. 4 that the amplitudes of transmitted wave are 99.82%, 99.02% and 98.83% of the incident wave. The results show that the calculation induced error is within 1.20%, which is acceptable in the rock dynamic engineering. 4. Stress wave propagation through incline cracked rock bar 4.1. Dynamic stress field distribution The existence of crack leads to stress concentration around the tips of crack which results in the cracking of rock mass further influences the stress wave propagation. To investigate the effects of cracking on dynamic stress field, stress wave propagation through incline cracked rock bar is simulated as shown in the Fig. 5. To show the advantages of NMM in the simulation of stress wave propagation, the incline angles of the crack were selected at random. In the present study, angles of 13° and 31° were illustrated. The tensile strength and compressive strength of the rock bar are 30.0 MPa and 200.0 MPa, respectively. The cohesion is c = 15 MPa and the friction angle is φ = 30°. The parameters of the numerical model are selected according to the previous researches (Fan et al., 2013; Wu and Wong, 2012; Wu et al., 2013). Three rows of measurements are arranged with 10 points at location E, 25 points at location F and 25 points at location G, respectively. The measurement points located at locations F and G are parallel to the crack while the measurement points at location E are vertical. All points are numbered from the bottom to top subsequently (1 to 10 at location E and 1 to 25 at locations F and G, respectively). The horizontal locations of the crack tips are x = 0.82 m and x = 1.08 m and the angle of the crack is 13°.
Fig. 2. Update of ME loop and PC loop when cracking: a update of ME loop, b update of PC loop.
4
Tunnelling and Underground Space Technology 92 (2019) 103032
L.F. Fan, et al.
B
0.50 m
D
C
0.50 m
0.50 m
0.50 m
2.00 m
0.10 m
A
Fig. 3. Scheme of an intact rock model.
Fig. 4. Illustration of stress histories of measurement points in rock mass under dynamic non-equilibrium stress field.
The stress histories for the wave propagation through pre-existing oblique cracked rock bar are monitored and shown in Fig. 6. It is observed from Fig. 6(a) that a slight rise was observed at about 0.45 ms in Fig. 6(a). Since the main purpose of the present research is to study the stress field around the crack tips, the stress histories at points E-01 to E10 were studied. It is seen from Fig. 6(a) that the stress histories of ten measurement points E-01 to E-10 of x = 0.1 m show a good agreement with the amplitude of 1.0 MPa, which reveals that the lateral effects can be ignored in the present study. Moreover, it is also shown that the NMM can be utilized to simulate the one-dimensional wave propagation through rock mass with a convergent result. Fig. 6(b) and (c) show the stress distribution histories in front of and behind of the oblique crack, respectively. According to the previous researches under static loading (Li, 2010), the stress distribution around the crack is constant and in the form of centrosymmetric to the center of the crack. While for the stress wave propagation, the stress distribution varies as time increases. It is seen from Fig. 6(b) that the peak stresses around the tip in front of the crack (e.g. F-25) are larger than the amplitude of the incident wave (1.0 MPa). Therefore, stress concentration around crack tips also can be observed under dynamic loading. However, different phenomenon from static loading, it is also seen from Fig. 6(b) that the stress history around the bottom tip in front of the crack (e.g. F-01) are with different waveforms, which is induced by the reflected stress waves. Moreover, obvious tensile stress field (e.g. F-01, 02, 03 and so on) can be observed around the bottom tip in front of the crack.
Fig. 6. Stress distribution inside the rock model without cracking under dynamic non-equilibrium stress field.
The stress history behind the crack shows similar waveforms as the incident wave as shown in Fig. 6(c). It is seen from Fig. 6(c) that distinct stress concentration also can be observed around both the bottom and top tips of the crack (e.g. G-01 and G-25). The concentration rate around the bottom tip behind the crack is larger than that around the top tip behind the crack, which is in accordance with the previous researches under static loading (Li, 2010). It should also be noted that unlike the stress wave propagation through crack normally, the times associated with the peak stresses are different, which are induced by the oblique geometrical characters of the pre-existing crack results in the different wave propagation time. To investigate the stress field variation under dynamic loading, the
Fig. 5. Scheme of a rock model with a long crack and oblique measurement points. 5
Tunnelling and Underground Space Technology 92 (2019) 103032
L.F. Fan, et al.
transmitted waveforms recorded at location H convergent before they were influenced by the cracking. It is seen from Fig. 10 that tensile stress was recorded, which was reflected by the crack. The phenomenon is in accordance with the previous researches. However, it is interesting that the amplitude of the tensile stress decreases for the measurement points from H-01 to H-10, which is caused by the cracking induced dynamic energy loss. Such results cannot be simulated in the previous researches using the traditional manifold methods (Fan et al., 2013). It is also seen from Fig. 10 that the crack propagation has significant effects on the stress field at location I. It is seen that the stress waves at ten measurement points I-01 to I-10 convergent to each other before cracking, which is at about 0.52 ms for the present study. The amplitude of stress waves at measurement points from I to 01 to I-10 decreases after cracking. The phenomenon is similar as the effects of cracking on the wave propagation, which also could not be simulated by the traditional NMM (Fan et al., 2013). 5. Case study
Fig. 7. Inhomogeneous stress field around the crack under dynamic non-equilibrium stress field.
To validate the engineering application potentials of the developed NMM, two cases of underground cavern explosion with symmetrical and unsymmetrical discontinuities under dynamic stress wave propagation were investigated.
stress distribution at four typical times are investigated and shown in Fig. 7. Four typical times are selected as the times of the stresses around top tip and bottom tip in front of the crack, and top tip and bottom tip behind the crack reach their peaks as shown in the Fig. 6(b) and (c), respectively. It is seen from Fig. 7 that the stress increases from the bottom tip to top tip in front of the crack, while decreases from the bottom tip to top tip behind the crack in general. The variation of stress distribution is similar as that under static loading at times at peaks of F25, G-01 and G-25. The inflection point of stress distribution is not located at the tips of the crack, while very close to the tips, such as measurement point 4 at times at peaks of F-01, F-25, G-01 and G-25, and measurement point 22 at times at peaks of F-25, G-01 and G-25, respectively. In addition, it is also seen from Fig. 7 that the stress distribution is not centrosymmetric, which can be observed at time at peak of F-01. The phenomenon is caused by the reflected stress wave from the inclined crack, which is agreed with the results from Fig. 6.
5.1. Symmetrical cracking Fig. 11 shows the cross-section of a test site. The tunnel is circular shaped with a radius of 5.0 m and located at 50 m under the ground. To decrease the boundary effects, the cavern is in the center of a square with a length of 100 m. Displacements at three boundaries are fixed as shown in the Fig. 11. The properties of rock are the same as mentioned in the Section 4. A crack with a length of 27.94 m is located at the ramp above the cavern. The distances of two ends of the crack to the center point of the cavern are equal, namely, l1 = l2. Therefore, the stress wave arrives at both endpoints of the crack simultaneously and the stress value is equal. In the present model, the measurement points of row A and row C are horizontally located while those of row B are vertically located. There are 14 measurement points in each row and the space between two neighboring points is 2.0 m. The measurement points in each row are numbered as 1–14 from the center to the border. A dynamic loading is applied on the internal boundary of the tunnel to simulate a normally incident wave. A stress wave with multiple peaks of different amplitudes is applied with a waveform of
4.2. Cracking mode The failure criteria as shown in Eqs. (1) and (2) is introduced in the present NMM code. Fig. 8 shows the numerical model to investigate the cracking mode and the stress distribution variation under dynamic loading. An initial inclined crack starts at x = 0.90 m and ends at x = 1.00 m with an incline angle of 31° is assumed. The measurement points located at xH = 0.3 m and xI = 1.3 m are set to analyze the effects of cracking on the wave propagation, respectively. Each location with 10 equal-spaced measurement points as in the previous sections. Fig. 9(a) and (b) show a visual result of dynamic crack propagation. The cracking process can be divided into two stages. The first stage is that the crack tip closer to the dynamic loading starts to propagate when the stress reaches its critical value while the farther crack tip maintains the original state as shown in Fig. 9(a). The second stage is that the other crack tip propagates when the stress wave arrives the right tip as shown in Fig. 9(b). Fig. 10 gives the measuring data at locations H and I. The
⎧ 0 ⩽ t ⩽ 1/(2f1 ) ⎧ σ1 sin(2πf1 × t ) σI = σ2 sin(2πf2 × t ) , when 1/(2f1 ) + Δ T⩽ t ⩽ 1/(2f2 ) + 1/(2f1 ) + ΔT ⎨ ⎨ ⎩0 ⎩ others (11) where ΔT denotes the time period between two pulses of the wave. σ1 and σ2 are the amplitudes of the wave with multiple peaks, respectively. f1 and f2 are the corresponding frequencies. The numerical result of two symmetrically wing cracks with equal distances from the crack tips could be seen in Fig. 12. For a symmetrical crack case, the two crack tips start to propagate at the same time and the crack paths are almost symmetrical. The new cracks both through
Fig. 8. Scheme of a rock model with a short crack and vertical measurement points. 6
Tunnelling and Underground Space Technology 92 (2019) 103032
L.F. Fan, et al.
Fig. 9. Crack propagation under dynamic non-equilibrium stress field: a crack propagation at left tip of the crack, b crack propagation at right tip of the crack.
Fig. 11. Scheme of an underground cavern with an equidistant crack.
incident wave, respectively. Different from the one-dimensional wave propagation through rock mass, which is caused by the micro-defect induced attenuation, the attenuation in the present study is mainly caused by the wave front circular expansion, which reduces the dynamic energy density further results in the wave attenuation. Moreover, it is also can be observed from Fig. 13 that tensile stress is introduced between two transmitted pulses, which is caused by the compression stress wave applied on the internal surface of the circle carven. Fig. 10. Stress histories of measurement points during the dynamic crack process.
5.2. Unsymmetrical cracking Just as the symmetrical cracking case, an unsymmetrical crack with a length of 30.04 m is established as shown in Fig. 14. The distances from two crack tips to the center point of the cavern are unequal (e.g. l3 < l4.). Therefore, the stress wave reaches the two tips in sequence. The boundary conditions, material properties and dynamic loading are the same as mentioned in the Section 5.1. Fig. 15 shows the numerical result of crack propagation. It is seen from Fig. 15 that the cracking initiates unsymmetrically, which is
the measurement points of 8th and 9th of row A and row C which influence the wave propagation, respectively. Fig. 13 shows the stress wave histories. The transmitted wave is the stress history recorded at the measurement point A-1 (the 1st measurement point of row A). It can be observed from Fig. 13 that both first pulse and second pulse attenuate. The amplitudes of the first pulse and second pulse are 58.54% and 50.49% comparing to amplitudes of the 7
Tunnelling and Underground Space Technology 92 (2019) 103032
L.F. Fan, et al.
Fig. 14. Scheme of an underground cavern with a non-equidistant crack.
Fig. 12. Simulation result of equidistant crack propagation under dynamic nonequilibrium stress field.
Fig. 13. Waveforms of the incident wave and a typical transmitted wave for underground cavern explosion.
significantly different from the case of cracking under static loading. The crack tip closer to the wave source propagates while the other crack tip keeps constant. The phenomenon can be explained by the dynamic loading attenuation and the interaction between the dynamic loading and cracking. Since the wave front circle enlarges as the stress wave propagates from the wave source, the stress amplitude decreases as the distance from the wave source increases. Therefore, the farther tip of crack may not satisfy the cracking criteria even the closer tip of crack propagates because of the wave attenuation. Moreover, the cracking of the closer tip absorbs the dynamic energy of stress wave, which further increases the wave attenuation. A stress wave as shown in Fig. 13 is applied. Fig. 16 shows the amplitudes of the recorded stress histories at the measurement locations D, E and F. It is seen from Fig. 16 that the new crack passes through the measurement points 8th and 9th of row D according to the simulation result. The second pulse at row D attenuated after propagation of the new crack. Due to the unilateral cracking, the measurement points at row E
Fig. 15. Simulation result of non-equidistant crack propagation under dynamic non-equilibrium stress field.
and row F are not influenced by the pre-existing crack and the new crack. They are agreed well with each other as shown in Fig. 16. In contrast, there is a distinct decrease of stress at row D when the stress wave propagates through the new crack. 6. Conclusions The present study investigates the interaction between dynamic stress wave propagation and the cracking behavior using the NMM. The stress field around the crack tips, cracking modes of pre-existing crack and the stress field after crack propagation are analyzed. The NMM has its advantage in investigating the cracking behavior and dynamic stress field variation inside the rock mass under stress wave. It can be concluded that a transient stress field would be informed around the crack tips during the stress wave propagation. The dynamic 8
Tunnelling and Underground Space Technology 92 (2019) 103032
L.F. Fan, et al.
adjacent tunnel explosion. Tunn. Undergr. Sp. Tech. 35, 227–234. Li, X., Li, X.F., Zhang, Q.B., et al., 2018. A numerical study of spalling and related rockburst under dynamic disturbance using a particle-based numerical manifold method (PNMM). Tunn. Undergr. Sp. Tech. 81, 438–449. Ma, G.W., An, X.M., He, L., 2010. The numerical manifold method: A review. Int. J. Comput. Methods 7 (1), 1–32. Ning, Y., An, X., Ma, G., 2011. Footwall slope stability analysis with the numerical manifold method. Int. J. Rock. Mech. Min. 48 (6), 964–975. Nejati, H.R., Ghazvinian, A., 2014. Brittleness effect on rock fatigue damage evolution. Rock Mech. Rock Eng. 47 (5), 1839–1848. Rannou, J., Gravouil, A., Baïetto-Dubourg, M.C., 2009. A local multigrid x-fem strategy for 3-d crack propagation. Int. J. Numer. Meth. Eng. 77 (4), 581–600. Rossmanith, H.P. (Ed.), 2014. RockOCK FRACT. MECH., vol. 275 Springer. Shi, G., 1992. Modeling rock joints and blocks by manifold method. In: Proceedings of the 33rd US Rock Mechanics Symposium. New Mexico, pp. 639–648. Shen, B., Stephansson, O., Einstein, H.H., et al., 1995. Coalescence of fractures under shear stresses in experiments. J. Geophy. Res. 100 (B4), 5975–5990. Strouboulis, T., Babuška, I., Copps, K., 2000. The design and analysis of the generalized finite element method. Comput. Method. Appl. M. 181 (1), 43–69. Sagong, M., Bobet, A., 2002. Coalescence of multiple flaws in a rock-model material in uniaxial compression. Int. J. Rock. Mech. Min. 39 (2), 229–241. Wang, T.J., 1991. A continuum damage model for ductile fracture of weld heat affected zone. Eng. Fracture Mech. 40 (6), 1075–1082. Wang, T.J., 1992a. Unified CDM model and local criterion for ductile fracture-I. Unified CDM model for ductile fracture. Eng. Fracture Mech. 42 (1), 177–183. Wang, T.J., 1992b. Unified CDM model and local criterion for ductile fracture-II. Ductile fracture local criterion. Eng. Fracture Mech. 42 (1), 185–192. Wong, R.H.C., Chau, K.T., 1998. Crack coalescence in a rock-like material containing two cracks. Int. J. Rock. Mech. Min. 35 (2), 147–164. Wong, L.N.Y., Einstein, H.H., 2009a. Crack coalescence in molded gypsum and Carrara marble: Part 1. Macroscopic observations and interpretation. Rock Mech. Rock Eng. 42 (3), 475–511. Wong, L.N.Y., Einstein, H.H., 2009b. Systematic evaluation of cracking behavior in specimens containing single flaws under uniaxial compression. Int. J. Rock Mech. Min. 46 (2), 239–249. Wu, Z.J., Fan, L.F., Liu, Q.S., Ma, G.W., 2017. Micro-mechanical modeling of the macromechanical response and fracture behavior of rock using the numerical manifold method. Eng. Geol. 225 (20), 49–60. Wu, Z., Wong, L.N.Y., 2012. Frictional crack initiation and propagation analysis using the numerical manifold method. Comput. Geotech. 39 (1), 38–53. Wu, Z., Wong, L.N.Y., Fan, L., 2013. Dynamic study on fracture problems in viscoelastic sedimentary rocks using the numerical manifold method. Rock Mech. Rock Eng. 46 (6), 1415–1427. Yang, S.Q., Jing, H.W., 2011. Strength failure and crack coalescence behavior of brittle sandstone samples containing a single fissure under uniaxial compression. Int. J. Fract. 168 (2), 227–250. Yang, S.Q., Huang, Y.H., Jing, H.W., et al., 2014a. Discrete element modeling on fracture coalescence behavior of red sandstone containing two unparallel fissures under uniaxial compression. Eng. Geol. 178, 28–48. Yang, L., Yang, R., Qu, G., et al., 2014b. Caustic study on blast-induced wing crack behaviors in dynamic–static superimposed stress field. Int. J. Min. Sci. Tech. 24 (4), 417–423. Yang, Y.T., Tang, X.H., Zheng, H., et al., 2018. Hydraulic fracturing modeling using the enriched numerical manifold method. Appl. Math. Model. 53, 462–486. Wu, Z.J., Zhang, P.L., Fan, L.F., 2019a. Numerical study of the effect of confining pressure on the rock breakage efficiency and fragment size distribution of a TBM cutter using a coupled FEM-DEM method. Tunn. Undergr. Sp. Tech. 88, 260–275. Wu, Z.J., Zhou, Y., Fan, L.F., 2019b. A fracture aperture dependent thermal-cohesive coupled model for modelling thermal conduction in fractured rock mass. Comput. Geotech. 114, 103108. Xia, X., Li, H.B., Lu, Y.Q., 2018. A case study on the cavity effect of a water tunnel on the ground vibrations induced by excavating blasts. Tunn. Undergr. Sp. Tech. 71, 292–297. Xie, L.X., Lu, W.B., Zhang, Q.B., et al., 2018. Analysis of damage mechanisms and optimization of cut blasting design under high in-situ stresses. Tunn. Undergr. Sp. Tech. 66, 19–33. Zhu, W., Tang, C., 2006. Numerical simulation of Brazilian disk rock failure under static and dynamic loading. Int. J. Rock Mech. Min. 43 (2), 236–252. Zhao, X., Zhao, J., Cai, J., et al., 2008. UDEC modelling on wave propagation across fractured rock masses. Comput. Geotech. 35 (1), 97–104. Zhang, H., Li, L., An, X., et al., 2010. Numerical analysis of 2-d crack propagation problems using the numerical manifold method. Eng. Anal. Bound. Elem. 34 (1), 41–50. Zhang, Y., Feng, X., 2011. Extended finite element simulation of crack propagation in fractured rock masses. Mater. Res. Innov. 15 (s1), s594–s596. Zhang, G., Xu, L., Haifeng, L., 2015. Simulation of hydraulic fracture utilizing numerical manifold method. Sci China Tech. Sci. 58 (9), 1542–1557. Zhao, G.F., Xia, K., 2018. A study of mode-I self-similar dynamic crack propagation using a lattice spring model. Comput. Geotech. 96, 215–225. Zhao, G.F., Hu, X., Li, Q., et al., 2018. On the four-dimensional lattice spring model for geomechanics. J. Rock Mech. Geotech. Eng. 10 (4), 661–668. Zhao, G.F., Yin, Q., Russell, A.R., et al., 2019. On the linear elastic responses of the 2D bonded discrete element model. Int. J. Number. Anal. Met. 43 (1), 166–182. Zhu, J.B., Deng, X.F., Zhao, X.B., Zhao, J., 2013. A numerical study on wave transmission across multiple intersecting joint sets in rock masses with UDEC. Rock Mech. Rock Eng. 46 (6), 1429–1442. Zhou, X.F., Fan, L.F., Wu, Z.J., 2017. Effects of microfracture on wave propagation through rock mass. Int. J. Geomech. 17 (9), 04017072.
Fig. 16. Peak stress variation before and after crack propagation of non-equidistant crack under dynamic non-equilibrium stress field.
stress field is different from that under static loading. The transient stress field is induced by both stress concentration and the reflection and transmission of stress wave. The characteristics of transient stress field result in special cracking modes. The crack tip closer to the loading points propagates earlier than the farther one during the stress wave propagation. In addition, the cracking process further influences the stress field distribution. Acknowledgements The research is supported by the National Natural Science Foundation of China (NSFC) (Nos. 51778021, 11572282 and 41772309) and the Science Fund for Creative Research Groups of the National Natural Science Foundation of China (No. 51421005). References An, X.M., 2010. Extended numerical manifold method for engineering failure analysis. Ph.D thesis. Technological University. Singapore, Nanyang. Biswas, N., Ding, J., Balla, V.K., et al., 2012. Deformation and fracture behavior of laser processed dense and porous ti6al4v alloy under static and dynamic loading. Mater. Sci. Eng. A. 549, 213–221. Deng, X.F., Zhu, J.B., Chen, S.G., Zhao, Z.Y., Zhou, Y.X., Zhao, J., 2014. Numerical study on tunnel damage subject to blast-induced shock wave in jointed rock masses. Tunn. Undergr. Sp. Tech. 43, 88–100. Dong, Q., Li, X., Zhao, H., 2014. Experimental research on ultrasonic p-wave velocity variation of fractured rock mass under different stress paths. Adv. Mater. Res. 1065–1069, 35–39. Fan, L.F., Wang, L.J., Wu, Z.J., 2018. Wave transmission across linearly jointed complex rock masses. Int. J. Rock Mech. Min. 112, 193–200. Fan, L.F., Wu, Z.J., Wan, Z., Gao, J.W., 2017. Experimental investigation of thermal effects on dynamic behavior of granite. Appl. Therm. Eng. 125, 94–103. Fan, L., Yi, X., Ma, G., 2013. Numerical manifold method (NMM) simulation of stress wave propagation through fractured rock mass. Int. J. Appl. Mech. 05 (02), 1350022. Grégoire, D., Maigre, H., Rethore, J., et al., 2007. Dynamic crack propagation under mixed-mode loading-comparison between experiments and X-FEM simulations. Int. J. Solids. Struct. 44 (20), 6517–6534. Hoek, E., Bieniawski, Z.T., 1965. Brittle fracture propagation in rock under compression. Int. J. Fract. Mech. 1 (3), 137–155. Hu, X.D., Zhao, G.F., Deng, X.F., et al., 2018. Application of the four-dimensional lattice spring model for blasting wave propagation around the underground rock cavern. Tunn. Undergr. Sp. Tech. 80, 135–147. Irwin, G.R., 1997. Analysis of stresses and strains near the end of a crack traversing a plate. Spie Milestone Series MS 137 (167–170), 16. Leucci, G., Giorgi, L.D., 2006. Experimental studies on the effects of fracture on the p and s wave velocity propagation in sedimentary rock (“calcarenite del salento”). Eng. Geol. 84 (3), 130–142. Li, S.Y., 2010. Introduction to Rock Fracture Mechanics. China University of Science and Technology. Lee, H., Jeon, S., 2011. An experimental and numerical study of fracture coalescence in pre-cracked specimens under uniaxial compression. Int. J. Solid. Stut. 48 (6), 979–999. Li, J.C., Li, H.B., Ma, G.W., et al., 2013. Assessment of underground tunnel stability to
9