An analytical solution of equivalent elastic modulus considering confining stress and its variables sensitivity analysis for fractured rock masses

Donghui Chen, Huie Chen, Wen Zhng,*, Junqing Lou, Bo Shn

a College of Construction Engineering, Jilin University, Changchun,130026, China

b Northeast Electric Power Design Institute Co., Ltd., Changchun,130021, China

Keywords:Equivalent elastic modulus Confining stress Rock masses Sensitivity

A B S T R A C T

1. Introduction

As a kind of natural material, rock masses are widely encountered in different engineering projects, such as rock slopes, dam foundations, and tunnels (Wu et al., 2020). The reasonable determination of the mechanical parameters of fractured rock masses is critically necessary for the safe design of engineering projects(Singh and Seshagiri Rao,2005).It is well known that rock masses are composed of intact rocks and various types of discontinuities such as fractures, bedding planes, and faults (Huang et al., 2019;Chen et al.,2020;Nie et al.,2020;Zhang et al.,2016).The existence of discontinuities makes it challenging to determine the macroscopic mechanical parameters of fractured rock masses under laboratory or in situ conditions for rock engineering projects(Laghaei et al., 2018; Yong et al., 2018; Zhang et al., 2020a,b).

Numerous studies have demonstrated that the macroscopic mechanical properties of fractured rock masses in large-scale analyses are continuous or pseudo-continuous (Yang et al., 2014).Thus, the equivalent continuum approach is a feasible means to obtain the mechanical properties of fractured rock masses (Jiang et al., 2017). This approach has the advantage of being more appropriate for representing the global behavior of fractured rock masses in large-scale projects, where the influences of discontinuity systems are implicitly expressed in the equivalent parameters(Min and Jing, 2003). In the continuum approach, the equivalent elastic modulus of fractured rock masses is an indispensable parameter in deformation and stability analyses (Hoek and Diederichs, 2006; Yang et al., 2014; Majdi and Beiki, 2019). Therefore, there is a tremendous need to develop a convenient estimation of the equivalent elastic modulus.

Currently, the equivalent elastic modulus of fractured rock masses is generally considered constant for analyzing deformation and stability. It is rarely reported that the equivalent elastic modulus is altered when the confining stress is considered in practical projects (Asef and Reddish, 2002). Nevertheless,qualitative experimental investigations have indicated that the equivalent elastic modulus generally increases with the increase in confining stress; that is, it is characterized by stress-dependency(Nasseri et al., 2003; Min and Jing, 2004). In the Earth’s crust, the confining stress is referred to as crustal stress.Crustal stress is the fundamental stress environment of natural rock masses, and the stress values are often diverse in different zones,which could lead to a discrepancy in the equivalent elastic modulus. Therefore, an intensive study of the equivalent elastic modulus of rock masses by giving priority to crustal stress(confining stress)is essential in rock mass projects.

Over the years,substantial efforts have been made to determine the equivalent elastic modulus of rock masses. Currently, the available methods for estimating the equivalent elastic modulus fall into two categories, i.e. direct and indirect methods (Kulatilake et al.,1993; Min and Jing, 2003; Khani et al., 2013).

Direct methods include laboratory and in situ tests(Yang et al.,2015;Panthee et al.,2018).There is no doubt that test methods are the most explicit means to measure the elastic modulus.The most common laboratory tests are uniaxial or triaxial compression tests(Laghaei et al., 2018), which are typically performed on collected small-scale rock samples containing micro-fractures (Huang et al.,2019). However, the disturbance of the original stress state of the rock samples,the effect of the rock samples size on the parameters as well as an inadequate representation of discontinuity systems make the results acquired from laboratory tests hardly representative of field conditions(Singh and Seshagiri Rao,2005;Yang et al.,2015;Fattahi and VarmazyariBabanouri,2019).In situ tests mainly include field plate loading and borehole expansion tests (Kayabasi and Gokceoglu, 2018). Although the size of the rock samples in in situ tests is substantially larger than that in laboratory tests, the samples still do not represent the entire area of interest(Goodman,1980). Besides, in situ tests are often time-consuming and costly(Kavur et al., 2015; Fattahi and VarmazyariBabanouri, 2019).Furthermore,in situ tests produce results with large dispersion due to intrinsic hidden fractures and imprecise boundary conditions(Yang et al., 2015; Laghaei et al., 2018).

