- Email: [email protected]
Probabilistic modeling of excavation-induced damage depth around rock-excavated tunnels
Journal Pre-proof Probabilistic Modeling of Excavation-Induced Damage Depth around Rock-Excavated Tunnels Mahdi Shadabfar, Hongwei Huang, Hadi Kordestani, Edmond V. Muho PII:
S2590-1230(19)30075-1
DOI:
https://doi.org/10.1016/j.rineng.2019.100075
Reference:
RINENG 100075
To appear in:
Results in Engineering
Received Date: 10 November 2019 Accepted Date: 20 November 2019
Please cite this article as: M. Shadabfar, H. Huang, H. Kordestani, E.V. Muho, Probabilistic Modeling of Excavation-Induced Damage Depth around Rock-Excavated Tunnels, Results in Engineering, https:// doi.org/10.1016/j.rineng.2019.100075. This is a PDF file of an article that has undergone enhancements after acceptance, such as the addition of a cover page and metadata, and formatting for readability, but it is not yet the definitive version of record. This version will undergo additional copyediting, typesetting and review before it is published in its final form, but we are providing this version to give early visibility of the article. Please note that, during the production process, errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain. © 2019 The Author(s). Published by Elsevier B.V.
Probabilistic Modeling of Excavation-Induced Damage Depth around Rock-Excavated Tunnels Mahdi Shadab Fara,∗, Hadi Kordestanib , Edmond V. Muhoc a
Department of Geotechnical Engineering, Tongji University, Shanghai, China Department of Civil Engineering, Qingdao University of Technology, Qingdao, China c Department of Disaster Mitigation for Structures, College of Civil Engineering, Tongji University, Shanghai, China b
Abstract The cracked zones caused by tunnel excavations are among the major reasons for structural instability in mechanized tunnel projects. Thus, various methods have been proposed to evaluate such cracks in terms of size and depth. However, deterministic estimations have been associated with a margin of error due to the inherent uncertainties of the respective problems as well as the approximation nature of the existing methods. Therefore, in this study, a probabilistic approach was adopted to assess the depth of a cracked zone induced around a tunnel. For this purpose, the problem was formulated as a reliability model, and the Monte Carlo method was utilized to calculate the exceedance probability for the crack depth around the tunnel. The exceedance probability was then calculated for any different depths of crack and presented as probability contours around the excavation zone. It was demonstrated that the exceedance probability is reduced sharply with an increase in the crack depth such that it fell below 5% for the cracks deeper than 2.28 times the tunnel radius. Effects of the maximum tangential stress at excavation boundary and crack initiation threshold on the exceedance probability were further calculated using the reliability sensitivity analysis. Eventually, the influence of in-plain stress ratio on the failure probability was investigated, and the results were reported for three rock types, including limestone, granite, and mudstone. Keywords: Excavation-induced damage (EDZ), uncertainty modelling, Monte Carlo sampling method, failure probability, reliability sensitivity analysis. 1. Introduction In the course of excavation, redistribution of the in-situ stresses in the surrounding rock mass leads to a gradual increase in the tangential stresses while decreasing the radial stresses [1–4]. This develops an extreme state at the opening surface, causing a stress concentration ∗
Corresponding author Email addresses: [email protected] (Mahdi Shadab Far ), [email protected] (Hadi Kordestani), [email protected] (Edmond V. Muho) Preprint submitted to Results in Engineering
November 10, 2019
around the excavation area [5]. When the resulting stress exceeds the failure strength of the rock mass, excavation-induced cracks are initiated and gradually penetrate the rock mass until a new stress equilibrium is reached [6, 7]. Ultimately, these single cracks get interconnected to form damage layers with different crack densities. The first damage layer, which serves as host for very dense and interconnected cracks, is called the high damage zone (HDZ) [8]. Moving outward, layers of lower crack density but connected crack systems become dominant. These layers are collectively referred to as the inner excavation damage zone (EDZi). Then, the cracks diverge into joints in an outer excavation damage zone (EDZo) [9]. Beyond the EDZs, there is an elastic and damage-free yet affected zone called the excavation disturbed zone (EdZ) or the excavation influence zone (EIZ). The EIZ is preferred to avoid possible confusions between the EDZ and the EdZ [10, 11]. The described zones are schematically shown in Figure 1. EIZ EDZo
EDZi Excavation area
HDZ
Figure 1: Different affected zones around an excavation area.
The generated damage zones threaten the stability and performance of underground structures and are reportedly the source of major geological hazards, such as spalling, desaturation, and roof collapse [12, 13]. Accordingly, numerous studies have been conducted to investigate the initiation, propagation, extension, dimensions, and induced permeability of the EDZs [14–21]. In this respect, Martini et al. [22] examined the formation and propagation of the EdZ by recording data around a large-scale underground opening. Fakhimi et al. [23] performed biaxial compression tests on sandstone specimens to simulate the failure around a circular opening in brittle rocks. Similarly, Perras and Diederichs [11] used damage initiation and spalling limit (DISL) and proposed a guideline to predict the depths of EDZs in brittle rocks. In line with the above-mentioned models, several other empirical methods and numerical models are available in the literature for addressing the excavation-induced damage. Considering the assumptions of the Hoek-Brown’s Failure Criterion and the conceptual model of Diederichs for rock behavior, Hrubesova and Duris [24] examined three damage zones, namely HDZ, EDZ, and EIZ, around an explosion site and applied the same methodology to a real case study. Dadashzadeh and Diederichs [25] employed cohesion weakening 2
and frictional strengthening model and developed a finite-difference method in FLAC 2D software to predict variations of the EDZ depth in brittle rocks. Tang et al. [26] presented a modified failure criterion using a series of cyclic loading tests to formulate an analytical solution for predicting different EDZs. Combining the discrete fracture networks and a twodimensional finite-discrete element method (FDEM) code, Vazaios et al. [27] investigated the effects of joints on the reduction of the rock strength and bearing capacity to accommodate damage zones. In another study [28], they adopted the FDEM to investigate rock mass behavior and examine the spalling-induced damage extent. Xie and Peng [29] used an ultrasonic velocity detector to assess the range of EDZs in a rock medium during a drilling process and evaluated the damage extension from the drilling data. Tao et al. [5] developed an analytical model for estimating the extent of EDZ and provided a closed-form solution for the established problem. Then, they conducted a parametric study to investigate the effects of different variables on the induced damage. The results revealed that the geological strength index and uniaxial compressive strength imposed the largest contributions to the damage rates. However, possible disturbances (e.g., mechanical drilling and the waves generated in different stages of the drilling operation) could significantly affect the problem, even dominating the effects of the parameters involved in mathematical modeling in some cases. To the authors’ best knowledge, most of the studies on the excavation-induced damage size have adopted a deterministic approach for formulating the problem. However, the uncertainties associated with the stress conditions and crack resistance of the rock are so influential that cannot be neglected. In fact, deterministic modeling is equivalent to assuming 100% certainty about the load and bearing capacity of the rock medium, which makes no sense when it comes to real cases. Therefore, an alternative method is required to take the uncertainty into consideration and express the problem probabilistically. In the present study, reliability modeling is used to address this issue. To this end, the input parameters are initially modeled as random variables with certain means and standard deviations. Then, a limit state function is constructed based on the model developed by Perras and Diederichs [11] and the depths of the EDZs are expressed as a reliability problem. The Monte Carlo sampling method is then utilized to solve the established problem and the results are reported as exceedance probability at different depths of the EDZs and presented as probability contour maps around the tunnel. Effects of different parameters on the failure probability are also investigated through reliability sensitivity analysis. Finally, the effect of in-plane stress ratio is demonstrated by exceedance probability curves and the results are reported for different rock types. 2. Modeling As previously outlined in the Introduction section, this study attempts to express the formation of excavation-induced damage zone as a probabilistic problem, because deterministic models are mainly limited to specific solutions even for real cases encountering high levels of uncertainty. More specifically, a deterministic model assumes 100% probability for the occurrence of a particular damage depth, which makes it highly unrealistic. A more ra3
tional solution to this problem is to calculate exceedance probabilities for different probable damage depths; i.e., any damage depth corresponds to a particular value of exceedance probability. For this purpose, one must begin with defining a deterministic model to calculate the crack depth via an explicit relationship. Then, the deterministic model can be converted to a probabilistic one by setting the input parameters to random variables, constructing a reliability problem. Later, a conventional reliability analysis method can be adopted to solve the established problem (here, we use the Monte Carlo sampling method for this purpose). To this end, the Perras and Diederichs model [11] was considered as the primary deterministic model. In the Perras and Diederichs study [11], the natural variations of brittle rock were considered explicitly by incorporating a laboratory model into numerical simulations with brittle constitutive behavior. The results of all numerical modeling (different rock types, stress conditions, etc.) were then collected, and the best average fit for each damage zone was developed using regression models. As a result, a simple but reasonable estimation of normalized damage depth in three different damage classes (i.e., EDZi, EDZo, and HDZ) was provided covering a wide range of possible scenarios: σmax rd = 1 + B( − 1)D , (1) R CI where rd represents the depth of inner EDZ (EDZi), outer EDZ (EDZo), or highly damaged zone (HDZ), R denotes the tunnel radius, σmax presents the maximum tangential stress at excavation boundary, and CI is regarded as the crack initiation threshold. B and D are constant regression parameters that are available for different rock types, in-plane stress ratios (KHh ), and types of damaged zones (i.e., EDZi, EDZo, or HDZ). For instance, the values of B and D are 0.66 and 0.63, respectively, when predicting the depth of EDZo in limestone with KHh = 1.5. Denoting the excavation influence on the rock as S and the bearing capacity of the rock as C, the limit state function (LSF ) can be expressed as Eq. (2). If the disturbance caused by the excavation overshadows the rock bearing capacity (i.e., C is less than S), then the model faces damage. Otherwise, if C is greater than S, then the model is assumed to be safe. LSF = C − S.
(2)
When the rock is highly resistant to excavation-induced disturbance, the depth of cracks is very limited. However, if the rock is too weak to withstand such disturbances, the generated cracks penetrate more into the rock. As a result, the crack depth can well represent the rock bearing capacity. Hence, assuming that S is the estimated normalized depth of damage generated around the excavation hole (e.g., the value obtained from the Perras and Diederichs model [11]) and C is a specified normalized depth of crack (e.g., rd /R = 21 ), which is taken as the boundary of the model failure, the LSF is rewritten as follows: 1
The other possible values will be considered in Section 3.2.
4
Table 1: Random variables assumed in the analysis.
Rock type
Variable Mean σmax (MPa) 72 Limestone CI (MPa) 41
Standard deviation 21 12
σmax − 1)D . (3) CI Here, a major objective is to determine how much probability is needed to face damages with depths beyond C. This is equivalent to the probability of S > C, which can be calculated as follows: LSF = C − 1 − B(
Failure Probability = P (LSF < 0) = P (S > C).
(4)
Various methods have been proposed to estimate this probability [30, 31]. Considering its advantages [32, 33], the Monte Carlo sampling method was used in the present paper. To use the Monte Carlo method, sufficient numbers of random values should be generated for all variables involved in the problem. The random values should correspond to the mean and standard deviation of the respective variable and follow the governing distribution function. By substituting the generated random values into the LSF , a so-called LSF random sample was achieved for each set of random values. Counting the number of less-than-zero samples and dividing it by the total number of samples, the desired exceedance probability was calculated. To identify the optimal number of the random values (sample size), the sample size was increased gradually and the resultant changes in the exceedance probability were evaluated. This was continued until the probability was no longer sensitive to the samples size and the number of samples was high enough to ensure the convergence. Then, the calculated failure probability was recorded as the output of the Monte Carlo method. The process of using the Monte Carlo method, from selecting the primary model to the calculation of the exceedance probability, is schematically shown in Figure 2. 3. Analysis 3.1. Reliability analysis To start using the Monte Carlo sampling method, the input parameters, including the maximum tangential stress and crack initiation threshold, were defined as random variables. In this way, uncertainties are incorporated into the model parameters by considering a relatively wide range of possible caseloads in the analysis process. Here, we assumed that the variables followed normal distribution functions with means and standard deviations presented in Table 1. This assumption is in line with the values presented in the Perras and Diederichs model [11]. Continuing with the algorithm, 20,000 uniformly-distributed random values were generated for each of the considered variables. Then, using the inverse CDF of the normal distribution function (Eq. 5), the generated random values were transferred from uniform 5
Establish a primary deterministic model
Define the random variables
Define the limit state function (LSF)
Calculate the LSF based on the generated random numbers
Calculate the exceedance probability
NO
Monte Carlo loop
Increase the number of random samples
Generate random numbers for each random variable
Is the analysis converged? YES Record the probability, Pf
Figure 2: Schematic view of the Monte Carlo algorithm.
distribution to normal distribution with the probabilistic characteristics specified in Table 1. Histograms of the generated random values are shown in Figure 3 (blue column bars). zi = Φ−1 (ui ).
(5)
Using the generated random values, the rd /R ratio was calculated based on Eq. 1 for limestone with KHh = 1.5, and the proposed LSF was estimated through Eq. 3. The exceedance probability was then calculated via the following formula: n , (6) N where n is the number of cases for which LSF < 0 and N represents the total number of random values. In the current analysis, there were 2,354 less-than-zero cases among the whole random samples. Therefore, the exceedance probability was calculated as 2354/20000 = 0.1177 = 11.77%. To verify the results and ensure the convergence of the analysis, sensitivity of the results to the sample values was investigated. For this purpose, exceedance probability was calculated for different sample sizes and plotted on a graph. The results are illustrated in Figure 4. As can be seen in Figure 4, the results highly fluctuate in smaller sample sizes, so that P =
6
15
Histogram of σmax (MPa) Cases where LSF<0 Fitted distribution
Relative frequency (%)
Relative frequency (%)
15
10
5
00
10
5
0
20 40 60 80 100 120 140 160 Maximum tangential stress at the excavation boundary, σmax (MPa)
Histogram of CI (MPa) Cases where LSF<0 Fitted distribution
0
10
20
30 40 50 60 70 Crack initiation threshold, CI (MPa)
(a)
80
90
100
(b)
Probability of exceedance (%) for rd/R=2
Figure 3: Histograms of the generated random values for (a) σmax and (b) CI.
24 22
Probability of exceedance Final value
20 18 16 14 12 10 8
60
0.2
0.4
0.6
0.8 1 1.2 Sample size
1.4
1.6
1.8
2 × 104
Figure 4: Exceedance probability versus sample size.
a slight change in the sample size could make a considerable change in the exceedance probability. However, as the sample size increased, the change in the exceedance probability became smaller and the probability converged to 11.77% with negligible variations for sample sizes beyond 10,000. 3.2. Exceedance probability curve Thus far, the exceedance probability was calculated for C = rd /R = 2. However, the purpose of the present study is to calculate the exceedance probability for any desired value of normalized damage depth. Therefore, the rd /R was defined as a decision variable rather than a constant value. Then, assuming different values for the rd /R, the established problem was solved and the exceedance probability was calculated accordingly. The results are shown in Table 2 for different values of rd /R. The exceedance probability was plotted versus different values of rd /R and a curve was fitted to the results (black line in Figure 5). This curve, known as the exceedance probability 7
Table 2: Exceedance probability for different decision variables.
rd /R 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 2 2.2 2.4 2.6 2.8 3
Less-than-zero number P = n/N 19996 0.9998 19996 0.9998 19996 0.9998 19921 0.9961 17997 0.8999 16730 0.8365 13536 0.6768 8881 0.4441 4824 0.2412 2354 0.1177 1154 0.0577 592 0.0296 325 0.0163 206 0.0103 132 0.0066
P (%) rd /R 99.98 3.2 99.98 3.4 99.98 3.6 99.61 3.8 89.99 4 83.65 4.2 67.68 4.4 44.41 4.6 24.12 4.8 11.77 5 5.77 5.2 2.96 5.4 1.63 5.6 1.03 5.8 0.66 6
Less-than-zero number P = n/N 87 0.0044 65 0.0033 54 0.0027 46 0.0023 35 0.0018 25 0.0013 21 0.0011 19 0.0010 16 0.0008 15 0.0008 11 0.0006 9 0.0005 8 0.0004 8 0.0004 8 0.0004
P (%) 0.44 0.33 0.27 0.23 0.18 0.13 0.11 0.10 0.08 0.08 0.06 0.05 0.04 0.04 0.04
curve, shows the probability of exceedance for different values of normalized damage depth. The obtained exceedance probability curve has practical applications in the reliability-based evaluation of real projects where a target probability level (allowable probability) is set and the corresponding crack depth is determined. A typical numerical analysis is then required to establish the crack dimensions in the real project. If the size of the crack is greater than the allowable value, the expected probability is violated and, hence, further measures must be taken to reduce failure probability. Otherwise, the project is safe and the implementation process can be continued. Moreover, this graph shows the rock capacity to withstand the disturbances caused by excavation. If the diagram is generally high, the cracks are more likely to grow deep into the rock and, consequently, the failure probability is high. On the contrary, if the diagram is generally low, the cracks cannot easily grow through the rock and, thus, the rock medium exhibits a larger capacity to accommodate the excavation effects. 3.3. Data filtering In the previous section, the convergence of exceedance probability was ensured and it was demonstrated that the given probability was not sensitive to the number of samples. Nevertheless, a brief overview of the exceedance probability diagram presented in Figure 5 by black line reveals an error in this graph. It is expected that the probability value is 100% for rd /R ≤ 1 and the exceedance probability diagram begins to decline from rd /R = 1; however, this trend is not observed on the diagram. Looking at the definition of LSF in Eq. 3, it can be inferred that nonpositive values of (σmax /CI − 1) lead to the introduction of non-real values into the random samples, thereby affecting the calculated final probability value. For this reason, it should be ensured that: σmax − 1 ≥ 0 ⇒ σmax ≥ CI. CI 8
(7)
100 Before filtering After filtering
Exceedance probability (%)
90 80 70 60
50 40 30 20 10 0
0
0.5
1 1.5 2 2.5 Normalized damage depth, rd/R
3
3.5
4
Figure 5: Exceedance probability curve for EDZo in limestone with KHh = 1.5 before and after filtering.
Indeed, knowing that no constraint was applied to the sample generation process to ensure that the σmax values are greater than the CI ones, there were some cases among the random samples that produced unrealistic values for rd /R and caused problems on the exceedance probability diagram. To resolve this issue, a filter should be defined to detect and eliminate the samples that did not observe the σmax ≥ CI constraint. For this purpose, the σmax /CI values were calculated for all generated random values and arranged in a vector u. The command u>=1 in MATLAB produces a new binary vector with only 0 and 1 entries so that the element is set to 1 if u ≥ 1 and set to 0 otherwise. The new binary vector was herein named vector v. Now, the vector t = u × v was established by keeping all the elements for which the constraint σmax /CI ≥ 1 was held true. In fact, vector t is the same as vector u but with the problematic elements filtered out. In this way, the set of generated random values was filtered via some simple steps. After filtering the data, the total number of data points (N ) was equal to the number of elements of the vector t. The following commands are given for the filtering process: %% to remove the complex numbers from data u=sigma./CI; v=u>=1; t=u(v); N=length(t); To substitute the random values into the LSF and calculate the random samples, the vector t was used instead of the σmax /CI: LSF=C-1-0.62*(t-1).^0.58; Following this data filtering, the exceedance probability diagram changed. As shown in the blue diagram in Figure 5, significantly higher probability values were obtained upon 9
4
100
3
90 80
2
70
1
60
0
50
-1
40
Probability (%)
Normalized distance from tunnel center, y
filtering, especially for the rd /R < 2. The modified diagram could be accepted as the final solution. As presented in Figure 5, the probability values decreased with increasing the damage depth; so that the probability value fell below 5% for rd /R ≥ 2.28, indicating that the given damage depth is very unlikely to occur. Next, since the exceedance probability corresponding to virtually any normalized radius is available, one could visualize the crack depth and the corresponding probability of the crack penetrating into the rock around the tunnel. This was achieved by adding an outer loop to the Monte Carlo simulation to calculate the probability for each point around the tunnel. The results are shown in Figure 6 in the form of probability contours. As can be observed, the obtained values of probability were very high in the regions near the wall of the tunnel. The main reason for this result is that this region hosted a dense set of interconnected cracks penetrated to this depth. By moving farther from the excavation hole, the probability decreased gradually and eventually reached zero in areas far from the tunnel.
30
-2
20
-3
10
-4 -4 -3 -2 -1 0 1 2 3 4 Normalized distance from tunnel center, x
Figure 6: Contour plot of exceedance probability around the tunnel excavation.
3.4. Exceedance probability curves for EDZi and HDZ In the previous section, the exceedance probability diagram was calculated for the EDZo. The present section addresses EDZi and HDZ. For this purpose, the LSF was modified to represent the damage levels for the new zones using the model developed by Perras and Diederichs [11]. In this model, the coefficients B and D in Eq. 1 are presented for three damage zones with different KHh values (Table 3). The values in this table also help update the LSF for new cases. Performing calculations similar to those in the previous section for cases 1 to 6, new exceedance probability diagrams were obtained for KHh values of 1.5 and 2 (Figures 7a and 7b, respectively). 10
Table 3: Different values of B and D for EDZo, EDZi, and HDZ damage zones in limestone.
Case Rock type KHh 1 1.5 2 1.5 3 1.5 Limestone 4 2 5 2 6 2
100
Exceedance probability curve for case 1 Exceedance probability curve for case 2 Exceedance probability curve for case 3
90 80 70
60 50
40 30
20 10
0
Exceedance probability curve for case 4 Exceedance probability curve for case 5 Exceedance probability curve for case 6
90 Probability of exceedance (%)
Probability of exceedance (%)
100
Zone B D EDZ0 0.66 0.63 EDZi 0.43 0.58 HDZ 0.18 0.34 EDZ0 0.58 0.58 EDZi 0.36 0.49 HDZ 0.12 0.33
80
70 60 50 40 30 20 10
1
1.2
1.4
1.6 1.8 2 2.2 2.4 Normalized damage depth, rd/R
2.6
2.8
3
(a) Cases 1 to 3
0
1
1.2
1.4
1.6 1.8 2 2.2 2.4 Normalized damage depth, rd/R
2.6
2.8
3
(b) Cases 4 to 6
Figure 7: Exceedance probability curve for (a) cases 1 to 3 (KHh = 1.5) and (b) cases 4 to 6 (KHh = 2).
Having obtained the exceedance probability curve, the probability contour around the tunnel could be calculated in each case independently. The results of these calculations are shown in Figure 8. As can be seen, in cases 1 and 4 (i.e., EDZo), the contour covered a wide area around the tunnel, implying that the cracks had penetrated deep into the rock. Subsequently, the area covered by probability contours in cases 2 and 5 (i.e., EDZi) was narrower, while the cases 3 and 6 corresponded to a much limited area around the outer boundary of the tunnel (i.e., HDZ). The HDZ hosted relatively dense yet shallow cracks formed due to severe disturbance of the medium around the excavation region. 4. Discussion 4.1. Reliability sensitivity analysis An important point to address in this study is the effects of σmax and CI on the resulting failure probability. To deal with this issue, a parametric analysis was performed to see how the probability varied with the mean of each variable. For this purpose, the proposed algorithm was repeated for case 1 and the established reliability problem was analyzed by changing the mean of one variable while keeping the other variable constant. As a result, different exceedance probability curves were obtained for different mean values of the considered variable. The results for σmax and CI are shown in Figures 9a and 9b, respectively. 11
4
3
3
3
2
2
2
1
1
1
60
0
0
0
50
-1
-1
-1
-2
-2
-2
-3
-3
-3
-4 -4
-4 -4 -3 -2 -1 0 1 2 3 4 Normalized distance from tunnel center, x (mm)
-4 -4 -3 -2 -1 0 1 2 3 4 Normalized distance from tunnel center, x (mm)
90
80 70
40
Probability (%)
Normalized distance from tunnel center, y
4
30 20
-3 -2 -1 0 1 2 3 Normalized distance from tunnel center, x
4
(a) case 1
(b) case 2
10
(c) case 3 100
4
4
4
3
3
3
2
2
2
1
1
1
60
0
0
0
50
-1
-1
-1
-2
-2
-2
-3
-3
-3
-4 -4
-4 -4 -3 -2 -1 0 1 2 3 4 Normalized distance from tunnel center, x (mm)
-4 -4 -3 -2 -1 0 1 2 3 4 Normalized distance from tunnel center, x (mm)
90 80
70
40
Probability (%)
Normalized distance from tunnel center, y
100
4
30
20
-3 -2 -1 0 1 2 3 Normalized distance from tunnel center, x
(d) case 4
4
(e) case 5
10
(f) case 6
Figure 8: Contour plots of exceedance probability around tunnel excavation for cases 1 to 6.
As can be seen, the probability values increased with σmax , rising the exceedance probability curve to higher levels. On the contrary, the probability dropped sharply with increasing the value of CI, declining the exceedance probability curve. To further examine this issue and determine the extent to which the changes in each parameter affect the crack depth, the difference between the obtained probability curves was calculated. For instance, the difference between exceedance probability (∆P ) corresponding to different maximum tangential stresses was calculated as follows: ∆P = P (σmax ) − P (σmax = 30),
(8)
where P (σmax = 30) is the exceedance probability corresponding to σmax = 30 and P (σmax ) refers to the exceedance probability corresponding to σmax = 50, 70, 90, 110, 130, or 150. The results for σmax and CI are presented in Figures 9c and 9d, respectively. In the following, the position of maximum effect for each variable was estimated by finding the maximum point of each graph and fitting a curve. As can be seen, the σmax affected the cracks with normalized depths ranging between 1.3 and 1.7, while the maximum effect of CI was broader and covered the normalized crack depths ranging between 1.2 and 2.1.
12
100
100
Probability of exceedance (%)
80 70
60
Probability of exceedance (%)
P(σmax=150 MPa) P(σmax=130 MPa) P(σmax=110 MPa) P(σmax=90 MPa) P(σmax=70 MPa) P(σmax=50 MPa) P(σmax=30 MPa)
90
50 40
30 20
80
70 60 50
40 30
20 10
10 0
P(CI=70 MPa) P(CI=60 MPa) P(CI=50 MPa) P(CI=40 MPa) P(CI=30 MPa) P(CI=20 MPa) P(CI=10 MPa)
90
0.5
1
1.5 2 2.5 3 3.5 Normalized damage depth, rd/R
4
0
4.5
0
100 P(σmax=150)-P(σmax=30) P(σmax=130)-P(σmax=30) P(σmax=110)-P(σmax=30) P(σmax=90)-P(σmax=30) P(σmax=70)-P(σmax=30) P(σmax=50)-P(σmax=30) Location of maximum points
90
80 70 60
50 40
30 20 10
0
0.5
1
1.5 2 2.5 3 3.5 Normalized damage depth, rd/R
4
2 3 4 Normalized damage depth, rd/R
5
6
(b) Probability exceedance curve for different values of CI
4.5
(c) Difference between exceedance probability curves corresponding to different values of σmax
Difference between exceedance probabilities (%)
Difference between exceedance probabilities (%)
(a) Probability exceedance curve for different values of σmax
1
100 P(CI=10)-P(CI=70) P(CI=20)-P(CI=70) P(CI=30)-P(CI=70) P(CI=40)-P(CI=70) P(CI=50)-P(CI=70) P(CI=60)-P(CI=70) Location of maximum points
90 80
70 60
50 40 30
20 10 0
0
1
2 3 4 Normalized damage depth, rd/R
5
6
(d) Difference between exceedance probability curves corresponding to different values of CI
Figure 9: Effects of σmax and CI on exceedance probability for case 1.
4.2. Exceedance probability curve for other types of rocks So far, the exceedance probability has been calculated and discussed for the normalized damage depth in limestone. In this section, other types of rocks are examined using the model developed by Perras and Diederichs [11]. Based on the samples reported in this study, the probabilistic characteristics were calculated for two new rock types, as outlined in Table 4. The coefficients B and D (Eq. 1) for granite and mudstone are given in Table 5. Following the substitution of the coefficients B and D into the LSF function (Eq. 3), a new reliability problem was developed. To solve this problem, a new set of random values was generated for σmax and CI according to the rock types, and the Monte Carlo algorithm was utilized to calculate the exceedance probability curves. The results are shown in Figure 10. As can be inferred from Figure 10, the obtained probability for mudstone was larger than that for granite. In fact, the cracks in the mudstone had a greater chance to infiltrate through the rock, while the granite exhibited more resistance against the crack growth. This point was visualized by plotting probability contours around the tunnel (Figure 11). As 13
Table 4: The random variables assumed in the analysis.
Rock type Granite Mudstone
Variable Mean σmax (MPa) 220 CI (MPa) 125 σmax (MPa) 35 CI (MPa) 19
Standard deviation 30 16 11 6
Table 5: Values of B and D in different cases.
Case 7 8 9 10 11 12 13 14 15 16 17 18
Rock type
Granite
Mudstone
KHh 1.5 1.5 1.5 2 2 2 1.5 1.5 1.5 2 2 2
Zone B EDZo 0.62 EDZi 0.41 HDZ 0.09 EDZo 0.58 EDZi 0.36 HDZ 0.11 EDZo 0.71 EDZi 0.49 HDZ 0.2 EDZo 0.66 EDZi 0.39 HDZ 0.15
D 0.58 0.53 0.62 0.65 0.62 0.85 0.59 0.55 0.52 0.59 0.59 0.68
can be seen, in cases 7, 10, 13, and 16 (i.e., EDZo), wide areas were covered by probability contours. In comparison, the areas were drastically limited in cases 8, 11, 14, and 17 (i.e., EDZi) as well as cases 9, 12, 15, and 18 (i.e., HDZ). Additionally, the probability contours in cases 13 to 18 (i.e., mudstone) covered relatively wider areas than cases 7 to 12 (i.e., granite) and showed a slower decreasing trend as one moved away from the outer wall of the tunnel. This observation can be explained by the soft structure of mudstone, as compared to the granite, which allowed deep cracks to grow into the rock. The influence of in-plane stress ratio on the failure probability was investigated based on the resulting exceedance probability curves. For this purpose, the difference in the exceedance probability between different damage zones with different KHh values was calculated as follows: ∆P = P (KHh = 1.5) − P (KHh = 2).
(9)
This gave different graphs for different rock types (Figure 12). In Figure 12, the black, red, and blue curves correspond to the EDZo, EDZi, and HDZ, respectively. The first point to note from Figure 12 is that the effect of KHh was limited in granite, while it was more pronounced in the limestone and maximal in the mudstone. This was, of course, predictable as the hard structure of granite prevents the KHh to highly affect the crack propagation process. The limestone and mudstone were found to be more sensitive to 14
100
Exceedance probability curve for case 7 Exceedance probability curve for case 8 Exceedance probability curve for case 9
90 80
Probability of exceedance (%)
Probability of exceedance (%)
100
70 60 50 40 30 20 10 1
1.1
1.2
100 Probability of exceedance (%)
1.3 1.4 1.5 1.6 1.7 1.8 Normalized damage depth, rd/R (a) Cases 7 to 9 ( K Hh = 1.5)
1.9
70 60 50 40 30 20
0
2
80
1
1.1
1.2
100
Exceedance probability curve for case 13 Exceedance probability curve for case 14 Exceedance probability curve for case 15
90
70
60 50
40 30
20 10
0
80
10
Probability of exceedance (%)
0
Exceedance probability curve for case 10 Exceedance probability curve for case 11 Exceedance probability curve for case 12
90
1.3 1.4 1.5 1.6 1.7 Normalized damage depth, rd/R (b) Cases 10 to 12 (K Hh = 2)
1.8
1.9
2
Exceedance probability curve for case 16 Exceedance probability curve for case 17 Exceedance probability curve for case 18
90
80 70
60 50
40 30
20 10
1
1.5
2 2.5 Normalized damage depth, rd/R (c) Cases 13 to 15 (K Hh = 1.5)
3
3.5
01
1.5
2 2.5 Normalized damage depth, rd/R (d) Cases 16 to 18 (K Hh = 2)
3
3.5
Figure 10: Exceedance probability curves for (a) cases 7 to 9, (b) cases 10 to 12, (c) cases 13 to 15, and (d) cases 16 to 18.
KHh because of their softer structures. As another important point, the maximum impact of KHh on the EDZo was found to occur at large normalized crack depths approximately ranging between 1.5 and 2, while the maximum impacts of the KHh in the EDZi and HDZ occurred at the normalized crack depths ranging within 1.3-1.5 and 1-1.3, respectively. This finding seems reasonable since the EDZo inherently covers the deep cracks rooted in the rock, while the EDZo and HDZ tend to cover shallow yet dense crack systems. The probability differences in Figure 12 also show that the impact of KHh on the HDZ is much greater than those on the EDZi and then the EDZo. This result is understandable because the earliest impact of a change in the in-plane stress is imposed on the neighboring medium of the excavation hole, followed by a decrease in the impact with distance from the excavation wall. However, despite the general sense that an increase in the KHh can reduce the failure probability, the case observed in granite shows that the probability of failure can even increase with increasing the KHh , making the difference presented in Figure 12b negative. This issue might be related to the coefficients reported in Perras and Diederichs model [11]. There should be a particular factor affecting the regression models of B and D so that the Eq. 1 reports higher damage when KHh increases. Therefore, when using the 15
4
3
3
3
2
2
2
1
1
1
60
0
0
0
50
-1
-1
-1
-2
-2
-2
-3
-3
-3
-4 -4
-4 -4 -3 -2 -1 0 1 2 3 4 Normalized distance from tunnel center, x (mm)
-4 -4 -3 -2 -1 0 1 2 3 4 Normalized distance from tunnel center, x (mm)
80 70
20
-3 -2 -1 0 1 2 3 Normalized distance from tunnel center, x
4
(b) Case 8
10
(c) Case 9 100
4
4
4
3
3
3
2
2
2
1
1
1
60
0
0
0
50
-1
-1
-1
-2
-2
-2
-3
-3
-3
-4 -4
-4 -4 -3 -2 -1 0 1 2 3 4 Normalized distance from tunnel center, x (mm)
-4 -4 -3 -2 -1 0 1 2 3 4 Normalized distance from tunnel center, x (mm)
90 80 70
40
Probability (%)
Normalized distance from tunnel center, y
40 30
30 20
-3 -2 -1 0 1 2 3 Normalized distance from tunnel center, x
4
(d) Case 10
(e) Case 11
10
(f) Case 12 100
4
4
4
3
3
3
2
2
2
1
1
1
60
0
0
0
50
-1
-1
-1
-2
-2
-2
-3
-3
-3
-4 -4
-4 -4 -3 -2 -1 0 1 2 3 4 Normalized distance from tunnel center, x (mm)
-4 -4 -3 -2 -1 0 1 2 3 4 Normalized distance from tunnel center, x (mm)
90 80 70
40
Probability (%)
Normalized distance from tunnel center, y
90
Probability (%)
Normalized distance from tunnel center, y
4
(a) Case 7
30 20
-3 -2 -1 0 1 2 3 Normalized distance from tunnel center, x
4
(g) Case 13
(h) Case 14
10
(i) case 15 100
4
4
4
3
3
3
2
2
2
1
1
1
60
0
0
0
50
-1
-1
-1
-2
-2
-2
-3
-3
-3
-4 -4
-4 -4 -3 -2 -1 0 1 2 3 4 Normalized distance from tunnel center, x (mm)
-4 -4 -3 -2 -1 0 1 2 3 4 Normalized distance from tunnel center, x (mm)
90
80 70
40 30 20
-3 -2 -1 0 1 2 3 Normalized distance from tunnel center, x
(j) Case 16
4
(k) Case 17
10
(l) Case 18
Figure 11: Contour plots of exceedance probability around the excavation for cases 7 to 18.
16
Probability (%)
Normalized distance from tunnel center, y
100
4
Difference between exceedance probabilities (%)
Difference between exceedance probabilities (%)
60 Difference between cases 1 and 4 Difference between cases 2 and 5 Difference between cases 3 and 6
50 40
30 20 10 0 -10
0
0.5
1
1.5 2 2.5 3 Normalized damage depth, rd/R
3.5
4
25 Difference between cases 7 and 10 Difference between cases 8 and 11 Difference between cases 9 and 12
20
15 10
5 0 -5
-10 -15
-20 -25
0
0.5
1 1.5 Normalized damage depth, rd/R
2
(b) Granite Difference between exceedance probabilities (%)
(a) Limestone 25 Difference between cases 13 and 16 Difference between cases 14 and 17 Difference between cases 15 and 18
20
15 10 5
0 -5
0
0.5
1
1.5 2 2.5 3 3.5 Normalized damage depth, rd/R
4
4.5
5
(c) Mudstone
Figure 12: Effects of in-plane stress ratio in (a) limestone, (b) granite, and (c) mudstone.
Perras and Diederichs model [11] in granite for estimating HDZ zone, further caution should be considered, with the coefficients B and D to be calibrated by field tests at the specific site under investigation wherever possible. 5. Conclusion In this study, excavation-induced failure probability was investigated in three different types of rock including granite, mudstone, and limestone. For this purpose, a probabilistic model was established based on the model proposed by Perras and Diederichs. Then, the Monte Carlo sampling method was adopted to approximate the exceedance probability for different induced crack depths, with the results presented in the form of exceedance probability curves. Effects of the maximum tangential stress and the crack initiation threshold on failure probability were also examined using reliability sensitivity analysis. The main conclusions drawn from this study are as follows: As the depth of the crack increases, the probability falls sharply; as a result, the exceedance probability is less than 5% for the normalized crack depths of more than 2.28.
17
The effect of σmax is maximal at normalized crack depths ranging between 1.5 and 1.7, while that of CI is maximal when the normalized crack depth falls between 1.2 and 2.1. A change of in-plane stress ratio shows the most significant impacts on the EDZo when the normalized crack depth ranges between 1.5 and 2, while its maximum effects on the EDZi and HDZ is estimated to occur at the normalized crack depths ranging within 1.3-1.5 and 1-1.3, respectively. Increasing the in-plane stress ratio reduces the failure probability and causes the exceedance probability curve to locate at lower positions.
It is noteworthy that the analyses presented in this study are based on the assumption that the involved variables are independent. Estimating the failure probability around the excavation hole considering possible correlations between the involved variables is an interesting subject that will be investigated by the authors in their future research. References [1] L. Ding, L. Ma, H. Luo, M. Yu, X. Wu, Wavelet analysis for tunneling-induced ground settlement based on a stochastic model, Tunnelling and Underground Space Technology 26 (2011) 619–628. [2] Y. Fan, W. Lu, Y. Zhou, P. Yan, Z. Leng, M. Chen, Influence of tunneling methods on the strainburst characteristics during the excavation of deep rock masses, Engineering Geology 201 (2016) 85–95. [3] K. Zhang, J. L. Chavez Torres, Z. Zang, Numerical analysis of pipelines settlement induced by tunneling, Advances in Civil Engineering 2019 (2019) 1–10. [4] R. Read, 20 years of excavation response studies at aecl’s underground research laboratory, International Journal of Rock Mechanics and Mining Sciences 41 (2004) 1251 – 1275. Rock Mechanics Results from the Underground Research Laboratory, Canada. [5] M. Tao, Z. Hong, K. Peng, P. Sun, M. Cao, K. Du, Evaluation of excavation-damaged zone around underground tunnels by theoretical calculation and field test methods, Energies 12 (2019) 1–18. [6] L. X. Xie, Q. B. Zhang, J. C. Gu, W. B. Lu, S. Q. Yang, H. W. Jing, Z. L. Wang, Damage evolution mechanism in production blasting excavation under different stress fields, Simulation Modelling Practice and Theory 97 (2019) 101969. [7] S. Wang, X. Li, J. Yao, F. Gong, X. Li, K. Du, M. Tao, L. Huang, S. Du, Experimental investigation of rock breakage by a conical pick and its application to non-explosive mechanized mining in deep hard rock, International Journal of Rock Mechanics and Mining Sciences 122 (2019) 104063. 18
[8] H. Wang, Y. Jiang, S. Xue, B. Shen, C. Wang, J. Lv, T. Yang, Assessment of excavation damaged zone around roadways under dynamic pressure induced by an active mining process, International Journal of Rock Mechanics and Mining Sciences 77 (2015) 265– 277. [9] G. Walton, M. Lato, H. Ansch¨ utz, M. A. Perras, M. S. Diederichs, Non-invasive detection of fractures, fracture zones, and rock damage in a hard rock excavation — ¨ o Hard Rock Laboratory in Sweden, Engineering Geology 196 Experience from the Asp¨ (2015) 210–221. [10] T. Siren, P. Kantia, M. Rinne, Considerations and observations of stress-induced and construction-induced excavation damage zone in crystalline rock, International Journal of Rock Mechanics and Mining Sciences 73 (2015) 165–174. [11] M. A. Perras, M. S. Diederichs, Predicting excavation damage zone depths in brittle rocks, Journal of Rock Mechanics and Geotechnical Engineering 8 (2016) 60–74. [12] M. A. Perras, M. S. Diederichs, A review of excavation damage zones in sedimentary rocks with emphasis on numerical modelling for edz definition, The Canadian Geotechnical Society’s Annual Conference - Geo2010, Alberta, Canada. [13] E. Ghazvinian, Fracture initiation and propagation in low porosity crystalline rocks: implications for excavation damage zone (EDZ) mechanics, Phd thesis, Queen’s University, 2015. [14] S. Kwon, C. S. Lee, S. J. Cho, S. W. Jeon, W. J. Cho, An investigation of the excavation damaged zone at the KAERI underground research tunnel, Tunnelling and Underground Space Technology 24 (2009) 1–13. [15] W.-B. Lu, Y.-G. Hu, J.-H. Yang, M. Chen, P. Yan, Spatial distribution of excavation induced damage zone of high rock slope, International Journal of Rock Mechanics and Mining Sciences 64 (2013) 181–191. [16] W.-J. Cho, J.-S. Kim, C. Lee, H.-J. Choi, Gas permeability in the excavation damaged zone at KURT, Engineering Geology 164 (2013) 222–229. [17] H. Q. Yang, Y. Y. Zeng, Y. F. Lan, X. P. Zhou, Analysis of the excavation damaged zone around a tunnel accounting for geostress and unloading, International Journal of Rock Mechanics and Mining Sciences 69 (2014) 59–66. [18] W. C. Zhu, J. Wei, J. Zhao, L. L. Niu, 2D numerical simulation on excavation damaged zone induced by dynamic stress redistribution, Tunnelling and Underground Space Technology 43 (2014) 315–326. [19] Z. Wu, Y. Jiang, Q. Liu, H. Ma, Investigation of the excavation damaged zone around deep TBM tunnel using a Voronoi-element based explicit numerical manifold method, International Journal of Rock Mechanics and Mining Sciences 112 (2018) 158–170. 19
[20] Y. Xue, F. Gao, X. Liu, X. Liang, Permeability and pressure distribution characteristics of the roadway surrounding rock in the damaged zone of an excavation, International Journal of Mining Science and Technology 27 (2017) 211–219. [21] A. A. Rodriguez, U. Kuhlmann, P. Marschall, 3D modelling of the Excavation Damaged Zone using a Marked Point Process technique, Geomechanics for Energy and the Environment 17 (2019) 29–46. [22] C. D. Martini, R. S. Read, J. B. Martino, Observations of brittle failure around a circular test tunnel, International Journal of Rock Mechanics and Mining Sciences 34 (1997) 1065–1073. [23] A. Fakhimi, F. Carvalho, T. Ishida, J. F. Labuz, Simulation of failure around a circular opening in rock, International Journal of Rock Mechanics and Mining Sciences 39 (2002) 507–515. [24] E. Hrubesova, L. Duris, Assessment of tunnel’s face support pressure, IOP Conference Series: Materials Science and Engineering 236 (2017) 1–7. [25] N. Dadashzadeh, M. S. Diederichs, Reliability of prediction for tunnel excavation damage zone depth in brittle rocks, in: ISRM International Symposium - 10th Asian Rock Mechanics Symposium, International Society for Rock Mechanics and Rock Engineering, Singapore, 2018. [26] B. Tang, H. Cheng, Y. Tang, Z. Yao, M. Li, T. Zheng, J. Lin, Excavation damaged zone depths prediction for TBM-excavated roadways in deep collieries, Environmental Earth Sciences 77 (2018) 1–14. [27] I. Vazaios, N. Vlachopoulos, M. Diederichs, Assessing fracturing mechanisms and evolution of excavation damaged zone of tunnels in interlocked rock masses at high stresses using a finite-discrete element approach, Journal of Rock Mechanics and Geotechnical Engineering 11 (2019) 701–722. [28] I. Vazaios, N. Vlachopoulos, M. S. Diederichs, Mechanical analysis and interpretation of excavation damage zone formation around deep tunnels within massive rock masses using hybrid finite–discrete element approach: case of Atomic Energy of Canada Limited (AECL) Underground Research Laboratory (URL) test tunnel, Canadian Geotechnical Journal 56 (2019) 35–59. [29] Q. Xie, K. Peng, Space-time distribution laws of tunnel excavation damaged zones (EDZs) in deep mines and EDZ prediction modeling by random forest regression, Advances in Civil Engineering 2019 (2019) 1–13. [30] J. C. Kok-Kwang Phoon, Risk and reliability in geotechnical engineering, CRC Press, Florida, 2014.
20
[31] G. Baecher, J. Christian, Reliability and statistics in geotechnical engineering, Wiley, 2003. [32] M. Shadab Far, Y. Wang, Approximation of the Monte Carlo sampling method for reliability analysis of structures, Mathematical Problems in Engineering 2016 (2016). [33] D. P. Kroese, T. Taimre, Z. I. Botev, Handbook of Monte Carlo methods, Wiley, 2011.
21
Highlights - The depth of excavation-induced damage zones is established as a reliability problem. - The exceedance probability curve is developed for three rock types, including limestone, granite, and mudstone. - Failure probability contours are presented around the excavation region. - The effects of maximum tangential stress and the crack initiation threshold are investigated.
Declaration of interests ☒ The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper. ☐The authors declare the following financial interests/personal relationships which may be considered as potential competing interests:
Declarations of interest: none
S2590-1230(19)30075-1
DOI:
https://doi.org/10.1016/j.rineng.2019.100075
Reference:
RINENG 100075
To appear in:
Results in Engineering
Received Date: 10 November 2019 Accepted Date: 20 November 2019
Please cite this article as: M. Shadabfar, H. Huang, H. Kordestani, E.V. Muho, Probabilistic Modeling of Excavation-Induced Damage Depth around Rock-Excavated Tunnels, Results in Engineering, https:// doi.org/10.1016/j.rineng.2019.100075. This is a PDF file of an article that has undergone enhancements after acceptance, such as the addition of a cover page and metadata, and formatting for readability, but it is not yet the definitive version of record. This version will undergo additional copyediting, typesetting and review before it is published in its final form, but we are providing this version to give early visibility of the article. Please note that, during the production process, errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain. © 2019 The Author(s). Published by Elsevier B.V.
Probabilistic Modeling of Excavation-Induced Damage Depth around Rock-Excavated Tunnels Mahdi Shadab Fara,∗, Hadi Kordestanib , Edmond V. Muhoc a
Department of Geotechnical Engineering, Tongji University, Shanghai, China Department of Civil Engineering, Qingdao University of Technology, Qingdao, China c Department of Disaster Mitigation for Structures, College of Civil Engineering, Tongji University, Shanghai, China b
Abstract The cracked zones caused by tunnel excavations are among the major reasons for structural instability in mechanized tunnel projects. Thus, various methods have been proposed to evaluate such cracks in terms of size and depth. However, deterministic estimations have been associated with a margin of error due to the inherent uncertainties of the respective problems as well as the approximation nature of the existing methods. Therefore, in this study, a probabilistic approach was adopted to assess the depth of a cracked zone induced around a tunnel. For this purpose, the problem was formulated as a reliability model, and the Monte Carlo method was utilized to calculate the exceedance probability for the crack depth around the tunnel. The exceedance probability was then calculated for any different depths of crack and presented as probability contours around the excavation zone. It was demonstrated that the exceedance probability is reduced sharply with an increase in the crack depth such that it fell below 5% for the cracks deeper than 2.28 times the tunnel radius. Effects of the maximum tangential stress at excavation boundary and crack initiation threshold on the exceedance probability were further calculated using the reliability sensitivity analysis. Eventually, the influence of in-plain stress ratio on the failure probability was investigated, and the results were reported for three rock types, including limestone, granite, and mudstone. Keywords: Excavation-induced damage (EDZ), uncertainty modelling, Monte Carlo sampling method, failure probability, reliability sensitivity analysis. 1. Introduction In the course of excavation, redistribution of the in-situ stresses in the surrounding rock mass leads to a gradual increase in the tangential stresses while decreasing the radial stresses [1–4]. This develops an extreme state at the opening surface, causing a stress concentration ∗
Corresponding author Email addresses: [email protected] (Mahdi Shadab Far ), [email protected] (Hadi Kordestani), [email protected] (Edmond V. Muho) Preprint submitted to Results in Engineering
November 10, 2019
around the excavation area [5]. When the resulting stress exceeds the failure strength of the rock mass, excavation-induced cracks are initiated and gradually penetrate the rock mass until a new stress equilibrium is reached [6, 7]. Ultimately, these single cracks get interconnected to form damage layers with different crack densities. The first damage layer, which serves as host for very dense and interconnected cracks, is called the high damage zone (HDZ) [8]. Moving outward, layers of lower crack density but connected crack systems become dominant. These layers are collectively referred to as the inner excavation damage zone (EDZi). Then, the cracks diverge into joints in an outer excavation damage zone (EDZo) [9]. Beyond the EDZs, there is an elastic and damage-free yet affected zone called the excavation disturbed zone (EdZ) or the excavation influence zone (EIZ). The EIZ is preferred to avoid possible confusions between the EDZ and the EdZ [10, 11]. The described zones are schematically shown in Figure 1. EIZ EDZo
EDZi Excavation area
HDZ
Figure 1: Different affected zones around an excavation area.
The generated damage zones threaten the stability and performance of underground structures and are reportedly the source of major geological hazards, such as spalling, desaturation, and roof collapse [12, 13]. Accordingly, numerous studies have been conducted to investigate the initiation, propagation, extension, dimensions, and induced permeability of the EDZs [14–21]. In this respect, Martini et al. [22] examined the formation and propagation of the EdZ by recording data around a large-scale underground opening. Fakhimi et al. [23] performed biaxial compression tests on sandstone specimens to simulate the failure around a circular opening in brittle rocks. Similarly, Perras and Diederichs [11] used damage initiation and spalling limit (DISL) and proposed a guideline to predict the depths of EDZs in brittle rocks. In line with the above-mentioned models, several other empirical methods and numerical models are available in the literature for addressing the excavation-induced damage. Considering the assumptions of the Hoek-Brown’s Failure Criterion and the conceptual model of Diederichs for rock behavior, Hrubesova and Duris [24] examined three damage zones, namely HDZ, EDZ, and EIZ, around an explosion site and applied the same methodology to a real case study. Dadashzadeh and Diederichs [25] employed cohesion weakening 2
and frictional strengthening model and developed a finite-difference method in FLAC 2D software to predict variations of the EDZ depth in brittle rocks. Tang et al. [26] presented a modified failure criterion using a series of cyclic loading tests to formulate an analytical solution for predicting different EDZs. Combining the discrete fracture networks and a twodimensional finite-discrete element method (FDEM) code, Vazaios et al. [27] investigated the effects of joints on the reduction of the rock strength and bearing capacity to accommodate damage zones. In another study [28], they adopted the FDEM to investigate rock mass behavior and examine the spalling-induced damage extent. Xie and Peng [29] used an ultrasonic velocity detector to assess the range of EDZs in a rock medium during a drilling process and evaluated the damage extension from the drilling data. Tao et al. [5] developed an analytical model for estimating the extent of EDZ and provided a closed-form solution for the established problem. Then, they conducted a parametric study to investigate the effects of different variables on the induced damage. The results revealed that the geological strength index and uniaxial compressive strength imposed the largest contributions to the damage rates. However, possible disturbances (e.g., mechanical drilling and the waves generated in different stages of the drilling operation) could significantly affect the problem, even dominating the effects of the parameters involved in mathematical modeling in some cases. To the authors’ best knowledge, most of the studies on the excavation-induced damage size have adopted a deterministic approach for formulating the problem. However, the uncertainties associated with the stress conditions and crack resistance of the rock are so influential that cannot be neglected. In fact, deterministic modeling is equivalent to assuming 100% certainty about the load and bearing capacity of the rock medium, which makes no sense when it comes to real cases. Therefore, an alternative method is required to take the uncertainty into consideration and express the problem probabilistically. In the present study, reliability modeling is used to address this issue. To this end, the input parameters are initially modeled as random variables with certain means and standard deviations. Then, a limit state function is constructed based on the model developed by Perras and Diederichs [11] and the depths of the EDZs are expressed as a reliability problem. The Monte Carlo sampling method is then utilized to solve the established problem and the results are reported as exceedance probability at different depths of the EDZs and presented as probability contour maps around the tunnel. Effects of different parameters on the failure probability are also investigated through reliability sensitivity analysis. Finally, the effect of in-plane stress ratio is demonstrated by exceedance probability curves and the results are reported for different rock types. 2. Modeling As previously outlined in the Introduction section, this study attempts to express the formation of excavation-induced damage zone as a probabilistic problem, because deterministic models are mainly limited to specific solutions even for real cases encountering high levels of uncertainty. More specifically, a deterministic model assumes 100% probability for the occurrence of a particular damage depth, which makes it highly unrealistic. A more ra3
tional solution to this problem is to calculate exceedance probabilities for different probable damage depths; i.e., any damage depth corresponds to a particular value of exceedance probability. For this purpose, one must begin with defining a deterministic model to calculate the crack depth via an explicit relationship. Then, the deterministic model can be converted to a probabilistic one by setting the input parameters to random variables, constructing a reliability problem. Later, a conventional reliability analysis method can be adopted to solve the established problem (here, we use the Monte Carlo sampling method for this purpose). To this end, the Perras and Diederichs model [11] was considered as the primary deterministic model. In the Perras and Diederichs study [11], the natural variations of brittle rock were considered explicitly by incorporating a laboratory model into numerical simulations with brittle constitutive behavior. The results of all numerical modeling (different rock types, stress conditions, etc.) were then collected, and the best average fit for each damage zone was developed using regression models. As a result, a simple but reasonable estimation of normalized damage depth in three different damage classes (i.e., EDZi, EDZo, and HDZ) was provided covering a wide range of possible scenarios: σmax rd = 1 + B( − 1)D , (1) R CI where rd represents the depth of inner EDZ (EDZi), outer EDZ (EDZo), or highly damaged zone (HDZ), R denotes the tunnel radius, σmax presents the maximum tangential stress at excavation boundary, and CI is regarded as the crack initiation threshold. B and D are constant regression parameters that are available for different rock types, in-plane stress ratios (KHh ), and types of damaged zones (i.e., EDZi, EDZo, or HDZ). For instance, the values of B and D are 0.66 and 0.63, respectively, when predicting the depth of EDZo in limestone with KHh = 1.5. Denoting the excavation influence on the rock as S and the bearing capacity of the rock as C, the limit state function (LSF ) can be expressed as Eq. (2). If the disturbance caused by the excavation overshadows the rock bearing capacity (i.e., C is less than S), then the model faces damage. Otherwise, if C is greater than S, then the model is assumed to be safe. LSF = C − S.
(2)
When the rock is highly resistant to excavation-induced disturbance, the depth of cracks is very limited. However, if the rock is too weak to withstand such disturbances, the generated cracks penetrate more into the rock. As a result, the crack depth can well represent the rock bearing capacity. Hence, assuming that S is the estimated normalized depth of damage generated around the excavation hole (e.g., the value obtained from the Perras and Diederichs model [11]) and C is a specified normalized depth of crack (e.g., rd /R = 21 ), which is taken as the boundary of the model failure, the LSF is rewritten as follows: 1
The other possible values will be considered in Section 3.2.
4
Table 1: Random variables assumed in the analysis.
Rock type
Variable Mean σmax (MPa) 72 Limestone CI (MPa) 41
Standard deviation 21 12
σmax − 1)D . (3) CI Here, a major objective is to determine how much probability is needed to face damages with depths beyond C. This is equivalent to the probability of S > C, which can be calculated as follows: LSF = C − 1 − B(
Failure Probability = P (LSF < 0) = P (S > C).
(4)
Various methods have been proposed to estimate this probability [30, 31]. Considering its advantages [32, 33], the Monte Carlo sampling method was used in the present paper. To use the Monte Carlo method, sufficient numbers of random values should be generated for all variables involved in the problem. The random values should correspond to the mean and standard deviation of the respective variable and follow the governing distribution function. By substituting the generated random values into the LSF , a so-called LSF random sample was achieved for each set of random values. Counting the number of less-than-zero samples and dividing it by the total number of samples, the desired exceedance probability was calculated. To identify the optimal number of the random values (sample size), the sample size was increased gradually and the resultant changes in the exceedance probability were evaluated. This was continued until the probability was no longer sensitive to the samples size and the number of samples was high enough to ensure the convergence. Then, the calculated failure probability was recorded as the output of the Monte Carlo method. The process of using the Monte Carlo method, from selecting the primary model to the calculation of the exceedance probability, is schematically shown in Figure 2. 3. Analysis 3.1. Reliability analysis To start using the Monte Carlo sampling method, the input parameters, including the maximum tangential stress and crack initiation threshold, were defined as random variables. In this way, uncertainties are incorporated into the model parameters by considering a relatively wide range of possible caseloads in the analysis process. Here, we assumed that the variables followed normal distribution functions with means and standard deviations presented in Table 1. This assumption is in line with the values presented in the Perras and Diederichs model [11]. Continuing with the algorithm, 20,000 uniformly-distributed random values were generated for each of the considered variables. Then, using the inverse CDF of the normal distribution function (Eq. 5), the generated random values were transferred from uniform 5
Establish a primary deterministic model
Define the random variables
Define the limit state function (LSF)
Calculate the LSF based on the generated random numbers
Calculate the exceedance probability
NO
Monte Carlo loop
Increase the number of random samples
Generate random numbers for each random variable
Is the analysis converged? YES Record the probability, Pf
Figure 2: Schematic view of the Monte Carlo algorithm.
distribution to normal distribution with the probabilistic characteristics specified in Table 1. Histograms of the generated random values are shown in Figure 3 (blue column bars). zi = Φ−1 (ui ).
(5)
Using the generated random values, the rd /R ratio was calculated based on Eq. 1 for limestone with KHh = 1.5, and the proposed LSF was estimated through Eq. 3. The exceedance probability was then calculated via the following formula: n , (6) N where n is the number of cases for which LSF < 0 and N represents the total number of random values. In the current analysis, there were 2,354 less-than-zero cases among the whole random samples. Therefore, the exceedance probability was calculated as 2354/20000 = 0.1177 = 11.77%. To verify the results and ensure the convergence of the analysis, sensitivity of the results to the sample values was investigated. For this purpose, exceedance probability was calculated for different sample sizes and plotted on a graph. The results are illustrated in Figure 4. As can be seen in Figure 4, the results highly fluctuate in smaller sample sizes, so that P =
6
15
Histogram of σmax (MPa) Cases where LSF<0 Fitted distribution
Relative frequency (%)
Relative frequency (%)
15
10
5
00
10
5
0
20 40 60 80 100 120 140 160 Maximum tangential stress at the excavation boundary, σmax (MPa)
Histogram of CI (MPa) Cases where LSF<0 Fitted distribution
0
10
20
30 40 50 60 70 Crack initiation threshold, CI (MPa)
(a)
80
90
100
(b)
Probability of exceedance (%) for rd/R=2
Figure 3: Histograms of the generated random values for (a) σmax and (b) CI.
24 22
Probability of exceedance Final value
20 18 16 14 12 10 8
60
0.2
0.4
0.6
0.8 1 1.2 Sample size
1.4
1.6
1.8
2 × 104
Figure 4: Exceedance probability versus sample size.
a slight change in the sample size could make a considerable change in the exceedance probability. However, as the sample size increased, the change in the exceedance probability became smaller and the probability converged to 11.77% with negligible variations for sample sizes beyond 10,000. 3.2. Exceedance probability curve Thus far, the exceedance probability was calculated for C = rd /R = 2. However, the purpose of the present study is to calculate the exceedance probability for any desired value of normalized damage depth. Therefore, the rd /R was defined as a decision variable rather than a constant value. Then, assuming different values for the rd /R, the established problem was solved and the exceedance probability was calculated accordingly. The results are shown in Table 2 for different values of rd /R. The exceedance probability was plotted versus different values of rd /R and a curve was fitted to the results (black line in Figure 5). This curve, known as the exceedance probability 7
Table 2: Exceedance probability for different decision variables.
rd /R 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 2 2.2 2.4 2.6 2.8 3
Less-than-zero number P = n/N 19996 0.9998 19996 0.9998 19996 0.9998 19921 0.9961 17997 0.8999 16730 0.8365 13536 0.6768 8881 0.4441 4824 0.2412 2354 0.1177 1154 0.0577 592 0.0296 325 0.0163 206 0.0103 132 0.0066
P (%) rd /R 99.98 3.2 99.98 3.4 99.98 3.6 99.61 3.8 89.99 4 83.65 4.2 67.68 4.4 44.41 4.6 24.12 4.8 11.77 5 5.77 5.2 2.96 5.4 1.63 5.6 1.03 5.8 0.66 6
Less-than-zero number P = n/N 87 0.0044 65 0.0033 54 0.0027 46 0.0023 35 0.0018 25 0.0013 21 0.0011 19 0.0010 16 0.0008 15 0.0008 11 0.0006 9 0.0005 8 0.0004 8 0.0004 8 0.0004
P (%) 0.44 0.33 0.27 0.23 0.18 0.13 0.11 0.10 0.08 0.08 0.06 0.05 0.04 0.04 0.04
curve, shows the probability of exceedance for different values of normalized damage depth. The obtained exceedance probability curve has practical applications in the reliability-based evaluation of real projects where a target probability level (allowable probability) is set and the corresponding crack depth is determined. A typical numerical analysis is then required to establish the crack dimensions in the real project. If the size of the crack is greater than the allowable value, the expected probability is violated and, hence, further measures must be taken to reduce failure probability. Otherwise, the project is safe and the implementation process can be continued. Moreover, this graph shows the rock capacity to withstand the disturbances caused by excavation. If the diagram is generally high, the cracks are more likely to grow deep into the rock and, consequently, the failure probability is high. On the contrary, if the diagram is generally low, the cracks cannot easily grow through the rock and, thus, the rock medium exhibits a larger capacity to accommodate the excavation effects. 3.3. Data filtering In the previous section, the convergence of exceedance probability was ensured and it was demonstrated that the given probability was not sensitive to the number of samples. Nevertheless, a brief overview of the exceedance probability diagram presented in Figure 5 by black line reveals an error in this graph. It is expected that the probability value is 100% for rd /R ≤ 1 and the exceedance probability diagram begins to decline from rd /R = 1; however, this trend is not observed on the diagram. Looking at the definition of LSF in Eq. 3, it can be inferred that nonpositive values of (σmax /CI − 1) lead to the introduction of non-real values into the random samples, thereby affecting the calculated final probability value. For this reason, it should be ensured that: σmax − 1 ≥ 0 ⇒ σmax ≥ CI. CI 8
(7)
100 Before filtering After filtering
Exceedance probability (%)
90 80 70 60
50 40 30 20 10 0
0
0.5
1 1.5 2 2.5 Normalized damage depth, rd/R
3
3.5
4
Figure 5: Exceedance probability curve for EDZo in limestone with KHh = 1.5 before and after filtering.
Indeed, knowing that no constraint was applied to the sample generation process to ensure that the σmax values are greater than the CI ones, there were some cases among the random samples that produced unrealistic values for rd /R and caused problems on the exceedance probability diagram. To resolve this issue, a filter should be defined to detect and eliminate the samples that did not observe the σmax ≥ CI constraint. For this purpose, the σmax /CI values were calculated for all generated random values and arranged in a vector u. The command u>=1 in MATLAB produces a new binary vector with only 0 and 1 entries so that the element is set to 1 if u ≥ 1 and set to 0 otherwise. The new binary vector was herein named vector v. Now, the vector t = u × v was established by keeping all the elements for which the constraint σmax /CI ≥ 1 was held true. In fact, vector t is the same as vector u but with the problematic elements filtered out. In this way, the set of generated random values was filtered via some simple steps. After filtering the data, the total number of data points (N ) was equal to the number of elements of the vector t. The following commands are given for the filtering process: %% to remove the complex numbers from data u=sigma./CI; v=u>=1; t=u(v); N=length(t); To substitute the random values into the LSF and calculate the random samples, the vector t was used instead of the σmax /CI: LSF=C-1-0.62*(t-1).^0.58; Following this data filtering, the exceedance probability diagram changed. As shown in the blue diagram in Figure 5, significantly higher probability values were obtained upon 9
4
100
3
90 80
2
70
1
60
0
50
-1
40
Probability (%)
Normalized distance from tunnel center, y
filtering, especially for the rd /R < 2. The modified diagram could be accepted as the final solution. As presented in Figure 5, the probability values decreased with increasing the damage depth; so that the probability value fell below 5% for rd /R ≥ 2.28, indicating that the given damage depth is very unlikely to occur. Next, since the exceedance probability corresponding to virtually any normalized radius is available, one could visualize the crack depth and the corresponding probability of the crack penetrating into the rock around the tunnel. This was achieved by adding an outer loop to the Monte Carlo simulation to calculate the probability for each point around the tunnel. The results are shown in Figure 6 in the form of probability contours. As can be observed, the obtained values of probability were very high in the regions near the wall of the tunnel. The main reason for this result is that this region hosted a dense set of interconnected cracks penetrated to this depth. By moving farther from the excavation hole, the probability decreased gradually and eventually reached zero in areas far from the tunnel.
30
-2
20
-3
10
-4 -4 -3 -2 -1 0 1 2 3 4 Normalized distance from tunnel center, x
Figure 6: Contour plot of exceedance probability around the tunnel excavation.
3.4. Exceedance probability curves for EDZi and HDZ In the previous section, the exceedance probability diagram was calculated for the EDZo. The present section addresses EDZi and HDZ. For this purpose, the LSF was modified to represent the damage levels for the new zones using the model developed by Perras and Diederichs [11]. In this model, the coefficients B and D in Eq. 1 are presented for three damage zones with different KHh values (Table 3). The values in this table also help update the LSF for new cases. Performing calculations similar to those in the previous section for cases 1 to 6, new exceedance probability diagrams were obtained for KHh values of 1.5 and 2 (Figures 7a and 7b, respectively). 10
Table 3: Different values of B and D for EDZo, EDZi, and HDZ damage zones in limestone.
Case Rock type KHh 1 1.5 2 1.5 3 1.5 Limestone 4 2 5 2 6 2
100
Exceedance probability curve for case 1 Exceedance probability curve for case 2 Exceedance probability curve for case 3
90 80 70
60 50
40 30
20 10
0
Exceedance probability curve for case 4 Exceedance probability curve for case 5 Exceedance probability curve for case 6
90 Probability of exceedance (%)
Probability of exceedance (%)
100
Zone B D EDZ0 0.66 0.63 EDZi 0.43 0.58 HDZ 0.18 0.34 EDZ0 0.58 0.58 EDZi 0.36 0.49 HDZ 0.12 0.33
80
70 60 50 40 30 20 10
1
1.2
1.4
1.6 1.8 2 2.2 2.4 Normalized damage depth, rd/R
2.6
2.8
3
(a) Cases 1 to 3
0
1
1.2
1.4
1.6 1.8 2 2.2 2.4 Normalized damage depth, rd/R
2.6
2.8
3
(b) Cases 4 to 6
Figure 7: Exceedance probability curve for (a) cases 1 to 3 (KHh = 1.5) and (b) cases 4 to 6 (KHh = 2).
Having obtained the exceedance probability curve, the probability contour around the tunnel could be calculated in each case independently. The results of these calculations are shown in Figure 8. As can be seen, in cases 1 and 4 (i.e., EDZo), the contour covered a wide area around the tunnel, implying that the cracks had penetrated deep into the rock. Subsequently, the area covered by probability contours in cases 2 and 5 (i.e., EDZi) was narrower, while the cases 3 and 6 corresponded to a much limited area around the outer boundary of the tunnel (i.e., HDZ). The HDZ hosted relatively dense yet shallow cracks formed due to severe disturbance of the medium around the excavation region. 4. Discussion 4.1. Reliability sensitivity analysis An important point to address in this study is the effects of σmax and CI on the resulting failure probability. To deal with this issue, a parametric analysis was performed to see how the probability varied with the mean of each variable. For this purpose, the proposed algorithm was repeated for case 1 and the established reliability problem was analyzed by changing the mean of one variable while keeping the other variable constant. As a result, different exceedance probability curves were obtained for different mean values of the considered variable. The results for σmax and CI are shown in Figures 9a and 9b, respectively. 11
4
3
3
3
2
2
2
1
1
1
60
0
0
0
50
-1
-1
-1
-2
-2
-2
-3
-3
-3
-4 -4
-4 -4 -3 -2 -1 0 1 2 3 4 Normalized distance from tunnel center, x (mm)
-4 -4 -3 -2 -1 0 1 2 3 4 Normalized distance from tunnel center, x (mm)
90
80 70
40
Probability (%)
Normalized distance from tunnel center, y
4
30 20
-3 -2 -1 0 1 2 3 Normalized distance from tunnel center, x
4
(a) case 1
(b) case 2
10
(c) case 3 100
4
4
4
3
3
3
2
2
2
1
1
1
60
0
0
0
50
-1
-1
-1
-2
-2
-2
-3
-3
-3
-4 -4
-4 -4 -3 -2 -1 0 1 2 3 4 Normalized distance from tunnel center, x (mm)
-4 -4 -3 -2 -1 0 1 2 3 4 Normalized distance from tunnel center, x (mm)
90 80
70
40
Probability (%)
Normalized distance from tunnel center, y
100
4
30
20
-3 -2 -1 0 1 2 3 Normalized distance from tunnel center, x
(d) case 4
4
(e) case 5
10
(f) case 6
Figure 8: Contour plots of exceedance probability around tunnel excavation for cases 1 to 6.
As can be seen, the probability values increased with σmax , rising the exceedance probability curve to higher levels. On the contrary, the probability dropped sharply with increasing the value of CI, declining the exceedance probability curve. To further examine this issue and determine the extent to which the changes in each parameter affect the crack depth, the difference between the obtained probability curves was calculated. For instance, the difference between exceedance probability (∆P ) corresponding to different maximum tangential stresses was calculated as follows: ∆P = P (σmax ) − P (σmax = 30),
(8)
where P (σmax = 30) is the exceedance probability corresponding to σmax = 30 and P (σmax ) refers to the exceedance probability corresponding to σmax = 50, 70, 90, 110, 130, or 150. The results for σmax and CI are presented in Figures 9c and 9d, respectively. In the following, the position of maximum effect for each variable was estimated by finding the maximum point of each graph and fitting a curve. As can be seen, the σmax affected the cracks with normalized depths ranging between 1.3 and 1.7, while the maximum effect of CI was broader and covered the normalized crack depths ranging between 1.2 and 2.1.
12
100
100
Probability of exceedance (%)
80 70
60
Probability of exceedance (%)
P(σmax=150 MPa) P(σmax=130 MPa) P(σmax=110 MPa) P(σmax=90 MPa) P(σmax=70 MPa) P(σmax=50 MPa) P(σmax=30 MPa)
90
50 40
30 20
80
70 60 50
40 30
20 10
10 0
P(CI=70 MPa) P(CI=60 MPa) P(CI=50 MPa) P(CI=40 MPa) P(CI=30 MPa) P(CI=20 MPa) P(CI=10 MPa)
90
0.5
1
1.5 2 2.5 3 3.5 Normalized damage depth, rd/R
4
0
4.5
0
100 P(σmax=150)-P(σmax=30) P(σmax=130)-P(σmax=30) P(σmax=110)-P(σmax=30) P(σmax=90)-P(σmax=30) P(σmax=70)-P(σmax=30) P(σmax=50)-P(σmax=30) Location of maximum points
90
80 70 60
50 40
30 20 10
0
0.5
1
1.5 2 2.5 3 3.5 Normalized damage depth, rd/R
4
2 3 4 Normalized damage depth, rd/R
5
6
(b) Probability exceedance curve for different values of CI
4.5
(c) Difference between exceedance probability curves corresponding to different values of σmax
Difference between exceedance probabilities (%)
Difference between exceedance probabilities (%)
(a) Probability exceedance curve for different values of σmax
1
100 P(CI=10)-P(CI=70) P(CI=20)-P(CI=70) P(CI=30)-P(CI=70) P(CI=40)-P(CI=70) P(CI=50)-P(CI=70) P(CI=60)-P(CI=70) Location of maximum points
90 80
70 60
50 40 30
20 10 0
0
1
2 3 4 Normalized damage depth, rd/R
5
6
(d) Difference between exceedance probability curves corresponding to different values of CI
Figure 9: Effects of σmax and CI on exceedance probability for case 1.
4.2. Exceedance probability curve for other types of rocks So far, the exceedance probability has been calculated and discussed for the normalized damage depth in limestone. In this section, other types of rocks are examined using the model developed by Perras and Diederichs [11]. Based on the samples reported in this study, the probabilistic characteristics were calculated for two new rock types, as outlined in Table 4. The coefficients B and D (Eq. 1) for granite and mudstone are given in Table 5. Following the substitution of the coefficients B and D into the LSF function (Eq. 3), a new reliability problem was developed. To solve this problem, a new set of random values was generated for σmax and CI according to the rock types, and the Monte Carlo algorithm was utilized to calculate the exceedance probability curves. The results are shown in Figure 10. As can be inferred from Figure 10, the obtained probability for mudstone was larger than that for granite. In fact, the cracks in the mudstone had a greater chance to infiltrate through the rock, while the granite exhibited more resistance against the crack growth. This point was visualized by plotting probability contours around the tunnel (Figure 11). As 13
Table 4: The random variables assumed in the analysis.
Rock type Granite Mudstone
Variable Mean σmax (MPa) 220 CI (MPa) 125 σmax (MPa) 35 CI (MPa) 19
Standard deviation 30 16 11 6
Table 5: Values of B and D in different cases.
Case 7 8 9 10 11 12 13 14 15 16 17 18
Rock type
Granite
Mudstone
KHh 1.5 1.5 1.5 2 2 2 1.5 1.5 1.5 2 2 2
Zone B EDZo 0.62 EDZi 0.41 HDZ 0.09 EDZo 0.58 EDZi 0.36 HDZ 0.11 EDZo 0.71 EDZi 0.49 HDZ 0.2 EDZo 0.66 EDZi 0.39 HDZ 0.15
D 0.58 0.53 0.62 0.65 0.62 0.85 0.59 0.55 0.52 0.59 0.59 0.68
can be seen, in cases 7, 10, 13, and 16 (i.e., EDZo), wide areas were covered by probability contours. In comparison, the areas were drastically limited in cases 8, 11, 14, and 17 (i.e., EDZi) as well as cases 9, 12, 15, and 18 (i.e., HDZ). Additionally, the probability contours in cases 13 to 18 (i.e., mudstone) covered relatively wider areas than cases 7 to 12 (i.e., granite) and showed a slower decreasing trend as one moved away from the outer wall of the tunnel. This observation can be explained by the soft structure of mudstone, as compared to the granite, which allowed deep cracks to grow into the rock. The influence of in-plane stress ratio on the failure probability was investigated based on the resulting exceedance probability curves. For this purpose, the difference in the exceedance probability between different damage zones with different KHh values was calculated as follows: ∆P = P (KHh = 1.5) − P (KHh = 2).
(9)
This gave different graphs for different rock types (Figure 12). In Figure 12, the black, red, and blue curves correspond to the EDZo, EDZi, and HDZ, respectively. The first point to note from Figure 12 is that the effect of KHh was limited in granite, while it was more pronounced in the limestone and maximal in the mudstone. This was, of course, predictable as the hard structure of granite prevents the KHh to highly affect the crack propagation process. The limestone and mudstone were found to be more sensitive to 14
100
Exceedance probability curve for case 7 Exceedance probability curve for case 8 Exceedance probability curve for case 9
90 80
Probability of exceedance (%)
Probability of exceedance (%)
100
70 60 50 40 30 20 10 1
1.1
1.2
100 Probability of exceedance (%)
1.3 1.4 1.5 1.6 1.7 1.8 Normalized damage depth, rd/R (a) Cases 7 to 9 ( K Hh = 1.5)
1.9
70 60 50 40 30 20
0
2
80
1
1.1
1.2
100
Exceedance probability curve for case 13 Exceedance probability curve for case 14 Exceedance probability curve for case 15
90
70
60 50
40 30
20 10
0
80
10
Probability of exceedance (%)
0
Exceedance probability curve for case 10 Exceedance probability curve for case 11 Exceedance probability curve for case 12
90
1.3 1.4 1.5 1.6 1.7 Normalized damage depth, rd/R (b) Cases 10 to 12 (K Hh = 2)
1.8
1.9
2
Exceedance probability curve for case 16 Exceedance probability curve for case 17 Exceedance probability curve for case 18
90
80 70
60 50
40 30
20 10
1
1.5
2 2.5 Normalized damage depth, rd/R (c) Cases 13 to 15 (K Hh = 1.5)
3
3.5
01
1.5
2 2.5 Normalized damage depth, rd/R (d) Cases 16 to 18 (K Hh = 2)
3
3.5
Figure 10: Exceedance probability curves for (a) cases 7 to 9, (b) cases 10 to 12, (c) cases 13 to 15, and (d) cases 16 to 18.
KHh because of their softer structures. As another important point, the maximum impact of KHh on the EDZo was found to occur at large normalized crack depths approximately ranging between 1.5 and 2, while the maximum impacts of the KHh in the EDZi and HDZ occurred at the normalized crack depths ranging within 1.3-1.5 and 1-1.3, respectively. This finding seems reasonable since the EDZo inherently covers the deep cracks rooted in the rock, while the EDZo and HDZ tend to cover shallow yet dense crack systems. The probability differences in Figure 12 also show that the impact of KHh on the HDZ is much greater than those on the EDZi and then the EDZo. This result is understandable because the earliest impact of a change in the in-plane stress is imposed on the neighboring medium of the excavation hole, followed by a decrease in the impact with distance from the excavation wall. However, despite the general sense that an increase in the KHh can reduce the failure probability, the case observed in granite shows that the probability of failure can even increase with increasing the KHh , making the difference presented in Figure 12b negative. This issue might be related to the coefficients reported in Perras and Diederichs model [11]. There should be a particular factor affecting the regression models of B and D so that the Eq. 1 reports higher damage when KHh increases. Therefore, when using the 15
4
3
3
3
2
2
2
1
1
1
60
0
0
0
50
-1
-1
-1
-2
-2
-2
-3
-3
-3
-4 -4
-4 -4 -3 -2 -1 0 1 2 3 4 Normalized distance from tunnel center, x (mm)
-4 -4 -3 -2 -1 0 1 2 3 4 Normalized distance from tunnel center, x (mm)
80 70
20
-3 -2 -1 0 1 2 3 Normalized distance from tunnel center, x
4
(b) Case 8
10
(c) Case 9 100
4
4
4
3
3
3
2
2
2
1
1
1
60
0
0
0
50
-1
-1
-1
-2
-2
-2
-3
-3
-3
-4 -4
-4 -4 -3 -2 -1 0 1 2 3 4 Normalized distance from tunnel center, x (mm)
-4 -4 -3 -2 -1 0 1 2 3 4 Normalized distance from tunnel center, x (mm)
90 80 70
40
Probability (%)
Normalized distance from tunnel center, y
40 30
30 20
-3 -2 -1 0 1 2 3 Normalized distance from tunnel center, x
4
(d) Case 10
(e) Case 11
10
(f) Case 12 100
4
4
4
3
3
3
2
2
2
1
1
1
60
0
0
0
50
-1
-1
-1
-2
-2
-2
-3
-3
-3
-4 -4
-4 -4 -3 -2 -1 0 1 2 3 4 Normalized distance from tunnel center, x (mm)
-4 -4 -3 -2 -1 0 1 2 3 4 Normalized distance from tunnel center, x (mm)
90 80 70
40
Probability (%)
Normalized distance from tunnel center, y
90
Probability (%)
Normalized distance from tunnel center, y
4
(a) Case 7
30 20
-3 -2 -1 0 1 2 3 Normalized distance from tunnel center, x
4
(g) Case 13
(h) Case 14
10
(i) case 15 100
4
4
4
3
3
3
2
2
2
1
1
1
60
0
0
0
50
-1
-1
-1
-2
-2
-2
-3
-3
-3
-4 -4
-4 -4 -3 -2 -1 0 1 2 3 4 Normalized distance from tunnel center, x (mm)
-4 -4 -3 -2 -1 0 1 2 3 4 Normalized distance from tunnel center, x (mm)
90
80 70
40 30 20
-3 -2 -1 0 1 2 3 Normalized distance from tunnel center, x
(j) Case 16
4
(k) Case 17
10
(l) Case 18
Figure 11: Contour plots of exceedance probability around the excavation for cases 7 to 18.
16
Probability (%)
Normalized distance from tunnel center, y
100
4
Difference between exceedance probabilities (%)
Difference between exceedance probabilities (%)
60 Difference between cases 1 and 4 Difference between cases 2 and 5 Difference between cases 3 and 6
50 40
30 20 10 0 -10
0
0.5
1
1.5 2 2.5 3 Normalized damage depth, rd/R
3.5
4
25 Difference between cases 7 and 10 Difference between cases 8 and 11 Difference between cases 9 and 12
20
15 10
5 0 -5
-10 -15
-20 -25
0
0.5
1 1.5 Normalized damage depth, rd/R
2
(b) Granite Difference between exceedance probabilities (%)
(a) Limestone 25 Difference between cases 13 and 16 Difference between cases 14 and 17 Difference between cases 15 and 18
20
15 10 5
0 -5
0
0.5
1
1.5 2 2.5 3 3.5 Normalized damage depth, rd/R
4
4.5
5
(c) Mudstone
Figure 12: Effects of in-plane stress ratio in (a) limestone, (b) granite, and (c) mudstone.
Perras and Diederichs model [11] in granite for estimating HDZ zone, further caution should be considered, with the coefficients B and D to be calibrated by field tests at the specific site under investigation wherever possible. 5. Conclusion In this study, excavation-induced failure probability was investigated in three different types of rock including granite, mudstone, and limestone. For this purpose, a probabilistic model was established based on the model proposed by Perras and Diederichs. Then, the Monte Carlo sampling method was adopted to approximate the exceedance probability for different induced crack depths, with the results presented in the form of exceedance probability curves. Effects of the maximum tangential stress and the crack initiation threshold on failure probability were also examined using reliability sensitivity analysis. The main conclusions drawn from this study are as follows: As the depth of the crack increases, the probability falls sharply; as a result, the exceedance probability is less than 5% for the normalized crack depths of more than 2.28.
17
The effect of σmax is maximal at normalized crack depths ranging between 1.5 and 1.7, while that of CI is maximal when the normalized crack depth falls between 1.2 and 2.1. A change of in-plane stress ratio shows the most significant impacts on the EDZo when the normalized crack depth ranges between 1.5 and 2, while its maximum effects on the EDZi and HDZ is estimated to occur at the normalized crack depths ranging within 1.3-1.5 and 1-1.3, respectively. Increasing the in-plane stress ratio reduces the failure probability and causes the exceedance probability curve to locate at lower positions.
It is noteworthy that the analyses presented in this study are based on the assumption that the involved variables are independent. Estimating the failure probability around the excavation hole considering possible correlations between the involved variables is an interesting subject that will be investigated by the authors in their future research. References [1] L. Ding, L. Ma, H. Luo, M. Yu, X. Wu, Wavelet analysis for tunneling-induced ground settlement based on a stochastic model, Tunnelling and Underground Space Technology 26 (2011) 619–628. [2] Y. Fan, W. Lu, Y. Zhou, P. Yan, Z. Leng, M. Chen, Influence of tunneling methods on the strainburst characteristics during the excavation of deep rock masses, Engineering Geology 201 (2016) 85–95. [3] K. Zhang, J. L. Chavez Torres, Z. Zang, Numerical analysis of pipelines settlement induced by tunneling, Advances in Civil Engineering 2019 (2019) 1–10. [4] R. Read, 20 years of excavation response studies at aecl’s underground research laboratory, International Journal of Rock Mechanics and Mining Sciences 41 (2004) 1251 – 1275. Rock Mechanics Results from the Underground Research Laboratory, Canada. [5] M. Tao, Z. Hong, K. Peng, P. Sun, M. Cao, K. Du, Evaluation of excavation-damaged zone around underground tunnels by theoretical calculation and field test methods, Energies 12 (2019) 1–18. [6] L. X. Xie, Q. B. Zhang, J. C. Gu, W. B. Lu, S. Q. Yang, H. W. Jing, Z. L. Wang, Damage evolution mechanism in production blasting excavation under different stress fields, Simulation Modelling Practice and Theory 97 (2019) 101969. [7] S. Wang, X. Li, J. Yao, F. Gong, X. Li, K. Du, M. Tao, L. Huang, S. Du, Experimental investigation of rock breakage by a conical pick and its application to non-explosive mechanized mining in deep hard rock, International Journal of Rock Mechanics and Mining Sciences 122 (2019) 104063. 18
[8] H. Wang, Y. Jiang, S. Xue, B. Shen, C. Wang, J. Lv, T. Yang, Assessment of excavation damaged zone around roadways under dynamic pressure induced by an active mining process, International Journal of Rock Mechanics and Mining Sciences 77 (2015) 265– 277. [9] G. Walton, M. Lato, H. Ansch¨ utz, M. A. Perras, M. S. Diederichs, Non-invasive detection of fractures, fracture zones, and rock damage in a hard rock excavation — ¨ o Hard Rock Laboratory in Sweden, Engineering Geology 196 Experience from the Asp¨ (2015) 210–221. [10] T. Siren, P. Kantia, M. Rinne, Considerations and observations of stress-induced and construction-induced excavation damage zone in crystalline rock, International Journal of Rock Mechanics and Mining Sciences 73 (2015) 165–174. [11] M. A. Perras, M. S. Diederichs, Predicting excavation damage zone depths in brittle rocks, Journal of Rock Mechanics and Geotechnical Engineering 8 (2016) 60–74. [12] M. A. Perras, M. S. Diederichs, A review of excavation damage zones in sedimentary rocks with emphasis on numerical modelling for edz definition, The Canadian Geotechnical Society’s Annual Conference - Geo2010, Alberta, Canada. [13] E. Ghazvinian, Fracture initiation and propagation in low porosity crystalline rocks: implications for excavation damage zone (EDZ) mechanics, Phd thesis, Queen’s University, 2015. [14] S. Kwon, C. S. Lee, S. J. Cho, S. W. Jeon, W. J. Cho, An investigation of the excavation damaged zone at the KAERI underground research tunnel, Tunnelling and Underground Space Technology 24 (2009) 1–13. [15] W.-B. Lu, Y.-G. Hu, J.-H. Yang, M. Chen, P. Yan, Spatial distribution of excavation induced damage zone of high rock slope, International Journal of Rock Mechanics and Mining Sciences 64 (2013) 181–191. [16] W.-J. Cho, J.-S. Kim, C. Lee, H.-J. Choi, Gas permeability in the excavation damaged zone at KURT, Engineering Geology 164 (2013) 222–229. [17] H. Q. Yang, Y. Y. Zeng, Y. F. Lan, X. P. Zhou, Analysis of the excavation damaged zone around a tunnel accounting for geostress and unloading, International Journal of Rock Mechanics and Mining Sciences 69 (2014) 59–66. [18] W. C. Zhu, J. Wei, J. Zhao, L. L. Niu, 2D numerical simulation on excavation damaged zone induced by dynamic stress redistribution, Tunnelling and Underground Space Technology 43 (2014) 315–326. [19] Z. Wu, Y. Jiang, Q. Liu, H. Ma, Investigation of the excavation damaged zone around deep TBM tunnel using a Voronoi-element based explicit numerical manifold method, International Journal of Rock Mechanics and Mining Sciences 112 (2018) 158–170. 19
[20] Y. Xue, F. Gao, X. Liu, X. Liang, Permeability and pressure distribution characteristics of the roadway surrounding rock in the damaged zone of an excavation, International Journal of Mining Science and Technology 27 (2017) 211–219. [21] A. A. Rodriguez, U. Kuhlmann, P. Marschall, 3D modelling of the Excavation Damaged Zone using a Marked Point Process technique, Geomechanics for Energy and the Environment 17 (2019) 29–46. [22] C. D. Martini, R. S. Read, J. B. Martino, Observations of brittle failure around a circular test tunnel, International Journal of Rock Mechanics and Mining Sciences 34 (1997) 1065–1073. [23] A. Fakhimi, F. Carvalho, T. Ishida, J. F. Labuz, Simulation of failure around a circular opening in rock, International Journal of Rock Mechanics and Mining Sciences 39 (2002) 507–515. [24] E. Hrubesova, L. Duris, Assessment of tunnel’s face support pressure, IOP Conference Series: Materials Science and Engineering 236 (2017) 1–7. [25] N. Dadashzadeh, M. S. Diederichs, Reliability of prediction for tunnel excavation damage zone depth in brittle rocks, in: ISRM International Symposium - 10th Asian Rock Mechanics Symposium, International Society for Rock Mechanics and Rock Engineering, Singapore, 2018. [26] B. Tang, H. Cheng, Y. Tang, Z. Yao, M. Li, T. Zheng, J. Lin, Excavation damaged zone depths prediction for TBM-excavated roadways in deep collieries, Environmental Earth Sciences 77 (2018) 1–14. [27] I. Vazaios, N. Vlachopoulos, M. Diederichs, Assessing fracturing mechanisms and evolution of excavation damaged zone of tunnels in interlocked rock masses at high stresses using a finite-discrete element approach, Journal of Rock Mechanics and Geotechnical Engineering 11 (2019) 701–722. [28] I. Vazaios, N. Vlachopoulos, M. S. Diederichs, Mechanical analysis and interpretation of excavation damage zone formation around deep tunnels within massive rock masses using hybrid finite–discrete element approach: case of Atomic Energy of Canada Limited (AECL) Underground Research Laboratory (URL) test tunnel, Canadian Geotechnical Journal 56 (2019) 35–59. [29] Q. Xie, K. Peng, Space-time distribution laws of tunnel excavation damaged zones (EDZs) in deep mines and EDZ prediction modeling by random forest regression, Advances in Civil Engineering 2019 (2019) 1–13. [30] J. C. Kok-Kwang Phoon, Risk and reliability in geotechnical engineering, CRC Press, Florida, 2014.
20
[31] G. Baecher, J. Christian, Reliability and statistics in geotechnical engineering, Wiley, 2003. [32] M. Shadab Far, Y. Wang, Approximation of the Monte Carlo sampling method for reliability analysis of structures, Mathematical Problems in Engineering 2016 (2016). [33] D. P. Kroese, T. Taimre, Z. I. Botev, Handbook of Monte Carlo methods, Wiley, 2011.
21
Highlights - The depth of excavation-induced damage zones is established as a reliability problem. - The exceedance probability curve is developed for three rock types, including limestone, granite, and mudstone. - Failure probability contours are presented around the excavation region. - The effects of maximum tangential stress and the crack initiation threshold are investigated.
Declaration of interests ☒ The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper. ☐The authors declare the following financial interests/personal relationships which may be considered as potential competing interests:
Declarations of interest: none