- Email: [email protected]
Developing an algorithm to estimate in situ stresses using a hybrid numerical method based on local stress measurement
International Journal of Rock Mechanics & Mining Sciences 55 (2012) 80–85
Contents lists available at SciVerse ScienceDirect
International Journal of Rock Mechanics & Mining Sciences journal homepage: www.elsevier.com/locate/ijrmms
Developing an algorithm to estimate in situ stresses using a hybrid numerical method based on local stress measurement Z. Khademian a, K. Shahriar a,n, M. Gharouni Nik b a b
Department of Mining and Metallurgy, Amirkabir University of Technology (Tehran Polytechnic), Hafez Ave., Tehran, Iran School of Railway Engineering, IUST University, Tehran, Iran
a r t i c l e i n f o
a b s t r a c t
Article history: Received 2 March 2011 Received in revised form 2 April 2012 Accepted 15 May 2012 Available online 20 July 2012
The main contribution of this paper is to develop an algorithm to estimate in situ stress using a hybrid numerical method based on data of local stress measurement. Suggested hybrid method, first, estimates field stresses using an adverse Boundary Element Method (BEM3D) and local stresses are then calculated by Finite Element (FEM3D) or Finite Difference Methods (FDM3D) under certain conditions. Also it is possible to reach arbitrary accuracy for stress estimation. To show the applicability of this algorithm, it is applied for the Seymare Dam (SD) site to determine the stress field. Based on the local stresses measured by Hydraulic Fracturing (HF), the stress state at SD is estimated by the proposed algorithm in which the far-field stress state has been calculated using a 3D boundary element code, after which the local in situ stresses are calculated using a 3D finite element code. Finally, the results are compared with those elicited from the HF method. The finite element results, considering 10% allowable error, roughly coincides with the measured results. & 2012 Elsevier Ltd. All rights reserved.
Keywords: In situ stress Back analysis Numerical methods Field measurement Simulation
1. Introduction Recent experience with underground research has highlighted the necessity of understanding the in situ stress state for designing and constructing underground and open excavations, especially with regard to excavation stability and hydraulic suitability. Such projects require stress evaluations for rock masses with widths of several kilometers and depths greater than 1 km [1]. Analytical solutions for stresses exist for a range of idealized twodimensional (2D) topographies, but accurate analytical solutions are not valuable for an elastic–plastic property with an irregular three-dimensional (3D) surface [2]. Such terrains may be analyzed with domain numerical methods, such as Finite Difference Methods (FDM) or Finite Element Methods (FEMs) [3]. A problematic aspect of both these methods, however, is that the model boundary conditions are applied to boundaries of arbitrarily defined geometry, and it is not clear what these should be in reality. Hence, the main question is how the boundary conditions of a numerical model of a rock mass should be determined. In contrast, Boundary Element Methods (BEMs) have the advantages that no artificial truncated boundaries are necessary and therefore, BEMs give us more accurate results on open regions or deep underground but they are not such applicable for cases with complex geologic structure [4].
n
Corresponding author. Tel.: þ98 216 4542972. E-mail address: [email protected] (K. Shahriar).
1365-1609/$ - see front matter & 2012 Elsevier Ltd. All rights reserved. http://dx.doi.org/10.1016/j.ijrmms.2012.05.019
In this paper, a numerical inverse method is presented to characterize the state of in situ stresses in a rock mass aiming at determining the boundary conditions considering tectonic stresses, topography, anisotropy and heterogeneity of rock mass. Unlike back analysis done by FEM or FDM in which strain components or displacement ones are hypothesized as unknown parameters on model boundaries [5], we carried out a back analysis using the BEM taking the far-field stress state unknown. The BEM back analysis results can then serve as the input boundary conditions for FEM or FDM analysis, two popular and powerful methods well suited for non-linear or large-strain problems. By determination of the model boundary conditions, we can easily study on a variety of key issues, such as excavationinduced stress or the stress perturbations caused by faults. Finally, the suggested method is applied to calculate the stress field of Seymare Dam district, one of the greatest dam projects in Ilam, south–east of Iran. The numerical model is 2.1 km 1.7 km in area. Note that in this site, stress measurements are obtained in fifteen different locations by Hydraulic Fracturing methods.
2. Overview of method 2.1. Basic relationship The in situ stress field can never be completely measured, and as a result, the methods developed for evaluating it include simplifying assumptions [6]. In practice, the most common
Z. Khademian et al. / International Journal of Rock Mechanics & Mining Sciences 55 (2012) 80–85
principal stresses, and szz symbolizes the vertical stress at depth h. Actual field measurements, however, show the horizontal stresses commonly do not follow Eq. (1), and are several times larger than the vertical stress in many places [1,8,9]. Furthermore, the horizontal stress state is not uniform, especially in a seismically active region [10]. The current stress state reflects a series of historic and geologic events, the stiffness of a rock mass as well as gravitational and tectonic stresses [11]. Also, rock mass heterogeneities and discontinuities lead to local stress variations [7]. The following equations provide a more general relationship with the constraint that one of the principal stress components is vertical and varies linearly with respect to depth h:
Start
Stage 1: Selecting set (N-T) of In situ stress measurement data
Stage 2: Modeling the district based on BEM3D
Stage 3: Induced Stresses Determination under different far field stress components
s0 ¼ S þ rgha þ 2
2 Re-Calculating a and b
2
No
Stage 8: Comparing the calculated local stress with those of set (T)
No
β ≤ 0.1 Yes End
Fig.1. Proposed algorithm for estimating in situ stresses.
supposition is to assume lateral confinement (i.e., there is no horizontal displacement anywhere due to gravitational loading). Based on the mentioned assumption, the following equation is often utilized to describe the initial state of stress in a uniform soil or rock mass beneath a horizontal free surface [7]:
s0zz ¼ rgh, s0xx ¼ s0yy ¼
n 1n
rgh
ð2aÞ
3
Syy
7 05
0
0
ð2bÞ
3
kxx
kxy
0
ð2cÞ
3
kyy
7 05
0
0
ð2dÞ
where s0 and S stand for the in situ measured and the tectonic stress states, respectively. Also, a is a 3-by-3 matrix of coefficients due to gravitational loading, and k is a 3-by-3 matrix of coefficients for a linear distribution of horizontal stresses [1,12]. The parameters a and b are constants determined by measured data using a simplex method that allow for a non-linear distribution of certain stresses [9]. Since there are many factors affecting the variation of the in situ stress, and owing to the fact that only four factors are here used in calculating the stress state, therefore, some other parameters, e.g. process of deposition and erosion, change in the earth’s rotation axis, tidal wave pressure, are certainly effective which are unknown. Now, two parameters a and b are considered as the representative of the unknown parameters and then obtained after several trial and errors such that the finally obtained values (for a and b) guarantee the accuracy of the proposed algorithm. In Eq. (2), {Sx,Sy,Sxy,Syx, ax,ay,az, kx,ky,kxy,kyx} are considered as coefficients of far-field stress components. Considering the symmetry of matrices S and k, there are the only nine different field stress coefficients. The right-hand side of Eq. (2) contains a constant term, a term that varies linearly with h, and a non-linear one with a hyperbolic coefficient of b/(hþa).
Yes
Second Provision
0
axx 0 0 ayy 0 7 5 azz 0
6 k ¼ 4 kyx 0
Stage 7: Comparing the calculated local stress with the measured ones from set (N-T)
γ ≤ 0.1
Sxy
b rghk hþa
a¼6 40 0
Stage6: Re-simulating the district by FEM or FDM and determining the local stress
First Provision
Sxx
6 S ¼ 4 Syx 0
Stage 4: Coefficient matrix establishment
Stage5: Boundary Conditions Calculation
81
ð1Þ
where r is the mass density of the materials, g is referred to as the gravitational acceleration, n stands for Poisson’s ratio, and h is the depth from surface [7]. The terms sxx and syy are two horizontal
2.2. Back analysis In the majority of forward analyses of stress using numerical methods, after modeling the district, boundary conditions are applied to the model and local stresses are then computed [5]; but, in this paper, it is shown how to calculate boundary conditions by in situ stresses measured by HF method in certain location using a back analysis. It is noted that the most important stage of this analysis is selection of the more compatible numerical method with the district structure.
82
Z. Khademian et al. / International Journal of Rock Mechanics & Mining Sciences 55 (2012) 80–85
3. Description of the proposed method
Table 2 Induced stress states under different far-field stress component.
In this section, the suggested algorithm is described step-bystep and its applicability is shown by highlighting the most important stages. Fig. 1 illustrates the flowchart of the method. 3.1. Stage 1 In this stage, a set (N–T) must be selected in which N ¼ {p1,p2,p3,...,pn}, TCN and {p1, p2,y,pn} are points whose stresses are measured by HF. In other words, T is the set of the points whose stresses will be compared with the corresponding stresses calculated in Stage 9. Therefore, N–T is referred to as the set of the points used as input data in the algorithm. Also, the numbers of the members of (N–T) set must satisfy the condition: nt Z 9
ð3Þ
where n–t is the number of the members of set (N–T). 3.2. Stage 2 In simulation, the first significant step is choosing the more consistent numerical method with the rock mass structure whereas being able to satisfy the simulation purpose [13,14]. Table 1 reports the basic properties of different numerical methods. According to this table, the more applicable method(s) are proposed based on certain properties. Note that the proposed methods can be applied using different valid software. The scope of this stage is modeling the certain district for calculating boundary conditions thereof. Therefore, according to Table 1, as the far-field stresses of the district are unknown, the most suitable method is BEM. Also, to apply the irregular topography of the district, the 3D modeling is preferred instead of 2D. In addition, to reduce the complexity of the calculations, it is better to impose the homogeneous BEM3D. Since heterogeneity extremely affects in situ stress values [15], it is tried to be considered in other stages. 3.3. Stage 3 In this stage, we try to apply hypothetical far-field stresses on the model. In other words, by exploiting numerical modeling in stage 2, we can compute the induced stresses at a certain point via applying each far-field stress component introduced in Eq. (2) to BEM3D model. First, each coefficient of far-field stress component is supposed to be equal to 1, such that each component acts on the model independently. The induced stresses state is consequently calculated at set (N–T) points based on BEM. Table 2 shows nine independently different field stress components for calculating the induced stresses. So, the stress tensor includes 3 stress components at each point. It is noted that sxi,Sx symbolizes the stress along the x-axis at point i under the unit tectonic stress component Sx and iA(N T).
Massive rock Rock mass containing joint set with high spacing Rock mass containing joint set with low spacing Jointed rock Complex geological structure District with unknown far field stresses
Applied far field stress component
Induced state of stresses at point i
Sx Sy Sxy
Sx ¼ 1 Sy ¼ 1 Sxy ¼ 1 (rgh)x (rgh)y (rgh)z
sxi,Sx sxi,Sy sxi,Sxy sxi,ax sxi,ay sxi,axy sxi,kx sxi,ky sxi,kxy
axrgh ayrgh azrgh b k rgh aþh x b k rgh aþh y b k rgh a þ h xy
ða þb h rghÞx ða þb h rghÞy ða þb h rghÞxy
syi,Sx syi,Sy syi,Sxy syi,ax syi,ay syi,axy syi,kx syi,ky syi,kxy
txyi,Sx txyi,Sy txyi,Sxy txyi,ax txyi,ay txyi,axy txyi,kx txyi,ky txyi,kxy
3.4. Stage 4 In Stage 4, the coefficient matrix must be established. To reach this purpose, the stress states at a certain point should be expressed as the superposition of each field stress component multiplied by the individual stress obtained from numerical modeling based on BEM. For example, the stress component sxi at a certain point i can be written as
sxi ¼ Sx sxi,Sx þSy sxi,Sy þ Sxy sxi,Sxy þ ax sxi, ax þ þ ay sxi, ay þ axy sxi, axy þkx sxi,kx þky sxi,ky þ kxy sxi,kxy
ð4Þ
Assuming the measured stress component at a point i as (sxi)m, we can substitute (sxi)m ¼ sxi in Eq. (2) to obtain: ðsxi Þm ¼ Sx sxi,Sx þ Sy sxi,Sy þ Sxy sxi,Sxy þ ax sxi, ax þ ay sxi, ay þ axy sxi, axy þkx sxi,kx þ ky sxi,ky þkxy sxi,kxy
ð5Þ
where iA(N T). To summarize: 0 B B B B B B B B B B B B B B B B B B B @
sm x1 sm y1 sm x1 ^ ^ ^ ^ ^ ^
sm xy1 0
1 C C C C C C C C C C C C C C C C C C C A
sx1,Sx
Bs B y1,Sx B B txy1,Sx B B ^ B B B ^ B ¼B B ^ B B ^ B B B ^ B B ^ @
1
sx1,Sy sx1,Sxy sx1,ax sx1,ay sx1,axy sx1,kx sx1,kx sx1=x sy1,Sy sy1,Sxy sy1,ax sy1,ay sy1,axy sy1,kx sy1,ky sy1,kxy C C C txy1,Sy txy1,Sxy txy1,ax txy1,ay txy1,axy txy1,kx txy1,ky txy1,kxy C C ^
^
^
^
^
^
^
^
^
^
^
^
^
^
^
^
^ ^
^ ^
^ ^
^ ^
^ ^
^ ^
^ ^
^ ^
^
^
^
^
^
^
^
^
^
^
^
^
^
^
^
^
txyn,Sx txyn,Sy txyn,Sxy txyn,ax txyn,ay txyn,axy txyn,kx txyn,ky txyn,kxy
Table 1 Guidance for selection of numerical methods.
Rock mass properties
Far-field stress component
C C C C C C C C C C C C C C A
0 Rock mass behavior Elastic
Plastic
BEM FEM, BEM DEM –
FDM, FEM FDM, FEM DEM FDM, FEM
FDM, FEM BEM
1 Sx B Sy C B C B C B Sxy C B C B Sx C B C B C B Sy C B C BS C B xy C B C B Sx C B C BS C @ y A Sxy
ð6Þ
Z. Khademian et al. / International Journal of Rock Mechanics & Mining Sciences 55 (2012) 80–85
3.5. Stage 5
To begin, we define
To solve the equations set, since the number of equations is usually more than that of the unknown coefficients {Sx,Sy,Sxy, ax,ay,axy, kx,ky,kxy}, so, a unique solution may not be obtained and a least-squares estimation is, therefore, applied here to find the best solution for set of equations.
g¼
3.6. Stage 6 The purpose of this stage is the determination of the local stress using boundary conditions calculated by the above-mentioned BEM3D back analysis. In previous stages, the effects of different parameters including topography, anisotropy and tectonic stresses are applied on in situ stress variations. Now, we try to exert the effects of heterogeneity of rock mass such as rock properties as well as major and minor geological structures on in situ stresses estimation. It is necessary to model the same area of BEM3D analysis using another numerical method that is able to simulate complex geological structure. According to Table 1, Finite Element (FEM3D) or Finite Difference Method (FDM3D) could be suitable for re-simulation under certain conditions.
3.7. Stage 7 To increase the workability of this algorithm in estimation of in situ stress, it is necessary to compare the measured data (set (N–T)) with those gained from Re-simulation stage (Stage 6).
P
2 2 P ðsci s2 mi Þ ,
i ¼ 1
i ¼ 1
sci
N ¼ p1 ,p2 ,p3 ,:::,pn ,
83
( i A NT g, T N
ð7Þ
where sci and smi are calculated and measured stresses in point i, respectively. According to Eq. (7), if g r0.1, the next stage will be investigated; otherwise, it is addressed to re-calculate a and b coefficients by Eq. (2). Note that parameters a and b indicate a non-linear distribution of certain stresses in terms of depth h and should be determined using a simplex method to minimize the sum of the differences between stress values, given by Eq. (2), and the measured stresses by HF. If the first provision is not satisfied after 10 times of trials and errors, the user is allowed to go to the next step aiming at expediting the calculations. 3.8. Stage 8 As described before, n t numbers of stress measured points are applied to start the algorithm. After determination of local stresses by FEM or FDM (i.e. passing the first provision and finding the final values for a and b), the remaining stress values (t points) should be compared with those obtained from the algorithm (Stage 6) in the same points: P 2 2 1 ðsci smi Þ b ¼ i ¼P , iAT ð8Þ 2 i¼1
sci
To be able to proceed along the algorithm, the second provision must be satisfied, too. If b r0.1, the process is terminated and we will have the values of in situ stresses in every arbitrary point; otherwise, according to the second provision, another set (N–T) of measured points must be chosen from set N to restart the algorithm. Similar to the stopping condition for the first loop, if the second provision is not satisfied after L times, the user should decrease the accuracy of the estimation (i.e. increase the values of g and b). Also, the increase of the number of in situ stress measurements by HF is another solution. The parameter L is given by L¼
n! ðntÞ!
ð9Þ
4. Algorithm applicability To show the applicability of this algorithm, it is applied for Seymare Dam (SD) to determine its stress field using the suggested method and based on the local stresses measured by HF.
Fig. 2. Seymare Dam district and its power tunnel.
Table 3 Input material properties for Hybrid method [16]. Description
Density (t/m3)
Poisson’s ratio
Young’s modulus (GPa)
Shear and normal stiffness (MPa/m) 60
Upper Assmary
1 2 3
2.72 2.78 2.54
0.31 0.32 0.31
1000 950 150
Middle Assmary
4 5
2.71 2.68
0.3 0.28
600 1200
Lower Assmary
6
2.69
0.29
1800
84
Z. Khademian et al. / International Journal of Rock Mechanics & Mining Sciences 55 (2012) 80–85
4.1. Introduction of Seymare district
4.3. Estimation of far-field stress state by BEM3D
Fig. 2 shows the district of SD and Power Tunnel. Local hills approach 200 m–300 m in height. In this district, there are three main geologic units named upper Asmary, middle Asmary and lower Asmary featuring a maximum thickness of 150 m as well as three main normal faults with a strike of N100E, a dip of 60 and an estimated throw of about 30 m. Rock mass geotectonic and joints set properties are also shown in Table 3 [16].
The BEM3D model simulating SD district by Examine3D software is shown in Fig. 3, which is 2.1 km in length and 1.7 km in width. Examine3D is a software based on three-dimensional Boundary Element method which is applied in order to simulate the area of SD site (Stage 2). Also, it is possible to program the BEM3D in Matlab software and analyze the model in Ansys instead of using Examine3D software. The grid points and zones’ numbers of the model are 155,900 and 77,760, respectively. According to Stage 3, the induced stresses at thirteen point (set (N–T)) are calculated and the coefficients matrix are then established (Stage 4). Next, a least-squares approach is utilized to control the distribution of the error and to find the best approximation solution for this over-determined set of equations. Table 5 presents the results calculated by BEM3D for far-field stress state (Stage 5), in which y direction shows the North.
4.2. Local stress measurements using Hydraulic Fracturing Rock stress measurements have been performed at fifteen points in Power Tunnel Trade using Hydraulic Fracturing technique (HF) in SD. Thus n ¼15 and according to Stage 1, t parameter is therefore 1rt r6. To start the algorithm, we consider t ¼2. Table 4 reports the in situ stresses values measured by HF at these points. Note that the tensile stress is taken negative while the compressive stress has a positive sign. Table 4 In situ measured stresses for the SD district. Measured location
sv (Mpa)
sh (Mpa)
sH (Mpa)
sxy (Mpa)
Azimuth sH degree
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15
2.17 2.2 2.26 2.1 2.28 2.31 2.298 2.32 2.343 2.365 2.387 2.409 2.432 2.454 2.46
2.2 2.23 2.2 2.35 2.48 2.29 2.28 3 3.54 2.31 2.4 2.24 2.45 2.57 2.68
3.9 3.87 4.04 3.86 4.11 4.26 4.17 4.22 4.01 4.31 4.56 4.41 4.46 4.51 4.52
4.475 4.485 4.496 4.498 4.515 4.524 4.533 4.22 4.42 4.48 4.6 4.3 4.7 4.6 4.8
43 52 51 45 56 52 44 56 50 43 55 46 55 57 44
4.4. Estimation of local stress state by FEM3D According to Table 1 and considering the existent discontinuities in SD district, the most compatible method for re-simulating, is FEM. In order to illustrate how the far-field stress coefficients can be adopted as input boundary conditions (Stage 6), the Plaxis3D software is applied to re-simulate the same area used in BEM3D analysis. As known, Plaxis3D is the software based on FEM3D to solve complex problems in rock mechanics, and in this simulation. The Plaxis3D model features 2.1 km long, 1.7 km wide and 900 m deep. The model is composed of 280,160 zones, 301,700 grid points and fifteen interfaces comprising seven faults, five layer and three joint sets as well. 4.5. Investigation of the first provision As described before, thirteen stress measurement points are applied to start the algorithm. In this stage, the values of a and b are considered equal to 904.6 m and 812.1 m, respectively, as initial values [1]. The first try results in g Z0.1; therefore, a and b must been modified using a simplex method. The new values of a and b, 890.34 and 786.45, cause g to be decreased to 0.08 and the first provision is therefore satisfied (Stage 7). 4.6. Investigation of the second provision After satisfying first provision, the remaining stress values (t ¼2) measured by HF should be compared with the estimated data in Stage 6. Following Stage 8 and using Eq. (8), b is obtained more than 0.1; thus, other set of measurement points must be selected. In case of SD district this loop repeats twice and finally by choosing t ¼3, b is reduced to 0.098. Table 6 shows the
Table 5 Determination of boundary condition by far field stresses.
Fig. 3. Discretization of the boundaries into triangular leaf elements.
Far field stress component
Coefficient of far field stress component
Calculated far field stress component
Sx Sy Sxy
Sx ¼1.68 Sy ¼3.82 Sxy ¼ 0.695 ax ¼0.778 ay ¼ 0.674 az ¼0.967 kx ¼ 0.667
1.68 3.82 0.695 0.778rgh 0.674rgh 0.967rgh 0:667 b þa h rgh
ky ¼ 0.536
0:536 b þa h rgh
kxy ¼ 0.236
0:236 b þa h rgh
axrgh ayrgh axyrgh kx b þa h rgh ky b þa h rgh kxy b þa h rgh
Z. Khademian et al. / International Journal of Rock Mechanics & Mining Sciences 55 (2012) 80–85
Table 6 Comparison of measured in situ stresses (M) with those calculated (Cal) using the proposed algorithm in SD district (MPa). Location
sh Cal
1 2 2 1.8 3 2.4 4 1.99 5 2.22 6 2.369 7 2.3 8 2.8 9 3.401 10 2.4 11 2.34 12 2 13 2.1 14 2.3 15 2.61 g ¼ 0.08, b ¼ 0.098
sH
sxy
M
Cal
M
Cal
M
2.2 2.235 2.2 2.35 2.48 2.295 2.277 3 3.54 2.318 2.4 2.244 2.458 2.571 2.68
3.6 4 4 3.7 4.039 4.1 3.7 4 3.8 4.1 4.3 4.6 4.5 4.5 4.4
3.9 3.87 4.04 3.863 4.108 4.26 4.174 4.222 4.01 4.318 4.567 4.415 4.463 4.511 4.521
4.117 4.376 4.66 4.201 4.68 4.74 4.6 4.44 4.3 4.35 4.7 4.5 4.6 4.357 4.67
4.475 4.485 4.496 4.498 4.515 4.524 4.533 4.22 4.42 4.48 4.6 4.3 4.7 4.6 4.8
comparison of measured in situ stresses by HF with the finally calculated values by the proposed algorithm in SD. This comparison shows how acceptable the results of the proposed algorithm are in the sense of calculating the in situ stresses.
5. Conclusions A new algorithm for determination of the far-field and local stresses based on stress measurements just in a few regions is developed in this paper. The algorithm considers the tectonic stresses, topography, anisotropy and heterogeneity of the rock mass as effective parameters. The measured in situ stress is decomposed into tectonic and gravitational components, and nine independent far-field stress coefficients are achieved by back analysis using BEM3D simulation and the local stresses are then obtained using FEM or FDM based on the back analysis results. The proposed hybrid method is, also, applied on SD site
85
and the results show how the calculated in situ stresses coincide well with the measured ones. Also, it is shown if the number of stress measurement data is adequate, the stress state will be calculated with the accuracy of 10% in an unknown location. Moreover, it is possible to achieve an estimated in situ stress with desired error by changing the provisions in the algorithm.
References [1] Li G, Mizuta Y, Ishida T, Li H, Nakama S, Sato T. Stress field determination from local stress measurements by numerical modeling. Int J Rock Mech Min Sci 2009;30:111–23. [2] Savage WZ, Swolfs HS, Powers PS. Gravitational stresses in long symmetric ridges and valleys. Int J Rock Mech Min Sci 1985;22:291–302. [3] Hakala M, Kuula H, Hudson JA. Estimating the transversely isotropic elastic intact rock properties for in situ stress measurement data reduction: A case study of the Olkiluoto mica gneiss, Finland. Int J Rock Mech Min Sci 2007;44:14–46. [4] Becker AA. The boundary element method in engineering. London: McGrawHill; 1992. [5] Aklin JE. Finite element analysis with error estimators. Oxford: Elsevier; 2005. [6] Peng S, Zhang J. Engineering geology for underground rocks. New York: Almas Schimmel; 2007. [7] Amadei B, Stephenson O. Rock mechanics for underground mining. London: Chapman & Hall; 1993. [8] Jaeger JC, Cook NGW, Zimmerman RW. Fundamentals of rock mechanics. 4th ed. Oxford: Wiley-Blackwell; 2007. [9] Brown ET, Hoek E. Trends in relationships between measured in situ stresses and depth. Int J Rock Mech Min Sci 1987;15:211–5. [10] Kanagawa T, Hibino S, Ishida T, Hayashi M, Kitahara Y. In situ stress measurements in the Japanese islands-over coring results from a multielement Gauge used at 23 sites. Int J Rock Mech Min Sci 1986;23:29–39. [11] Brady B, Brown ET. Rock mechanics for underground mining. London: Chapman & Hall; 1993. [12] Cornet FH, Valette B. In situ stress determination from Hydraulic Injection test data. J Geophys Res 1984;89:527–37. [13] Craig D, Richard A, Regueiro I, Borja Ronaldo I. Modelling brittle fracture, slip weakening and variable friction in geomaterials with an embedded strong discontinuity Finite Element. US Department of Energy: Sandia National Laboratories; 2006. [14] Muir D. Geotechnical modeling. London: Chapman & Hall; 2004. [15] Amadei B, Savage WZ. Effect of joints on rock mass strength and deformability. In: Hudson JA, editor. Comprehensive Rock Engineering. Oxford: Pergamon Press; 1993. [16] Iran Water & Power Resources Development Co. The in situ stress measurement result in Seymare Dam site by Hydraulic Fracturing. 3rd Report. Iran; 2010.
Contents lists available at SciVerse ScienceDirect
International Journal of Rock Mechanics & Mining Sciences journal homepage: www.elsevier.com/locate/ijrmms
Developing an algorithm to estimate in situ stresses using a hybrid numerical method based on local stress measurement Z. Khademian a, K. Shahriar a,n, M. Gharouni Nik b a b
Department of Mining and Metallurgy, Amirkabir University of Technology (Tehran Polytechnic), Hafez Ave., Tehran, Iran School of Railway Engineering, IUST University, Tehran, Iran
a r t i c l e i n f o
a b s t r a c t
Article history: Received 2 March 2011 Received in revised form 2 April 2012 Accepted 15 May 2012 Available online 20 July 2012
The main contribution of this paper is to develop an algorithm to estimate in situ stress using a hybrid numerical method based on data of local stress measurement. Suggested hybrid method, first, estimates field stresses using an adverse Boundary Element Method (BEM3D) and local stresses are then calculated by Finite Element (FEM3D) or Finite Difference Methods (FDM3D) under certain conditions. Also it is possible to reach arbitrary accuracy for stress estimation. To show the applicability of this algorithm, it is applied for the Seymare Dam (SD) site to determine the stress field. Based on the local stresses measured by Hydraulic Fracturing (HF), the stress state at SD is estimated by the proposed algorithm in which the far-field stress state has been calculated using a 3D boundary element code, after which the local in situ stresses are calculated using a 3D finite element code. Finally, the results are compared with those elicited from the HF method. The finite element results, considering 10% allowable error, roughly coincides with the measured results. & 2012 Elsevier Ltd. All rights reserved.
Keywords: In situ stress Back analysis Numerical methods Field measurement Simulation
1. Introduction Recent experience with underground research has highlighted the necessity of understanding the in situ stress state for designing and constructing underground and open excavations, especially with regard to excavation stability and hydraulic suitability. Such projects require stress evaluations for rock masses with widths of several kilometers and depths greater than 1 km [1]. Analytical solutions for stresses exist for a range of idealized twodimensional (2D) topographies, but accurate analytical solutions are not valuable for an elastic–plastic property with an irregular three-dimensional (3D) surface [2]. Such terrains may be analyzed with domain numerical methods, such as Finite Difference Methods (FDM) or Finite Element Methods (FEMs) [3]. A problematic aspect of both these methods, however, is that the model boundary conditions are applied to boundaries of arbitrarily defined geometry, and it is not clear what these should be in reality. Hence, the main question is how the boundary conditions of a numerical model of a rock mass should be determined. In contrast, Boundary Element Methods (BEMs) have the advantages that no artificial truncated boundaries are necessary and therefore, BEMs give us more accurate results on open regions or deep underground but they are not such applicable for cases with complex geologic structure [4].
n
Corresponding author. Tel.: þ98 216 4542972. E-mail address: [email protected] (K. Shahriar).
1365-1609/$ - see front matter & 2012 Elsevier Ltd. All rights reserved. http://dx.doi.org/10.1016/j.ijrmms.2012.05.019
In this paper, a numerical inverse method is presented to characterize the state of in situ stresses in a rock mass aiming at determining the boundary conditions considering tectonic stresses, topography, anisotropy and heterogeneity of rock mass. Unlike back analysis done by FEM or FDM in which strain components or displacement ones are hypothesized as unknown parameters on model boundaries [5], we carried out a back analysis using the BEM taking the far-field stress state unknown. The BEM back analysis results can then serve as the input boundary conditions for FEM or FDM analysis, two popular and powerful methods well suited for non-linear or large-strain problems. By determination of the model boundary conditions, we can easily study on a variety of key issues, such as excavationinduced stress or the stress perturbations caused by faults. Finally, the suggested method is applied to calculate the stress field of Seymare Dam district, one of the greatest dam projects in Ilam, south–east of Iran. The numerical model is 2.1 km 1.7 km in area. Note that in this site, stress measurements are obtained in fifteen different locations by Hydraulic Fracturing methods.
2. Overview of method 2.1. Basic relationship The in situ stress field can never be completely measured, and as a result, the methods developed for evaluating it include simplifying assumptions [6]. In practice, the most common
Z. Khademian et al. / International Journal of Rock Mechanics & Mining Sciences 55 (2012) 80–85
principal stresses, and szz symbolizes the vertical stress at depth h. Actual field measurements, however, show the horizontal stresses commonly do not follow Eq. (1), and are several times larger than the vertical stress in many places [1,8,9]. Furthermore, the horizontal stress state is not uniform, especially in a seismically active region [10]. The current stress state reflects a series of historic and geologic events, the stiffness of a rock mass as well as gravitational and tectonic stresses [11]. Also, rock mass heterogeneities and discontinuities lead to local stress variations [7]. The following equations provide a more general relationship with the constraint that one of the principal stress components is vertical and varies linearly with respect to depth h:
Start
Stage 1: Selecting set (N-T) of In situ stress measurement data
Stage 2: Modeling the district based on BEM3D
Stage 3: Induced Stresses Determination under different far field stress components
s0 ¼ S þ rgha þ 2
2 Re-Calculating a and b
2
No
Stage 8: Comparing the calculated local stress with those of set (T)
No
β ≤ 0.1 Yes End
Fig.1. Proposed algorithm for estimating in situ stresses.
supposition is to assume lateral confinement (i.e., there is no horizontal displacement anywhere due to gravitational loading). Based on the mentioned assumption, the following equation is often utilized to describe the initial state of stress in a uniform soil or rock mass beneath a horizontal free surface [7]:
s0zz ¼ rgh, s0xx ¼ s0yy ¼
n 1n
rgh
ð2aÞ
3
Syy
7 05
0
0
ð2bÞ
3
kxx
kxy
0
ð2cÞ
3
kyy
7 05
0
0
ð2dÞ
where s0 and S stand for the in situ measured and the tectonic stress states, respectively. Also, a is a 3-by-3 matrix of coefficients due to gravitational loading, and k is a 3-by-3 matrix of coefficients for a linear distribution of horizontal stresses [1,12]. The parameters a and b are constants determined by measured data using a simplex method that allow for a non-linear distribution of certain stresses [9]. Since there are many factors affecting the variation of the in situ stress, and owing to the fact that only four factors are here used in calculating the stress state, therefore, some other parameters, e.g. process of deposition and erosion, change in the earth’s rotation axis, tidal wave pressure, are certainly effective which are unknown. Now, two parameters a and b are considered as the representative of the unknown parameters and then obtained after several trial and errors such that the finally obtained values (for a and b) guarantee the accuracy of the proposed algorithm. In Eq. (2), {Sx,Sy,Sxy,Syx, ax,ay,az, kx,ky,kxy,kyx} are considered as coefficients of far-field stress components. Considering the symmetry of matrices S and k, there are the only nine different field stress coefficients. The right-hand side of Eq. (2) contains a constant term, a term that varies linearly with h, and a non-linear one with a hyperbolic coefficient of b/(hþa).
Yes
Second Provision
0
axx 0 0 ayy 0 7 5 azz 0
6 k ¼ 4 kyx 0
Stage 7: Comparing the calculated local stress with the measured ones from set (N-T)
γ ≤ 0.1
Sxy
b rghk hþa
a¼6 40 0
Stage6: Re-simulating the district by FEM or FDM and determining the local stress
First Provision
Sxx
6 S ¼ 4 Syx 0
Stage 4: Coefficient matrix establishment
Stage5: Boundary Conditions Calculation
81
ð1Þ
where r is the mass density of the materials, g is referred to as the gravitational acceleration, n stands for Poisson’s ratio, and h is the depth from surface [7]. The terms sxx and syy are two horizontal
2.2. Back analysis In the majority of forward analyses of stress using numerical methods, after modeling the district, boundary conditions are applied to the model and local stresses are then computed [5]; but, in this paper, it is shown how to calculate boundary conditions by in situ stresses measured by HF method in certain location using a back analysis. It is noted that the most important stage of this analysis is selection of the more compatible numerical method with the district structure.
82
Z. Khademian et al. / International Journal of Rock Mechanics & Mining Sciences 55 (2012) 80–85
3. Description of the proposed method
Table 2 Induced stress states under different far-field stress component.
In this section, the suggested algorithm is described step-bystep and its applicability is shown by highlighting the most important stages. Fig. 1 illustrates the flowchart of the method. 3.1. Stage 1 In this stage, a set (N–T) must be selected in which N ¼ {p1,p2,p3,...,pn}, TCN and {p1, p2,y,pn} are points whose stresses are measured by HF. In other words, T is the set of the points whose stresses will be compared with the corresponding stresses calculated in Stage 9. Therefore, N–T is referred to as the set of the points used as input data in the algorithm. Also, the numbers of the members of (N–T) set must satisfy the condition: nt Z 9
ð3Þ
where n–t is the number of the members of set (N–T). 3.2. Stage 2 In simulation, the first significant step is choosing the more consistent numerical method with the rock mass structure whereas being able to satisfy the simulation purpose [13,14]. Table 1 reports the basic properties of different numerical methods. According to this table, the more applicable method(s) are proposed based on certain properties. Note that the proposed methods can be applied using different valid software. The scope of this stage is modeling the certain district for calculating boundary conditions thereof. Therefore, according to Table 1, as the far-field stresses of the district are unknown, the most suitable method is BEM. Also, to apply the irregular topography of the district, the 3D modeling is preferred instead of 2D. In addition, to reduce the complexity of the calculations, it is better to impose the homogeneous BEM3D. Since heterogeneity extremely affects in situ stress values [15], it is tried to be considered in other stages. 3.3. Stage 3 In this stage, we try to apply hypothetical far-field stresses on the model. In other words, by exploiting numerical modeling in stage 2, we can compute the induced stresses at a certain point via applying each far-field stress component introduced in Eq. (2) to BEM3D model. First, each coefficient of far-field stress component is supposed to be equal to 1, such that each component acts on the model independently. The induced stresses state is consequently calculated at set (N–T) points based on BEM. Table 2 shows nine independently different field stress components for calculating the induced stresses. So, the stress tensor includes 3 stress components at each point. It is noted that sxi,Sx symbolizes the stress along the x-axis at point i under the unit tectonic stress component Sx and iA(N T).
Massive rock Rock mass containing joint set with high spacing Rock mass containing joint set with low spacing Jointed rock Complex geological structure District with unknown far field stresses
Applied far field stress component
Induced state of stresses at point i
Sx Sy Sxy
Sx ¼ 1 Sy ¼ 1 Sxy ¼ 1 (rgh)x (rgh)y (rgh)z
sxi,Sx sxi,Sy sxi,Sxy sxi,ax sxi,ay sxi,axy sxi,kx sxi,ky sxi,kxy
axrgh ayrgh azrgh b k rgh aþh x b k rgh aþh y b k rgh a þ h xy
ða þb h rghÞx ða þb h rghÞy ða þb h rghÞxy
syi,Sx syi,Sy syi,Sxy syi,ax syi,ay syi,axy syi,kx syi,ky syi,kxy
txyi,Sx txyi,Sy txyi,Sxy txyi,ax txyi,ay txyi,axy txyi,kx txyi,ky txyi,kxy
3.4. Stage 4 In Stage 4, the coefficient matrix must be established. To reach this purpose, the stress states at a certain point should be expressed as the superposition of each field stress component multiplied by the individual stress obtained from numerical modeling based on BEM. For example, the stress component sxi at a certain point i can be written as
sxi ¼ Sx sxi,Sx þSy sxi,Sy þ Sxy sxi,Sxy þ ax sxi, ax þ þ ay sxi, ay þ axy sxi, axy þkx sxi,kx þky sxi,ky þ kxy sxi,kxy
ð4Þ
Assuming the measured stress component at a point i as (sxi)m, we can substitute (sxi)m ¼ sxi in Eq. (2) to obtain: ðsxi Þm ¼ Sx sxi,Sx þ Sy sxi,Sy þ Sxy sxi,Sxy þ ax sxi, ax þ ay sxi, ay þ axy sxi, axy þkx sxi,kx þ ky sxi,ky þkxy sxi,kxy
ð5Þ
where iA(N T). To summarize: 0 B B B B B B B B B B B B B B B B B B B @
sm x1 sm y1 sm x1 ^ ^ ^ ^ ^ ^
sm xy1 0
1 C C C C C C C C C C C C C C C C C C C A
sx1,Sx
Bs B y1,Sx B B txy1,Sx B B ^ B B B ^ B ¼B B ^ B B ^ B B B ^ B B ^ @
1
sx1,Sy sx1,Sxy sx1,ax sx1,ay sx1,axy sx1,kx sx1,kx sx1=x sy1,Sy sy1,Sxy sy1,ax sy1,ay sy1,axy sy1,kx sy1,ky sy1,kxy C C C txy1,Sy txy1,Sxy txy1,ax txy1,ay txy1,axy txy1,kx txy1,ky txy1,kxy C C ^
^
^
^
^
^
^
^
^
^
^
^
^
^
^
^
^ ^
^ ^
^ ^
^ ^
^ ^
^ ^
^ ^
^ ^
^
^
^
^
^
^
^
^
^
^
^
^
^
^
^
^
txyn,Sx txyn,Sy txyn,Sxy txyn,ax txyn,ay txyn,axy txyn,kx txyn,ky txyn,kxy
Table 1 Guidance for selection of numerical methods.
Rock mass properties
Far-field stress component
C C C C C C C C C C C C C C A
0 Rock mass behavior Elastic
Plastic
BEM FEM, BEM DEM –
FDM, FEM FDM, FEM DEM FDM, FEM
FDM, FEM BEM
1 Sx B Sy C B C B C B Sxy C B C B Sx C B C B C B Sy C B C BS C B xy C B C B Sx C B C BS C @ y A Sxy
ð6Þ
Z. Khademian et al. / International Journal of Rock Mechanics & Mining Sciences 55 (2012) 80–85
3.5. Stage 5
To begin, we define
To solve the equations set, since the number of equations is usually more than that of the unknown coefficients {Sx,Sy,Sxy, ax,ay,axy, kx,ky,kxy}, so, a unique solution may not be obtained and a least-squares estimation is, therefore, applied here to find the best solution for set of equations.
g¼
3.6. Stage 6 The purpose of this stage is the determination of the local stress using boundary conditions calculated by the above-mentioned BEM3D back analysis. In previous stages, the effects of different parameters including topography, anisotropy and tectonic stresses are applied on in situ stress variations. Now, we try to exert the effects of heterogeneity of rock mass such as rock properties as well as major and minor geological structures on in situ stresses estimation. It is necessary to model the same area of BEM3D analysis using another numerical method that is able to simulate complex geological structure. According to Table 1, Finite Element (FEM3D) or Finite Difference Method (FDM3D) could be suitable for re-simulation under certain conditions.
3.7. Stage 7 To increase the workability of this algorithm in estimation of in situ stress, it is necessary to compare the measured data (set (N–T)) with those gained from Re-simulation stage (Stage 6).
P
2 2 P ðsci s2 mi Þ ,
i ¼ 1
i ¼ 1
sci
N ¼ p1 ,p2 ,p3 ,:::,pn ,
83
( i A NT g, T N
ð7Þ
where sci and smi are calculated and measured stresses in point i, respectively. According to Eq. (7), if g r0.1, the next stage will be investigated; otherwise, it is addressed to re-calculate a and b coefficients by Eq. (2). Note that parameters a and b indicate a non-linear distribution of certain stresses in terms of depth h and should be determined using a simplex method to minimize the sum of the differences between stress values, given by Eq. (2), and the measured stresses by HF. If the first provision is not satisfied after 10 times of trials and errors, the user is allowed to go to the next step aiming at expediting the calculations. 3.8. Stage 8 As described before, n t numbers of stress measured points are applied to start the algorithm. After determination of local stresses by FEM or FDM (i.e. passing the first provision and finding the final values for a and b), the remaining stress values (t points) should be compared with those obtained from the algorithm (Stage 6) in the same points: P 2 2 1 ðsci smi Þ b ¼ i ¼P , iAT ð8Þ 2 i¼1
sci
To be able to proceed along the algorithm, the second provision must be satisfied, too. If b r0.1, the process is terminated and we will have the values of in situ stresses in every arbitrary point; otherwise, according to the second provision, another set (N–T) of measured points must be chosen from set N to restart the algorithm. Similar to the stopping condition for the first loop, if the second provision is not satisfied after L times, the user should decrease the accuracy of the estimation (i.e. increase the values of g and b). Also, the increase of the number of in situ stress measurements by HF is another solution. The parameter L is given by L¼
n! ðntÞ!
ð9Þ
4. Algorithm applicability To show the applicability of this algorithm, it is applied for Seymare Dam (SD) to determine its stress field using the suggested method and based on the local stresses measured by HF.
Fig. 2. Seymare Dam district and its power tunnel.
Table 3 Input material properties for Hybrid method [16]. Description
Density (t/m3)
Poisson’s ratio
Young’s modulus (GPa)
Shear and normal stiffness (MPa/m) 60
Upper Assmary
1 2 3
2.72 2.78 2.54
0.31 0.32 0.31
1000 950 150
Middle Assmary
4 5
2.71 2.68
0.3 0.28
600 1200
Lower Assmary
6
2.69
0.29
1800
84
Z. Khademian et al. / International Journal of Rock Mechanics & Mining Sciences 55 (2012) 80–85
4.1. Introduction of Seymare district
4.3. Estimation of far-field stress state by BEM3D
Fig. 2 shows the district of SD and Power Tunnel. Local hills approach 200 m–300 m in height. In this district, there are three main geologic units named upper Asmary, middle Asmary and lower Asmary featuring a maximum thickness of 150 m as well as three main normal faults with a strike of N100E, a dip of 60 and an estimated throw of about 30 m. Rock mass geotectonic and joints set properties are also shown in Table 3 [16].
The BEM3D model simulating SD district by Examine3D software is shown in Fig. 3, which is 2.1 km in length and 1.7 km in width. Examine3D is a software based on three-dimensional Boundary Element method which is applied in order to simulate the area of SD site (Stage 2). Also, it is possible to program the BEM3D in Matlab software and analyze the model in Ansys instead of using Examine3D software. The grid points and zones’ numbers of the model are 155,900 and 77,760, respectively. According to Stage 3, the induced stresses at thirteen point (set (N–T)) are calculated and the coefficients matrix are then established (Stage 4). Next, a least-squares approach is utilized to control the distribution of the error and to find the best approximation solution for this over-determined set of equations. Table 5 presents the results calculated by BEM3D for far-field stress state (Stage 5), in which y direction shows the North.
4.2. Local stress measurements using Hydraulic Fracturing Rock stress measurements have been performed at fifteen points in Power Tunnel Trade using Hydraulic Fracturing technique (HF) in SD. Thus n ¼15 and according to Stage 1, t parameter is therefore 1rt r6. To start the algorithm, we consider t ¼2. Table 4 reports the in situ stresses values measured by HF at these points. Note that the tensile stress is taken negative while the compressive stress has a positive sign. Table 4 In situ measured stresses for the SD district. Measured location
sv (Mpa)
sh (Mpa)
sH (Mpa)
sxy (Mpa)
Azimuth sH degree
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15
2.17 2.2 2.26 2.1 2.28 2.31 2.298 2.32 2.343 2.365 2.387 2.409 2.432 2.454 2.46
2.2 2.23 2.2 2.35 2.48 2.29 2.28 3 3.54 2.31 2.4 2.24 2.45 2.57 2.68
3.9 3.87 4.04 3.86 4.11 4.26 4.17 4.22 4.01 4.31 4.56 4.41 4.46 4.51 4.52
4.475 4.485 4.496 4.498 4.515 4.524 4.533 4.22 4.42 4.48 4.6 4.3 4.7 4.6 4.8
43 52 51 45 56 52 44 56 50 43 55 46 55 57 44
4.4. Estimation of local stress state by FEM3D According to Table 1 and considering the existent discontinuities in SD district, the most compatible method for re-simulating, is FEM. In order to illustrate how the far-field stress coefficients can be adopted as input boundary conditions (Stage 6), the Plaxis3D software is applied to re-simulate the same area used in BEM3D analysis. As known, Plaxis3D is the software based on FEM3D to solve complex problems in rock mechanics, and in this simulation. The Plaxis3D model features 2.1 km long, 1.7 km wide and 900 m deep. The model is composed of 280,160 zones, 301,700 grid points and fifteen interfaces comprising seven faults, five layer and three joint sets as well. 4.5. Investigation of the first provision As described before, thirteen stress measurement points are applied to start the algorithm. In this stage, the values of a and b are considered equal to 904.6 m and 812.1 m, respectively, as initial values [1]. The first try results in g Z0.1; therefore, a and b must been modified using a simplex method. The new values of a and b, 890.34 and 786.45, cause g to be decreased to 0.08 and the first provision is therefore satisfied (Stage 7). 4.6. Investigation of the second provision After satisfying first provision, the remaining stress values (t ¼2) measured by HF should be compared with the estimated data in Stage 6. Following Stage 8 and using Eq. (8), b is obtained more than 0.1; thus, other set of measurement points must be selected. In case of SD district this loop repeats twice and finally by choosing t ¼3, b is reduced to 0.098. Table 6 shows the
Table 5 Determination of boundary condition by far field stresses.
Fig. 3. Discretization of the boundaries into triangular leaf elements.
Far field stress component
Coefficient of far field stress component
Calculated far field stress component
Sx Sy Sxy
Sx ¼1.68 Sy ¼3.82 Sxy ¼ 0.695 ax ¼0.778 ay ¼ 0.674 az ¼0.967 kx ¼ 0.667
1.68 3.82 0.695 0.778rgh 0.674rgh 0.967rgh 0:667 b þa h rgh
ky ¼ 0.536
0:536 b þa h rgh
kxy ¼ 0.236
0:236 b þa h rgh
axrgh ayrgh axyrgh kx b þa h rgh ky b þa h rgh kxy b þa h rgh
Z. Khademian et al. / International Journal of Rock Mechanics & Mining Sciences 55 (2012) 80–85
Table 6 Comparison of measured in situ stresses (M) with those calculated (Cal) using the proposed algorithm in SD district (MPa). Location
sh Cal
1 2 2 1.8 3 2.4 4 1.99 5 2.22 6 2.369 7 2.3 8 2.8 9 3.401 10 2.4 11 2.34 12 2 13 2.1 14 2.3 15 2.61 g ¼ 0.08, b ¼ 0.098
sH
sxy
M
Cal
M
Cal
M
2.2 2.235 2.2 2.35 2.48 2.295 2.277 3 3.54 2.318 2.4 2.244 2.458 2.571 2.68
3.6 4 4 3.7 4.039 4.1 3.7 4 3.8 4.1 4.3 4.6 4.5 4.5 4.4
3.9 3.87 4.04 3.863 4.108 4.26 4.174 4.222 4.01 4.318 4.567 4.415 4.463 4.511 4.521
4.117 4.376 4.66 4.201 4.68 4.74 4.6 4.44 4.3 4.35 4.7 4.5 4.6 4.357 4.67
4.475 4.485 4.496 4.498 4.515 4.524 4.533 4.22 4.42 4.48 4.6 4.3 4.7 4.6 4.8
comparison of measured in situ stresses by HF with the finally calculated values by the proposed algorithm in SD. This comparison shows how acceptable the results of the proposed algorithm are in the sense of calculating the in situ stresses.
5. Conclusions A new algorithm for determination of the far-field and local stresses based on stress measurements just in a few regions is developed in this paper. The algorithm considers the tectonic stresses, topography, anisotropy and heterogeneity of the rock mass as effective parameters. The measured in situ stress is decomposed into tectonic and gravitational components, and nine independent far-field stress coefficients are achieved by back analysis using BEM3D simulation and the local stresses are then obtained using FEM or FDM based on the back analysis results. The proposed hybrid method is, also, applied on SD site
85
and the results show how the calculated in situ stresses coincide well with the measured ones. Also, it is shown if the number of stress measurement data is adequate, the stress state will be calculated with the accuracy of 10% in an unknown location. Moreover, it is possible to achieve an estimated in situ stress with desired error by changing the provisions in the algorithm.
References [1] Li G, Mizuta Y, Ishida T, Li H, Nakama S, Sato T. Stress field determination from local stress measurements by numerical modeling. Int J Rock Mech Min Sci 2009;30:111–23. [2] Savage WZ, Swolfs HS, Powers PS. Gravitational stresses in long symmetric ridges and valleys. Int J Rock Mech Min Sci 1985;22:291–302. [3] Hakala M, Kuula H, Hudson JA. Estimating the transversely isotropic elastic intact rock properties for in situ stress measurement data reduction: A case study of the Olkiluoto mica gneiss, Finland. Int J Rock Mech Min Sci 2007;44:14–46. [4] Becker AA. The boundary element method in engineering. London: McGrawHill; 1992. [5] Aklin JE. Finite element analysis with error estimators. Oxford: Elsevier; 2005. [6] Peng S, Zhang J. Engineering geology for underground rocks. New York: Almas Schimmel; 2007. [7] Amadei B, Stephenson O. Rock mechanics for underground mining. London: Chapman & Hall; 1993. [8] Jaeger JC, Cook NGW, Zimmerman RW. Fundamentals of rock mechanics. 4th ed. Oxford: Wiley-Blackwell; 2007. [9] Brown ET, Hoek E. Trends in relationships between measured in situ stresses and depth. Int J Rock Mech Min Sci 1987;15:211–5. [10] Kanagawa T, Hibino S, Ishida T, Hayashi M, Kitahara Y. In situ stress measurements in the Japanese islands-over coring results from a multielement Gauge used at 23 sites. Int J Rock Mech Min Sci 1986;23:29–39. [11] Brady B, Brown ET. Rock mechanics for underground mining. London: Chapman & Hall; 1993. [12] Cornet FH, Valette B. In situ stress determination from Hydraulic Injection test data. J Geophys Res 1984;89:527–37. [13] Craig D, Richard A, Regueiro I, Borja Ronaldo I. Modelling brittle fracture, slip weakening and variable friction in geomaterials with an embedded strong discontinuity Finite Element. US Department of Energy: Sandia National Laboratories; 2006. [14] Muir D. Geotechnical modeling. London: Chapman & Hall; 2004. [15] Amadei B, Savage WZ. Effect of joints on rock mass strength and deformability. In: Hudson JA, editor. Comprehensive Rock Engineering. Oxford: Pergamon Press; 1993. [16] Iran Water & Power Resources Development Co. The in situ stress measurement result in Seymare Dam site by Hydraulic Fracturing. 3rd Report. Iran; 2010.