Because of the above-mentioned difficulties,indirect methods have been of great interest to researchers in recent years.Indirect methods include empirical, numerical and analytical methods(Laghaei et al., 2018). These methods have their special application fields and have gained widespread prevalence in practical rock engineering. Their characteristics and disadvantages have been reviewed and summarized in the literature (Kulatilake et al., 1993; Min and Jing, 2003; Cui et al., 2016; Laghaei et al.,2018).

There has been a strong focus on the analytical method,which is a concise and effective means for obtaining the equivalent elastic modulus of fractured rock masses,because it produces results that highlight the influences of the most crucial issues or variables to solve a problem (Min and Jing, 2003). The endeavors to obtain analytical solutions have a long history, and several analytical solutions have been derived for simple fracture systems (Salamon,1968; Singh, 1973; Gerrard, 1982; Fossum, 1985; Hu and Huang,1993, Li, 2000; Ebadi et al., 2011). Researchers have determined the global elastic modulus of rock masses with continuous persistent fractures by considering the stress-strain relation of the constituent materials (the intact rocks and fractures).This is based on the assumption that the global deformation of fractured rock masses is the sum of the responses of the components based on the elastic theory. Confining stress is an important factor affecting the elastic modulus,because it not only affects the deformation of the fractures, but also the deformation of the intact rocks (i.e. the Poisson’s effect caused by confining stress). Ebadi et al. (2011)proposed an analytical solution for rock masses containing fracture sets considering the effect of confining stress on the fractures,but the effect of confining stress on the intact rocks is still indistinct from the perspective of analytical solution. Consequently, it is necessary to propose an analytical solution of equivalent elastic modulus comprehensively taking into account the influence of confining stress on each component of rock masses.

In this paper, we develop an analytical solution to determine the equivalent elastic modulus of fractured rock masses by considering the influence of confining stress on constituent materials(the intact rocks and fractures).First,the analytical solution of the equivalent elastic modulus of fractured rock masses with random discrete fractures(RDFs)or regular fracture sets(RFSs)is deduced based on equivalent continuum theory (Section 2).Subsequently, the performance of the derived analytical solution is evaluated (Section 3) through the three-dimensional (3D)distinct element code (3DEC). Finally, the sensitivity of the variables in the analytical solution is analyzed (Section 4) using a global sensitivity analysis approach, i.e. extended Fourier amplitude sensitivity test (EFAST).

2. Description of the proposed analytical solution

The equivalent continuum theory holds that the elastic deformation of rock masses is the superposition of the elastic deformation of the intact rocks and fractures. Based on the above principles,the proposed analytical solution of the equivalent elastic modulus of fractured rock masses containing RDF (Fig.1a) or RFS(Fig.1b) are introduced in detail.

2.1. Elastic deformation of the rock masses

In the proposed analytical approach, it is assumed that the persistence of all fractures is infinite and the aperture is zero. A 3D rectangular coordinate system is utilized, where theX-axis points to the east, theY-axis to the north and theZ-axis to the vertical direction. First, a rock mass containing one fracture is chosen as a primary model to interpret the detailed derivation of the analytical solution of the equivalent elastic modulus, as shown in Fig. 2.

Letni(cosαi,cosβi,cosγi)stands for the normal vector of theith fracture plane, where αi, βiand γidenote the angles between the normal vector and the axesX,YandZ,respectively.Then,the angles αi,βiand γiare expressed with respect to the dip angle(φi)and dip direction (λi) of theith fracture as follows:The increments of the applied confining stresses are Δσ1, Δσ2and Δσ3, respectively, as shown in Fig. 2. The extensively applied conventional triaxial compression test is suitable to investigate the influence of the confining stress on the elastic modulus.Hence,let the confining stress Δσ2be equal to Δσ3. The components of the total stress of theith fracture plane in theX-,Y-andZ-directions are obtained based on the elastic theory:

Fig. 2. The 3D rectangular coordinate system XYZ.

Then, the increments of the normal (Δσni) and shear (Δτi)stresses on theith fracture plane can be expressed as

Fig.1. Rock mass samples containing (a) random discrete fractures and (b) regular fracture sets.

The increments of the normal (Δdni) and shear (Δdsi) displacements on this fracture under loading,Δσ1and Δσ2,respectively,are obtained as

whereKniandKsiare the normal and shear stiffnesses of theith fracture, respectively.

Finally, the displacement increment of theith fracture in the applied stress(Δσ1) direction,Δdi, is obtained as

Assuming that the intact rocks are characteristized by elastic and isotropic properties,the elastic displacement increment of the intact rocks is obtained based on the generalized Hooke’s law:

where ν andErrepresent the Poisson’s ratio and the elastic modulus of the intact rocks,respectively;andLrepresents the height of rock masses sample, as shown in Fig. 2. The total elastic deformation(Δd)of a rock mass containing one fracture is the superposition of the deformation of the intact rocks and the deformation of one fracture, i.e.

2.2. Equivalent elastic modulus of the rock masses

2.2.1. Rock masses containing RDF

For the rock masses containing RDF (Fig. 1a), the total elastic deformation (Δd) can be expressed as

whereMstands for the number of fractures within rock masses.The corresponding strain of rock masses is defined as

2.2.2. Rock masses containing RFS

For the rock masses containing RFS (Fig. 1b), the total elastic deformation(Δd) can be expressed as

whereNis the number of fracture sets;andn(=Lcos γi/Si)is the number of fractures of theith fracture set in the length ofL, as shown in Fig. 2. Then, the corresponding strain of rock masses is defined as

Finally, the equivalent elastic modulus of rock masses containing RFS is obtained as

Therefore, the elastic modulus of rock masses in the applied stress(Δσ1) direction,Em, can be obtained as follows:

Here,we introduce the variable stress ratio(k),which is defined as the ratio of Δσ2to Δσ1(i.e.k=Δσ2/Δσ1). The range ofkis from 0 to 1. Then, Eq. (12) can be simplified as follows:

“No; do not do that,” replied the nightingale; “the bird did very well as long as it could. Keep it here still. I cannot live in the palace, and build my nest; but let me come when I like. I will sit on a bough outside your window, in the evening, and sing to you, so that you may be happy, and have thoughts full of joy. I will sing to you of those who are happy, and those who suffer; of the good and the evil, who are hidden around you. The little singing bird flies far from you and your court to the home of the fisherman and the peasant’s cot. I love your heart better than your crown; and yet something holy lingers round that also. I will come, I will sing to you; but you must promise me one thing.”

For the rock masses containing multiple fractures or fracture sets,the stress state on the fracture planes may be disturbed due to the interactions between the intact rocks and fractures when subjected to an external load. In the current proposed model, it should be noted that interactions and disturbance on stress distribution are not considered like all analytical methods for fractured rock masses (Min and Jing, 2003).

For the case without considering the confining stress, i.e.k= 0, the analytical solution of the equivalent elastic modulus of rock masses containing two fracture sets was derived by Yoshinaka and Yamabe (1986) and Li (2000), as shown in Eq.(17). This equation has been extensively applied in the literature(Park and Min, 2015; Cui et al., 2016; Zoorabadi, 2016; Alshkane et al., 2017).

LetN=2 andk=0,Eq.(16)proposed in this study is identical to Eq. (17).

3. Verification of the proposed analytical solution

The equivalent elastic modulus of fractured rock masses containing one, two and three fracture sets are obtained using numerical tests to verify the performance of the proposed analytical solution in this study.

3.1. Establishment of the numerical models

In the proposed analytical solution, the equivalent elastic modulus of fractured rock masses is determined based on the geometrical data of the fractures as well as the elastic parameters of the intact rocks and the fractures. Therefore, the 3DEC, which can well represent the variables in the analytical solution,is utilized for the numerical tests(Itasca, 2018).

A compression test is the most common method to determine the elastic modulus. Because of equal confining stress in the analytical solution, the conventional triaxial compression test is adopted here.Most fracture measurements in rock masses present a mixture of multiple sets(Hudson and Priest,1983;Liu et al.,2017).Here, three types of numerical models of fractured rock masses containing one, two and three fracture sets are established,respectively,as shown in Fig.3.The geometrical parameters of the fracture sets are listed in Table 1.

The stress-controlled method is a practical way for applying loads in a numerical model of the triaxial compression test(Alshkane et al., 2017). Therefore, the axial load is applied using a stress-controlled method at the upper boundary, while the lower boundary is fixed in the vertical direction.In numerical models,the confining stress is applied on both sides of the model (Fig. 3d).

The axial displacements are monitored by monitoring grid points on the upper face, as shown in Fig. 3d, and the axial strains are calculated by averaging the grid strains using a FISH script.Note that FISHisaprogramminglanguageembeddedin3DECthatenablesusers to define new variables and functions.The axial stress is measured by monitoring the average zone stress σzzfor all zones using a FISH script.

Fig. 3. Configuration of the numerical models of fractured rock masses containing (a)one fracture set,(b)two fracture sets,and(c)three fracture sets.(d)Loading scheme in compression test and location of the measurement points (a fractured rock mass containing one fracture set is used as an example).

A linear elastic constitutive model is selected to represent the elastic deformation behavior of the intact rocks and fractures. The hypothetical geometrical and mechanical properties of the intact rocks and fracture sets are used as an example for the verification,and are listed in Table 1. These parameters can be obtained using existing methods and techniques. Based on the fracture data measured in the field,clustering methods proposed by scholars are capable of grouping fractures sets,and obtain the average dip angle and dip direction of each set(Shanley and Mahtab,1976;Hammah and Curran,1998; Jimenez, 2008). The spacing of each set can be determined by the scanline method(Priest and Hudson,1981).The deformation parameters listed in Table 1 can be obtained by tests or empirical relations, and the specific procedures can be referred to Ulusay(2015) and Yoshinaka and Yamabe (1986).

Table 1 The geometrical and deformation parameters of fracture sets.

3.2. Comparison of the results

A series of conventional triaxial compression tests is performed under different confining stresses according to the cases shown in Fig. 3. The slope of an approximately straight line of the stressstrain curve is used as the magnitude of the equivalent elastic modulus.

The results of the analytical solution and numerical tests are presented in Fig.4.The equivalent elastic modulus of fractured rock masses nonlinearly increases with an increase in the stress ratiok,i.e.the confining stress.The results demonstrate that the equivalent elastic modulus is characterized by stress-dependency.

Fig.4. Variation of the equivalent elastic modulus of fractured rock masses with stress ratio k: (a) Rock masses containing one fracture set, (b) Rock masses containing two fracture sets, and (c) Rock masses containing three fracture sets.

The comparison of the results indicates that the relative errors between the analytical solution and the numerical tests are approximately 5%,representing a good agreement.The interactions between the fractures and the intact rocks separated by the fractures cannot be taken into account in an analytical solution (Cui et al., 2016). However, the so-called intersection generally results in local stress concentration(Min and Jing,2003),which may lead to a slight deviation between the analytical solution and numerical results to some extent.

4. Sensitivity analysis

The sensitivity analysis is capable of quantifying the impact of each input variable on the model output (Varella et al., 2010).Deformation analysis is an essential step in practical engineering projects, and the accuracy of the analysis results depends on the selected mechanical parameters (Majdi and Beiki, 2019). It can be concluded that multiple variables influence the equivalent elastic modulus of fractured rock masses. The structure of natural rock masses is invariably affected by the complex geological environment, resulting in significant spatial variability of the variables(Majdi and Beiki, 2019). As a result, a sensitivity analysis is indispensable to identify the dominant variables for future applications in engineering projects.

4.1. EFAST

The EFAST is a quantitative and global sensitivity analysis algorithm proposed by Saltelli et al. (2000) that combines the advantages of the Fourier amplitude sensitivity test and Sobol’s algorithm. The EFAST is suitable for complex nonlinear and nonmonotonic models (Saltelli et al., 2010), such as the proposed analytical solution of the equivalent elastic modulus.Therefore,the EFAST is employed in this study to evaluate the sensitivity of the variables in the analytical solution.

Based on the variance decomposition, the EFAST approach provides two types of sensitivity indices for each variable:the firstorder sensitivity index,i.e.the main sensitivity index(MSI),and the total sensitivity index (TSI). The MSI represents the influence of a single variable on the model output. In contrast, the TSI measures the overall contribution of each variable to the model output considering the interactions between the variable of interest and all other variables.This method takes into account the coupling effect among the variables and is suitable for practical problems. The following is a brief introduction to the EFAST.

A model withnvariables can be written as

wherexi(i= 1, 2, …,n) represents the input variable. The total varianceVYof the model output is decomposed into the following components using the variance decomposition of Sobol’s method:

whereViis the variance of input variablexi;Vijis the variance of the interaction between variablesxiandxj;VijmandV12…nare the higher-order variances of the interactions between multiple variablesxi,xj,xm, …,xn. Then, the MSI for the input variablexiis defined as

Similarly, the second and higher-order sensitivity indices are defined as

The TSI of variablexiis the sum of each order of the sensitivity indices:SiandSTirange between 0 and 1, with higher values indicating more essential effects. If there are no interactions among the variables in the model,the other terms equal zero except forSi,i.e.Si=STi. The detailed process will be introduced in Section 4.2. For a more detailed description of the EFAST algorithm, readers are referred to the literature (Saltelli et al.,1999, 2004).

4.2. Design of the experiments for sensitivity analysis

As mentioned above, fractures within rock masses generally occur in several sets. Therefore, a sensitivity analysis of the analytical solution of the equivalent elastic modulus of fractured rock masses with fracture sets is conducted. We take the rock masses with one fracture set as an example to conduct the sensitivity analysis. Since the interactions between fractures in the analytical solution are not considered,the same sensitivity analysis can be obtained for rock masses with multiple fracture sets.

Several specialized tools for conducting sensitivity analysis are available. In the current study, the software package SimLab2.2 published by the Joint Research Centre of the European Commission is employed. The program is freely available and efficiently implements the EFAST algorithm. Many models have complicated equations that can hardly be coded using simple mathematical functions in SimLab,thus SimLab allows users to link to an external model via executable files.The flowchart of the sensitivity analysis is shown in Fig.5.The experimental design consists of the following five steps:

Fig. 5. The flowchart of the sensitivity analysis.

(1) Step 1: Selection of variables. As shown in Eq. (16), the contributions of seven variables(i.e.k,ν,Er,S,Kn,Ksand γ)in the analytical solution need to be analyzed.

(2) Step 2:Selection of the distribution types and ranges of each variable. Determining each variable’s distribution types and ranges is necessary, which might influence the sensitivity indices in a global sensitivity analysis (Helton, 1993; Wang et al., 2013). The selection of the distributions and ranges of the variables depends on the research objectives and the available information (Helton, 1993). First, uniform distribution is assigned to seven variables to easily determine the lower and upper boundary values of the distribution when the initial information on the distribution types is limited(Vanuytrecht et al., 2014; Vazquez-Cruz et al., 2014). Subsequently, normal and exponential distributions are assigned to all variables to explore the influence of the distribution types on the sensitivity analysis results.

Generally, there are two feasible approaches to determine the ranges of the variables: (i) The values of variables with explicit physical definition should cover the entire physical range as much as possible; and (ii) For variables with no physical definition, researchers should refer to the literature or expert experience (Ma et al., 2015). In this study, all variables in Eq. (16) have explicit physical meaning. The ranges of the seven variables are summarized in Table 2,which represent typical conditions.Especially,the spacing of the fractures is theoretically arbitrary. The objective of this study is to determine the equivalent elastic modulus of fractured rock masses. To date, no uniform standard has been developed for the spacing ranges in which fractured rock masses can be equivalent to the continuous medium. Therefore, we select a conservative range of less than 2 m according to the classification by the ISRM (1978).

However,the ranges of the variables may be varied for different geological conditions. Thus, experiments are conducted to investigate the influences of variable ranges on the sensitivity indices.First,the ranges of all variables are set to the default values listed in Table 2 to observe the MSI and TSI.Then,only the dip angle ranges are separated into 0°-30°, 30°-60°and 60°-90°, and the others are consistent with Table 2. Subsequently, only the spacing ranges are set to 0-0.67 m,0.67-1.34 m and 1.34-2 m,and the others are consistent with Table 2.

Table 2 Variables and statistical settings used for sensitivity analysis.

(3) Step 3: Generation of input samples. The EFAST algorithm needs to sample for each variable. It is suggested to be effective when the sample size of each variable is 65 times larger than the number of variables(Luan et al.,2017;Marino et al.,2008).This study generates a sample of sizeN=2002(286 × 7) for each variable based on the specified distribution types in Step 2 to ensure that the sensitivity indices reach a stable value.

(4) Step 4: Calculation of model output. A total of 2002 calculations to obtain the equivalent elastic modulus are carried out using Eq. (16)and the sample values obtained in Step 3.

(5) Step 5:Analysis of the model outputs.The results obtained in Step 4 are imported into SimLab in the format specified by SimLab. Finally, both the first order and total sensitivity indices(SiandSTi)are calculated,as described in Section 4.1.

4.3. Results of sensitivity analysis

The sensitivity analysis results obtained by EFAST for the default distribution type (i.e. uniform distribution) and ranges of the variables are presented in Fig. 6a. The sensitivity indices exhibit marked discrepancies. The rock masses can be regarded as the assemblage of the intact rocks separated by fractures. When subjected to external compression loads,the intact rocks play the most critical role in the elastic compression. Therefore, the elastic modulus of the intact rocks is the most sensitive variable, and its MSI and TSI are 0.28 and 0.49, respectively. The dip angle and spacing are critical geometric parameters of the fractures that change the stress state of the fracture planes, leading to the deterioration of the mechanical properties of rock masses to some extent. Thus, the dip angle and spacing are the second and third most influential variables, respectively. The other variables have negligible or minor influence on the equivalent elastic modulus.The distinct discrepancies in the sensitivity indices of the variables account for the different contributions of each variable to the equivalent elastic modulus.Fig.6a indicates that the TSI values are higher than that of the MSI.This means that the interactions of the seven variables on the model output are much more significant than the effect of the individual variable.

Fig.6. Sensitivity indices of the variables with different distribution types:(a)Uniform distribution, (b) Normal distribution, and (c) Exponential distribution.

Fig. 6b and c present the sensitivity analysis results for the variables with the normal and exponential distributions, respectively.A comparison of Fig.6a-c shows that the MSI and TSI of each variable are slightly different, but the sensitivity rankings of the variables with varying types of distribution are identical. These investigations demonstrate that the sensitivity indices generally depend less on the distributions types of the variables, consistent with previous findings in other research fields (Helton, 1993;Vanuytrecht et al., 2014). In addition, it also proves that the assumption of uniform distribution is justified in this study.

An analysis of the sensitivity indices of a variable involves the estimation of the range of the variable for a specific background.The range of a variable can be separated into several finite intervals with identical width, and some available information in terms of the influence of variable ranges on the sensitivity indices can be acquired.The results of the sensitivity indices for different ranges of the dip angle and spacing are presented in Fig. 7, which indicate significant differences in the sensitivity indices and rankings of the variables. As shown in Fig. 7a and b, the dip angle with a smaller range results in smaller MSI and TSI compared to the results for the default range (Fig. 6). Besides, the MSI and TSI increase with an increase in the dip angle,which indicates that the equivalent elastic modulus is more sensitive to a steep dip angle than a gentle dip angle. The increase in the dip angle leads to a decrease in the sensitivity to spacing and normal stiffness and an increase in the sensitivity to the stress ratio.When the spacing range is separated into smaller ranges, the MSI and TSI of spacing also decrease compared to the default case.Furthermore,it is observed in Fig.7c and d that the equivalent elastic modulus is more sensitive to close spacing than wide spacing. The findings mentioned above show that the variable range significantly affects both the sensitivity indices and their relative importance.

In summary, the sensitivity analysis demonstrates that the sensitivity indices of the variables are more dependent on the ranges than the distribution types of the variables. In practical projects,the assumption of uniform distribution for the variables is appropriate due to the convenience of determining the upper and lower boundaries in a global sensitivity analysis.

Fig.7. Sensitivity indices values of the variables:(a)MSI for different ranges of dip angle,(b)TSI for different ranges of dip angle,(c)MSI for different ranges of spacing,and(d)TSI for different ranges of spacing.

5. Discussion

The equivalent elastic modulus is a primary input parameter for the deformation and stability analyses of fractured rock masses in the equivalent continuous approach (Yang et al., 2014). As described in Section 3, the equivalent elastic modulus is characterized by stress-dependency (i.e. dependency of the modulus on the confining stress). The confining stress around natural rock masses exists in the form of crustal stress. In practical large-scale rock engineering projects, crustal stress (confining stress) is inevitable due to the combined effect of gravitational stress and tectonic stress (Zhang et al., 2017). Therefore, the discrepancy in deformationcaused by the difference in the equivalent elastic modulus is analyzed by comparing the cases with and without considering confining stress.

Here, a deep rock tunnel excavation project is taken as a hypothetical example to illustrate the discrepancy in deformation quantitatively.We assume that the rock masses around the tunnel contain two fracture sets (i.e. Sets 1 and 2 listed in Table 1). The equivalent continuous approach is adopted to analyze the deformation caused by the excavation of the tunnel. The dimensions of the established numerical model of the deep tunnel based on FLAC3D are shown in Fig. 8a. The initial crustal stress field is specified according to the case used by Xing et al.(2019);that is,the stress ratio is set to 1 for the two horizontal stresses. The case where confining stress is not considered is marked as Case 1, and the case considering the confining stress is marked as Case 2.

The hypothetical geometrical and deformation parameters of the fracture sets(Sets 1 and 2)and rock blocks are shown in Table 1.The equivalent elastic moduli of Cases 1 and 2 are calculated by assigningk= 0 and 1 to Eq. (16), whose values are 15.5 GPa and 20.3 GPa, respectively. A linear elastic constitutive model is employed. The sides and bottom boundaries of the model are constrained,and the top boundary is free.As shown in Fig.8a,four monitoring points(A,B,CandD)around the tunnel are installed to obtain the deformation quantitatively. PointsAandBare used to monitor the horizontal displacement, and pointsCandDare used to monitor the vertical displacement.

The deformation results of the four monitoring points for the two cases after the tunnel excavation are presented in Fig.8b.Note that only the magnitude of the deformation value is of interest here.The deformation of the tunnel roof caused by excavation is the largest, which should be paid more attention to in practical engineering. Moreover, it is found that the displacement values of all monitoring points in Case 2 are approximately 30%lower than that in Case 1. The reason is that the equivalent elastic modulus increases when the confining stress is considered. The induced discrepancy in the deformation due to confining stress is manifest and cannot be neglected in practical deformation and stability analyses.

Fig. 8. A rock engineering project about the tunnel excavation: (a) Numerical model of the deep tunnel, and (b) Comparison of the results of vertical displacement caused by the excavation of the tunnel for two cases.

In large-scale engineering projects, the equivalent continuum approach is an alternative means for the advantage of implicitly incorporating the influence of discontinuities into the equivalent mechanical behavior (Wu and Kulatilake, 2012). The abovementioned simple example quantitatively illustrates the discrepancy in the deformation resulting from the effect of the initial crustal stress (confining stress) on the equivalent elastic modulus,given that the fundamental information on the intact rocks and fractures is available in engineering projects.For practical projects,the application of the proposed analytical solution of the equivalent elastic modulus considering confining stress is summarized as the following four steps:

(1) The initial crustal stress field is of great importance for engineering construction and should be estimated first by combining measured data and inversion approaches.

(2) The stress ratiokis calculated based on the obtained crustal stress field, and then the analysis domains can be divided into different zones according tok.

(3) The equivalent elastic modulus can be predicted for different zones using the proposed analytical solution.

(4) The equivalent elastic modulus is assigned to the corresponding zones for subsequent deformation and stability analyses.

This work provides a detailed description of an analytical method to determine the equivalent elastic modulus of fractured rock masses while considering the confining stress for the case of σ2= σ3. In terms of current technology, the conventional triaxial compression test is commonly used to describe the effect of the confining stress(Asef and Reddish,2002;Ulusay,2015;Yang,2016;Laghaei et al.,2018).However,the development of new equipment for a true triaxial compression test(i.e.σ2≠σ3)has always been a prevalent topic in the field of rock mechanics to describe the influence of the confining stress objectively.Therefore,the analytical solution for the equivalent elastic modulus of fractured rock masses containing RDF or RFS for the case of σ2≠σ3is given in Eqs. (23)and (24), which can provide a reference for future research.Limited by the length of the article,the similar derivation process is not described in detail here.

wherek1= Δσ2/Δσ1,k2= Δσ3/Δσ1,A=k22sin2αcos2α+k12sin 2 βcos 2 β+ sin 2 γcos 2 γ- 2k1k2cos 2 αcos 2 β- 2k2∙cos2αcos2γ- 2k1cos2βcos2γ,B=k2cos2α+k1cos2β+

cos2γ,and the other symbols are the same as in Section 2.

6. Conclusions

In this study, a comprehensive description of an analytical approach for obtaining the equivalent elastic modulus of fractured rock masses was presented. Specifically, the confining stress (i.e. crustal stress) is considered in the analytical solution to better understand its influence on the equivalent elastic modulus. The sensitivity of the variables in the analytical solution was also analyzed. The main conclusions are summarized as follows:

(1) An analytical solution to determine the equivalent elastic modulus of fractured rock masses containing RFS or RDF was proposed based on the equivalent continuum theory. The analytical solution is not only associated with the geometrical data of the fracture systems but also the deformation parameters of the intact rocks and fractures and the confining stress. It has the advantage of conveniently obtaining the equivalent elastic modulus in applying the equivalent continuum approach.

(2) The stress-dependency(i.e. dependency of the modulus on the confining stresses) of the equivalent elastic modulus was represented in the form of an analytical solution. The results indicated that the confining stress had a notable influence on the equivalent elastic modulus, which nonlinearly increased with an increase in the confining stress.Therefore, applications for large-scale engineering projects in various crustal stress environments should fully consider the effect of crustal stress(confining stress)on deformation and stability.

(3) The ranges of the variables affect both the sensitivity indices and their relative importance. A wider range of a given variable will result in a more extensive sensitivity index. Variables that possess the interval of equal width but span different ranges will bring about different sensitivity indices.Therefore,the determination of the ranges of the variables is of great significance before sensitivity analysis. In turn, the sensitivity analysis results can provide a basis for estimating the variable values in the inversion analysis in practical engineering projects.

Declaration of competing interest

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.

Acknowledgments

This research was financially supported by the National Nature Science Foundation of China(Grant Nos.42022053 and 41877220).