# Parametric Evaluation of Simultaneous Effects of Damaged Zone Parameters and Rock Strength Properties on GRC

Authors: Ali Ghorbani; Hadi Hasanzadehshooiili; Amin Eslami
Journal: Advances in Civil Engineering
 (2021)

---

## Abstract

The convergence-confinement method via the ground reaction curve (GRC) is used as the common practice of tunnel design which demands accurate determination of the stress state and material strength behavior in different zones around the tunnel section. Besides, formation of the excavation/blast-induced damaged zone (EDZ/BDZ) adds more complexity to the problem due to variation of elasticity modulus of the rock mass in this zone. As a result, advanced numerical methods via finite element/difference commercial packages or user-coded, semi-numerical techniques are required to develop the GRC, which demands a high degree of proficiency and knowledge of computational plasticity and geomechanics. In this study, a new, simple, and accurate method is proposed for prediction of GRC of circular tunnels constructed in the damaged, elastoplastic rock masses obeying softening in the plastic zone. The effects of deterioration caused by the drilling/blast in the EDZ were taken into account by assuming a reduced and varying Young’s modulus using the disturbance factor, in the form of Hoek–Brown failure criterion and the Geological Strength Index (GSI). Besides, effects of intermediate principal stress and the exponential decaying dilation parameter are taken into account thanks to adoption of the unified strength criterion (USC) as the material strength criteria. To do so, genetic algorithm (GA) via the method of evolutionary polynomial regression (EPR) is used to find a relationship between a number of 19 affecting parameters on the GRC as the input, and the internal support pressure as the target of prediction. Verification analysis was performed to verify the validity of the results using field measurements data as well as other advanced numerical studies found in the literature. Lastly, variation of the support pressure with simultaneous changes in the affecting input parameters was investigated using multivariable parametric study.


---

## Body

## 1. Introduction

The ground reaction curve (GRC) in the framework of convergence-confinement method (CCM) can be assumed as the most common method of tunnel design due to its simplicity while keeping the accuracy in engineering practice. The main advantages of this method lie in turning the 3-dimensional problem of tunnel’s face advancement into a 2-dimensional problem by relating the distance from tunnel face to the inner pressure in GRC. Besides, a factor of safety can be estimated using the ratio between the support capacity and the load demand. The calculations can be made using closed-form solutions, numerical methods, or some combinations of each to develop the appropriate curves. However, there are some limitations in using the GRC which mainly stem from its assumptions: homogenous rock mass, isotropic stress field, circular cross section of the tunnel, and plain-strain condition [1]. In order to accurately calculate the appropriate time for the installation of support system of underground tunnels, also, to restrict the occurring deformations, the ground reaction curve (GRC) of underground openings should be correctly characterized. Development of GRC is conducted based on the required support pressure corresponding to the formation of an arbitrary convergence value. Hence, the state of stress and the material strength play a key role in the determination of such curve [2–13]. In this framework, the behavior of tunnels constructed in the rock mass is generally assessed based on the elastoplasticity theory and the principles of convergence-confinement method [5–7, 11–14]. There are different strength criteria for the rock mass material, such as the Mohr–Coulomb [3, 8, 15, 16], the Hoek–Brown [2, 7, 10, 17], and the unified strength criterion [18–22]. To make the calculations even more accurate, the softening behavior is also taken into the consideration accounting for the variation of material strength parameters from the peak strength parameters to the residual ones (referring to the distance from the tunnel center) [4, 5, 8–10, 13, 23].In addition of the formation of a softening zone around the tunnels, depending on the amount of blasting/excavation damage induced to the surrounding rock mass, a damaged zone will be formed. True characterization of this blast/excavated induced damage zone (BIDZ/EIDZ) is carried out based on the deterioration of material strength parameters, consideration of the variation of Young’s modulus, and the weight of this zone into the problem. The effect of EDZ has been studied numerously [24–36]. However, there is limited, available solutions, which consider the effects of both softening material behavior in the plastic zone and the existence of BIDZ/EIDZ around the tunnel [26, 37–40]. To this end, complicated finite element/difference commercial packages should be adopted. User-coded, semi-numerical techniques are other possible solution approaches to the problem [4, 5, 9, 12, 15–17, 25, 37, 41–45]. Described methods need professional numerical modeling proficiency and the knowledge of computational plasticity and geomechanics. On the other hand, design engineers need a simple, efficient, and timely method for accurate prediction of GRC to make a preliminary design for their underground tunneling projects. Hence, the availability of an accurate, timely, and simple method will be a significant aid to the practitioners and design engineers.As described, the objective of the research is to propose a new, simple, and accurate method for the prediction of GRC of circular tunnels constructed in the damaged, elastoplastic rock masses obeying softening in the plastic zone. It should be described that the adoption of appropriate material strength criteria and dilation functions are two other important parameters affecting on the developed curves. In this study, the unified strength criterion (USC) (selected to consider the effects of intermediate principal stress (b) on the strength of rock mass) and an exponential decaying dilation parameter for showing a better performance are used for the modeling purpose.The proposed method is based on the development of an evolutionary polynomial regression (EPR) technique, which uses the best features of the genetic algorithm (GA) and the least squared method to predict new target functions based on their previous learning of the relationship between input and output parameters. In this study, as the first step, the results of previous solutions for the GRC development were gathered to generate a database using the algorithm proposed by Ghorbani and Hasanzadehshooiili [26]. Material strength parameters in elastic, plastic, and damaged states, the radius of tunnel, and all other affecting parameters along with arbitrary convergence values were considered as the input parameters. The corresponding internal support pressure, which is the target of predictions, is the output parameter. Hence, first, the problem was solved and validated for different strength and geometry parameters. As the second step, an evolutionary polynomial regression (EPR) code, available in the MATLAB environment, was used to train the neural network based on the gathered database, in order to find a mathematical model for forecasting the internal support pressure (Pi) values.In this study, the EPR models are proposed based on the well-validated solutions of the problem of GRC development for tunnels constructed in the damaged, elastoplastic rock mass. The accuracy of the proposed model is shown through the high value of coefficient of determination (COD) and the acceptable value for the root mean squared error (RMSE). The developed model is then applied as the basis for the complementary studies. Sensitivity analysis based on the cosine amplitude method (CAM) has then shown the strength of relationship between the internal support pressure and concerning input parameters. Besides, variation of the support pressure with the simultaneous changes in the affecting input parameters is investigated using the multivariable parametric study.The proposed evolutionary-based approach is introduced as a capable, timely, and accurate method for the prediction of the internal support pressures for the tunnels constructed in the BID/EID, elastoplastic rock mass obeying the USC criterion and showing the softening behavior in the plastic zone. In addition, the importance of material strength parameters, damaged zone’s factors (D), and the tunnel’s radius and their roles in the variation of support pressure is fully characterized using both sensitivity and a three-dimensional, multivariable, graphical parametric study on the input parameters.


## 2. Problem Definition and Methodology

A full description of the proposed algorithm along with the formulation and governing equations is presented in Ghorbani and Hasanzadehshooiili [26]. However, in this part, a summary of the procedure of stress-strain formulation, required to derive the GRC, is presented. Figure 1 illustrates definition of the problem’s geometry. The tunnel has a circular section with radius of bi, circumscribed by an isotropic elastic-plastic rock mass including the excavation/blast-induced damage zone (EDZ/BDZ). The tunnel section is under hydrostatic pressure of s0 while its boundaries are under pressure of the internal support pressure (Pi). The radius of the plastic zone and that of the excavation damaged zone are represented by Rp and REDZ, respectively. Decreasing of the internal support pressure (Pi) results in initiation of the radial displacements. The deformations can be of either elastic or plastic type, distinguished by the critical value of the internal support pressure (Pic). The deformations are considered as elastic, provided that they are formed in a pressure higher than the critical value of the support pressure (Pic) [9, 10]. More decreasing of the internal support pressure leads to formation of the plastic stresses and strains, which are initiated with softening behavior and continue to a residual one. The transition between softening and residual behavior is governed using a parameter called the critical deviatoric plastic strain parameter (γp∗) [2, 4, 9, 10, 17, 25, 46].Figure 1
Problem definition.On the other hand, adoption of a proper plastic potential function is crucial for the progress of the plastic strains. In addition, a flow rule, which can be of either the associated or nonassociated flow rules, governs the material’s plastic behavior in the plastic zone. In fact, considering an associated or nonassociated flow rule can significantly affect the material’s plastic behavior in the plastic zone as the relationship between radial and tangential plastic strain increments is defined by the flow rule. In this paper, the nonassociated flow rule is put into practice to model the relationship between radial and tangential plastic strain increments.Determination of the Pic is dependent on the utilized strength criteria, such that it can be obtained by a closed-form solution in the case of linear Mohr–Coulomb strength criterion [5] or by numerical solution via the Newton–Raphson algorithm in the case of Hoek–Brown criterion [9, 10, 47] as well as the unified strength criterion (USC) [26].In this study, the unified strength criterion (USC) is put into practice to obtain thePic using the Newton–Raphson method. Adoption of USC was due to its advantage in tracking the principal stress-strain incremental paths in a nonlinear way, which yields more accurate results in prediction of the rock mass behavior [2, 4, 10, 47, 48]. Moreover, using USC makes it possible to take the effect of the intermediate principal stress (b), which is believed to have a significant effect on the radius of the EDZ [19–21].The EDZ is developed by approaching the tunnel’s face. The radial stresses within the plastic-EDZ boundary are found numerically using the Runge-Kutta-Fehlberg (RKF) method [49] which also governs the stages of the solution. The problem is considered in the plastic zone provided that the radial stress is larger than the plastic-EDZ boundary’s radial stress, (srP−E). However, for the radial stress values lower than the (srP−E), the output of the last step is used as the initial condition to obtain a new solution for the EDZ.Calculation of the stresses and strains in different points around the tunnel boundary is made using different methods depending on which zone/boundary they are located. Figure2 presents the different zones formed around the tunnel section area. In this problem, there are three main types of boundaries with different calculation approaches:(1)
The elastic-plastic boundary (=outer plastic boundary) which is represented by the radius of plastic zone (Rp).(2)
The outer boundary of the EDZ (=inner plastic boundary).(3)
The inner boundary of the EDZ (=tunnel surface).Figure 2
Stresses at different regions and boundaries, after Ghorbani and Hasanzadehshooiili [26].In order to calculate the stresses and strains of the points located in the elastic-plastic boundary, the Newton–Raphson (NR) method of root finding is used to calculate the radial stress (σr = sR), which is also equal to the critical support pressure (Pic). The stresses located on the outer boundary of the EDZ are calculated using the Runge-Kutta-Fehlberg (RKF) method. On the elastic-plastic boundary, i.e., the external boundary of the plastic zone, the radial stress (sr) is equal to sR and is also equal to the critical support pressure (Pic). Besides, on the internal plastic boundary, i.e., the plastic-EDZ boundary, the stresses are calculated using the Runge-Kutta-Fehlberg (RKF) method. Lastly, the radial stress on the inner boundary of the EZD is equal to the tunnel’s support pressure (Pi).The stresses are first calculated on the boundaries to be used as the initial values for the rest of areas between the boundaries. Lee and Pietruszczak [9] calculated the radial stress within the plastic zone area assuming a linear variation for the radial stress from its initial value at r = Rp, sr = sR to its value on the tunnel surface r = bi, sr = Pi, considering a rock mass with elastic-plastic material.In this study, in summary, the following procedure has been taken for analysis of the plastic zone: first, the Runge-Kutta-Fehlberg method is put into practice to determine the radial stress on the on the plastic-EDZ boundary (srP−E). The region is separated into two zones, namely, the plastic zone and the EDZ. The stresses within the plastic zone can be calculated using the predetermined stresses on the boundaries as the priori, i.e., sR (the radial stresses on the elastic-plastic boundary) and srP−E (the radial stress on the plastic-EDZ boundary). Besides, the stresses and strains within the EDZ are calculated using the known radial stress on the plastic-EDZ (srP−E) and the tunnel surface (Pi). In order to perform the calculations, the plastic and EDZ are divided into m and n number of zones, respectively. Thickness of each annulus is determined so that the corresponding governing equation of equilibrium is satisfied.The most important distinction of consideration of the EDZ into the problem is to assume a reduced and varying Young’s modulus for this zone rather than a constant one, to take the effects of gradual deterioration caused by the drilling/blast around the tunnel. Herein, the relationship by Hoek and Diederichs [50] was used to calculate Young’s modulus on the tunnel surface:(1)Em′=1000001−D′/21+e75+25D′−GSI′/11,where D′ is the disturbance factor, in the form of Hoek–Brown failure criterion and GSI is the Geological Strength Index. The “prime” superscript denotes that the parameter belongs to the excavation damaged zone (EDZ). The disturbance factor was assumed to be zero (D′ = 0) for calculation of Young’s modulus in the plastic-EDZ boundary [51]. Assuming a linear interpolation between the values of Young’s modulus at the plastic-EDZ boundary and the tunnel surface, Young’s modulus within the different radiuses of the EDZ was calculated using the following equation [51]:(2)Ej′=E0′+Em′−E0′REDZρj′×REDZ−b,where Ej′, E0′Em′ represent Young’s modulus of jth annulus in the EDZ, on the plastic-EDZ boundary (D′ = 0), and on the tunnel’s surface (with a proper value for the disturbance factor), respectively. ρj′ is the radius of the jth annulus normalized by REDZ and b is the parameter of the intermediate principal stress. Figure 2 illustrates a summary of the procedure to find the stresses and strains in different zones/boundaries around the tunnel. Furthermore, Figure 3 represents the algorithm proposed by Ghorbani and Hasanzadehshooiili [26] which was used in the calculation of the support pressure (Pi) to generate the database of GRCs used in this study.Figure 3
The algorithm used in generation of the database of GRCs used in this study. After Ghorbani and Hasanzadehshooiili [26].


## 3. Developing a Relationship

### 3.1. Gathering Database for Model Development

The first step in finding a relationship for the GRC is to gather a database among the available sources to be used as training data for the neural network. Each GRC is constructed by a number of input-output datasets. The value of the internal support pressure (Pi) was taken as the target of prediction, which is function of different input parameters. Ghorbani and Hasanzadehshooiili [26] proposed a comprehensive algorithm for determination of GRC for circular tunnels in elastoplastic rock mass considering the EDZ and softening behavior. Besides, they utilized the unified strength criteria (USC) so that the effect of intermediate principal stress (b) and the exponential decaying dilation behavior can be taken into consideration. In order to verify and validate the proposed algorithm, they used a wide range of different sources including different GRCs under various conditions. In this study, the same data sources used in Ghorbani and Hasanzadehshooiili [26] were gathered as the database and then were used as the training data in the evolutionary polynomial regression (EPR) via the framework of genetic algorithm (GA). The summary of the data sources is as follows: a case of field measurements regarding the Hanlingjie tunnel in Hunan, China, at a depth of 146 m [17], two cases of Elastic-plastic rock mass considering strain softening [4, 9], a case of strain softening rock mass considering exponential decaying post-peak dilation parameter [52], and a case of elastic-perfect-plastic behavior of a rock mass considering EDZ [43]. Besides, Ghorbani and Hasanzadehshooiili [26] utilized the problem originally solved by Zareifard and Fahimifar [43] in order to cover the effects of new features of their proposed algorithm, such as the effects of the intermediate principle stress (b), the weight of the damaged rock mass and existence of a softening zone, and the exponential decaying dilation behavior. The data regarding these newly solved problems were also added to the database. All in all, a list of the input parameters included in the database along with their range of variation is presented in Tables 1 and 2, respectively.Table 1
List of affecting parameters on the GRC included in the database.
biRadius of the tunnelm′Material constant of the EDZΝPoisson ratio of the rock masss'Material constant of the EDZmiMaterial constant of intact rocka′Material constant of the EDZσ0In situ hydrostatic confining pressureγpStrain softening parameterGSIpPeak Geological Strength Index of the rock massEi′Young’s modulus in the inner annulus of the EDZDDisturbance factorE0′Young’s modulus in the outer annulus of the EDZEYoung’s modulus of rock massREDZRadius of the excavation-induced damaged zoneɸpPeak friction coefficientbIntermediate principal stress parameterɸrResidual friction coefficientUiConvergence (displacement)γp∗Critical deviatoric plastic strain parameterPiInternal support pressureTable 2
Range of variation of input and output parameters of the database.
Parameterbi(m)νmiσ0(MPA)GSIpDE(MPa)ɸpɸrγp∗Max.5.50.3517.773630500.655001251Min.50.25104.841.880.64000000STD.0.140.032.257.302.350.00434.793.511.630.25Mean5.050.2610.7227.6749.250.605361.111.510.820.07Parameterm′s'a′γpEi′(MPA)E0′(MPA)REDZ(m)bUi(m)Pi(MPA)Max.0.80.0010.5273038005500810.4391847Min.0.4283.95E − 050.50277.181257.186.500.0210STD.0.110.000.0114.161053.54986.460.370.350.062.60Mean0.770.000.5010.433313.983784.927.140.210.093.00


### 3.2. Sensitivity Analysis

The cosine amplitude method (CAM) can be used to find the strength of relationship between two parameters [25, 53]. Herein, this method is used to determine the degree of importance of each of the 19 input parameters in variation of the target (output) parameter, which is the internal support pressure (Pi).All the data pairs (each of the input parameters pairing with the target parameter) would form anX-array, such as X=x1,x2,x3,…,xi,xn where each element, xi, is a vector of the length of m:(3)xi=xi1,xi2,xi3,…,xim.By defining each parameter as dimensionless normalized vectors, the strength of relationship between the input and target parameters is given by(4)rij=∑k=1mxikxjk∑k=1mxik∑k=1mxjk,wherei,j=1,2,…,n,where xik and xjk are the normalized vectors of the input and target parameter and rij denotes the strength of relationship which would take a value between 0 and 1, where 0 implies no relationship and 1 shows the highest relationship. The name of this method, CAM, is taken from the fact that the result of the inner product of two vectors would be 0 if they are orthogonal vectors (meaning there is no relationship between them), while the result would be 1 if the two vectors are codirectional (meaning full relationship). It should be noted that the convergence (displacement) value was omitted from this analysis since it is as a fictional input parameter to the GRC.The results of CAM analysis are presented in Figure4. As can be seen, there are numerous parameters with relatively high contribution in the GRC. According to the results, the in situ hydrostatic pressure (s0) as well as the geometry of the tunnel (radius of the tunnel bi) presents a high strength of relationship in GRC. Interestingly, the properties of the excavation damaged zone make equal or even relatively higher contribution in the internal support pressure (Pi) than those of the intact rock mass. In fact, m′, s′, and a′ (the constant parameters of the Hoek–Brown failure criterion taken for the EDZ) with 77.9, 78.9, and 75.1% of contribution, the range of Young’s modulus assumed for the EDZ (Ei′, E0′) with 78.3 and 77.8%, and the extent of formation of the EDZ (REDZ) with 75.4 as well as its disturbance factor (D) with 75.5% of contribution are placed amongst the most important parameters in the GRC. This highlights the importance of considering the EDZ in design of tunnel as well as taking proper values for the material’s stiffness and strength of this zone.Figure 4
Results of sensitivity analysis using CAM.On the other hand, the properties of the intact rock mass (E, ν, mi, GSIp) have shown 77.1, 72.4, 69.1, and 76.5% of strength relationship, respectively. In contrast, the softening parameters (φp, φr, γp, and γp∗) made the least strength of relationship with 9.8–37.2%. Besides, it can be inferred that the role of material properties parameters of the rock mass and the EDZ is shown to be more meaningful in the value of support pressure, compared to that of the intermediate principal stress (b).


### 3.3. Evolutionary Polynomial Regression (EPR)

In this study, the genetic algorithm (GA) via the evolutionary polynomial regression (EPR) method is used to make a correlation between a number of 19 affecting parameters on the ground reaction curve (GRC), as the input parameters, and the internal support pressure of the tunnel (Pi), as the output parameter, i.e., the target of prediction. This method has successfully been used in developing mathematical models in different applications, especially in civil engineering [25, 54–61].Evolutionary polynomial regression (EPR), developed by Giustolisi and Savic [62], is known as a hybrid data-driven method. This method takes advantage of optimization in fitness to training data, yet by keeping the efficiency of mathematical expression, which is thanks to using a multiobjective search algorithm to establish multiple models. Hence, various optimal models are compared to pick the best fit such that the governing pattern can be created in a timely manner [63]. EPR applies different functions, components, forms, number of terms, and gens to present models in logarithmic, exponential, trigonometric, and inverse trigonometric sentences simultaneously by minimizing the error of prediction.In this study, the EPR was implemented in MATLAB environment. EPR provides a multiobjective genetic algorithm (MOGA) strategy in order to increase the regression’s parsimony.This strategy performs optimization of the model based on minimization of the number of inputs, number of terms, the error of prediction, or any combination of them. Herein, in order to maximize the model’s parsimony, the optimization was made with the target of least possible number of inputs and number of terms, along with the least sum of squared errors (SSE). Therefore, amongst the 19 input parameters, a number of 7 parameters (mi, GSIp, E, s′, a′, γp∗, Ei′) were automatically excluded from the model. This is probably because such parameters are interdependent variables where their effect can be incorporated by other parameters. The program automatically removes the less important or interdependent variables to increase the model’s parsimony via a procedure of trial and error to find the most efficient parameters to be used in the model. Also, as some parameters have similar effects on the GRC, it is believed that the program relies on the variation of Pi with some of them aiming at accurate predictions. Although some of the deleted parameters are still important to the problem, the remaining parameters with similar natures were enough to predict different values of Pi based on the concerning input parameters for the whole range of studied input parameters. As the proposed relationship is efficient for the development of GRCs of rock masses with input parameters, which lie in the range of parameters studied in this paper, some other deleted parameters may be still required for the representation of a new formula for other input ranges (useful for other stress, and material states/strength values). Hence, it is suggested to employ the proposed relationship for cases with adequate consistencies to the problem defined in this paper. Numerous functions in various forms were examined to find a relationship with the highest parsimony. The final relationship to predict the internal support pressure (Pi) of the tunnel is presented by the following equation:(5)PiMPa=1.3746×10−11ν−2σ03D−1m′REDZ3Ui−3−3.431×10−9φp−0.5+0.00056024×φr−0.25REDZ−0.25−3.9274×10−6φr−0.25b0.75Ui−1+1.704×10−13γp−0.5bUi−3−3.2565×10−10biν−1σ03E0′0.5REDZ−1Ui−3+0.16224.Figure5 illustrates a comparison of the predicted values of the support pressure verses the original values existing in the database. The quality of regression can be assessed by statistical measurements such as the coefficient of determination (COD = R2), the sum of squared errors (SSE), and the mean of squared errors (MSE). The COD was obtained as 90.33% which shows the good performance of the obtained mathematical model. Moreover, the SSE and MSE are obtained as 0.63% and 0.58%, respectively, verifying the decent precision of the proposed model. However, as Figure 5 suggests, the error of prediction mostly occurred in the range of support pressure less than 3 MPa. Therefore, it can be inferred that the model gives the best results in the range of 3–8 MPa, which is also the most common range of support pressure in practice.Figure 5
Comparison of the predicted support pressure with those of the database.


### 3.4. Verification and Comparison

In order to verify the validity of the obtained equation (equation (5)), a verification analysis was performed using two sets of available data in the literature. The first case was the results of field measurements at the crown of Hanlingjie tunnel in Hunan, China, reported by Zou et al. [17]. The tunnel is located at the depth of 146 m with the radius of 5.5 meters. It should be noted that some of the required data were back-calculated using proper relationships and correlations found in the literature, due to lack of data, especially regarding the EDZ. Moreover, the effect of intermediate principal stress was also investigated by varying parameter (b) in the range of 0–1. Table 3 presents a summary of the used parameters and Figure 6 depicts the results of calculations in the form of variation of radial displacements (Ui) against the internal support pressure. As can be seen, a meaningful consistency is observed as the required support pressure to restrict the displacement to a defined value can be acceptably predicted by the proposed relationship. It should be noted that some damaged parameters were not specified in this real case and the required values are calculated based on the available relationships. Therefore, the accuracy of predictions can be even higher while accessing to all the required in-field/laboratory measured parameters.Table 3
Values of the parameters used in the verification analysis.
ParameterZou et al. [17]Lee and Pietruszczak [9]bi (m)5.53ν0.350.25mi17.775.88s0 (MPa)4.815GSIp41.8855.25D0.60.6E (MPa)40005700ɸp1215ɸr55γp0.0080.01m′0.4280.6s′3.95E− 050.002a′0.5270.5γ′ (kN/m3)28.652.5Ei′ (MPa)277.185700E0′ (MPa)1257.185700REDZ7.54b0, 0.5, 10.5Figure 6
Verification analysis using field measurements data reported by Zou et al. [17].The second case is the study by Lee and Pietruszczak [9], in which the effect of strain softening was taken into account for the problem of circular tunnel in an elastic-plastic rock mass. The effects of EDZ and the intermediate principal stress were not considered in Lee and Pietruszczak [9]. Therefore, the proposed relationship is more practical for the cases of elastic-plastic-EDZ rock masses with strain softening in the plastic region. Hence, equivalent damaged parameters must be calculated prior to the analysis since the original studied case did not consider the occurred damaged zones. In this regard, damaged zone parameters are assumed to be identical to the residual zone parameters and the radius of the damaged zone is assumed to be 1 m (according to the values proposed in the literature). Other affecting and missing parameters are calculated using the relationships presented in Ghorbani and Hasanzadehshooiili [26]. The parameters used in this analysis are presented via Table 3 while Figure 7 depicts the results of comparison. As can be seen, a good consistency is observed between two studies which shows the strength and validity of the obtained equation by evolutionary polynomial regression (EPR) in this study while keeping simplicity and ease of use. This is important as the obtained equation can be used as a robust tool in the preliminary stages of design and construction in engineering practice by considering the state-of-the-art features of the problem such as strain softening, the intermediate principal stress, and the degradation caused by excavation damaged zone in the problem.Figure 7
Comparison of the results in this study with Lee and Pietruszczak [9].


## 3.1. Gathering Database for Model Development

The first step in finding a relationship for the GRC is to gather a database among the available sources to be used as training data for the neural network. Each GRC is constructed by a number of input-output datasets. The value of the internal support pressure (Pi) was taken as the target of prediction, which is function of different input parameters. Ghorbani and Hasanzadehshooiili [26] proposed a comprehensive algorithm for determination of GRC for circular tunnels in elastoplastic rock mass considering the EDZ and softening behavior. Besides, they utilized the unified strength criteria (USC) so that the effect of intermediate principal stress (b) and the exponential decaying dilation behavior can be taken into consideration. In order to verify and validate the proposed algorithm, they used a wide range of different sources including different GRCs under various conditions. In this study, the same data sources used in Ghorbani and Hasanzadehshooiili [26] were gathered as the database and then were used as the training data in the evolutionary polynomial regression (EPR) via the framework of genetic algorithm (GA). The summary of the data sources is as follows: a case of field measurements regarding the Hanlingjie tunnel in Hunan, China, at a depth of 146 m [17], two cases of Elastic-plastic rock mass considering strain softening [4, 9], a case of strain softening rock mass considering exponential decaying post-peak dilation parameter [52], and a case of elastic-perfect-plastic behavior of a rock mass considering EDZ [43]. Besides, Ghorbani and Hasanzadehshooiili [26] utilized the problem originally solved by Zareifard and Fahimifar [43] in order to cover the effects of new features of their proposed algorithm, such as the effects of the intermediate principle stress (b), the weight of the damaged rock mass and existence of a softening zone, and the exponential decaying dilation behavior. The data regarding these newly solved problems were also added to the database. All in all, a list of the input parameters included in the database along with their range of variation is presented in Tables 1 and 2, respectively.Table 1
List of affecting parameters on the GRC included in the database.
biRadius of the tunnelm′Material constant of the EDZΝPoisson ratio of the rock masss'Material constant of the EDZmiMaterial constant of intact rocka′Material constant of the EDZσ0In situ hydrostatic confining pressureγpStrain softening parameterGSIpPeak Geological Strength Index of the rock massEi′Young’s modulus in the inner annulus of the EDZDDisturbance factorE0′Young’s modulus in the outer annulus of the EDZEYoung’s modulus of rock massREDZRadius of the excavation-induced damaged zoneɸpPeak friction coefficientbIntermediate principal stress parameterɸrResidual friction coefficientUiConvergence (displacement)γp∗Critical deviatoric plastic strain parameterPiInternal support pressureTable 2
Range of variation of input and output parameters of the database.
Parameterbi(m)νmiσ0(MPA)GSIpDE(MPa)ɸpɸrγp∗Max.5.50.3517.773630500.655001251Min.50.25104.841.880.64000000STD.0.140.032.257.302.350.00434.793.511.630.25Mean5.050.2610.7227.6749.250.605361.111.510.820.07Parameterm′s'a′γpEi′(MPA)E0′(MPA)REDZ(m)bUi(m)Pi(MPA)Max.0.80.0010.5273038005500810.4391847Min.0.4283.95E − 050.50277.181257.186.500.0210STD.0.110.000.0114.161053.54986.460.370.350.062.60Mean0.770.000.5010.433313.983784.927.140.210.093.00


## 3.2. Sensitivity Analysis

The cosine amplitude method (CAM) can be used to find the strength of relationship between two parameters [25, 53]. Herein, this method is used to determine the degree of importance of each of the 19 input parameters in variation of the target (output) parameter, which is the internal support pressure (Pi).All the data pairs (each of the input parameters pairing with the target parameter) would form anX-array, such as X=x1,x2,x3,…,xi,xn where each element, xi, is a vector of the length of m:(3)xi=xi1,xi2,xi3,…,xim.By defining each parameter as dimensionless normalized vectors, the strength of relationship between the input and target parameters is given by(4)rij=∑k=1mxikxjk∑k=1mxik∑k=1mxjk,wherei,j=1,2,…,n,where xik and xjk are the normalized vectors of the input and target parameter and rij denotes the strength of relationship which would take a value between 0 and 1, where 0 implies no relationship and 1 shows the highest relationship. The name of this method, CAM, is taken from the fact that the result of the inner product of two vectors would be 0 if they are orthogonal vectors (meaning there is no relationship between them), while the result would be 1 if the two vectors are codirectional (meaning full relationship). It should be noted that the convergence (displacement) value was omitted from this analysis since it is as a fictional input parameter to the GRC.The results of CAM analysis are presented in Figure4. As can be seen, there are numerous parameters with relatively high contribution in the GRC. According to the results, the in situ hydrostatic pressure (s0) as well as the geometry of the tunnel (radius of the tunnel bi) presents a high strength of relationship in GRC. Interestingly, the properties of the excavation damaged zone make equal or even relatively higher contribution in the internal support pressure (Pi) than those of the intact rock mass. In fact, m′, s′, and a′ (the constant parameters of the Hoek–Brown failure criterion taken for the EDZ) with 77.9, 78.9, and 75.1% of contribution, the range of Young’s modulus assumed for the EDZ (Ei′, E0′) with 78.3 and 77.8%, and the extent of formation of the EDZ (REDZ) with 75.4 as well as its disturbance factor (D) with 75.5% of contribution are placed amongst the most important parameters in the GRC. This highlights the importance of considering the EDZ in design of tunnel as well as taking proper values for the material’s stiffness and strength of this zone.Figure 4
Results of sensitivity analysis using CAM.On the other hand, the properties of the intact rock mass (E, ν, mi, GSIp) have shown 77.1, 72.4, 69.1, and 76.5% of strength relationship, respectively. In contrast, the softening parameters (φp, φr, γp, and γp∗) made the least strength of relationship with 9.8–37.2%. Besides, it can be inferred that the role of material properties parameters of the rock mass and the EDZ is shown to be more meaningful in the value of support pressure, compared to that of the intermediate principal stress (b).


## 3.3. Evolutionary Polynomial Regression (EPR)

In this study, the genetic algorithm (GA) via the evolutionary polynomial regression (EPR) method is used to make a correlation between a number of 19 affecting parameters on the ground reaction curve (GRC), as the input parameters, and the internal support pressure of the tunnel (Pi), as the output parameter, i.e., the target of prediction. This method has successfully been used in developing mathematical models in different applications, especially in civil engineering [25, 54–61].Evolutionary polynomial regression (EPR), developed by Giustolisi and Savic [62], is known as a hybrid data-driven method. This method takes advantage of optimization in fitness to training data, yet by keeping the efficiency of mathematical expression, which is thanks to using a multiobjective search algorithm to establish multiple models. Hence, various optimal models are compared to pick the best fit such that the governing pattern can be created in a timely manner [63]. EPR applies different functions, components, forms, number of terms, and gens to present models in logarithmic, exponential, trigonometric, and inverse trigonometric sentences simultaneously by minimizing the error of prediction.In this study, the EPR was implemented in MATLAB environment. EPR provides a multiobjective genetic algorithm (MOGA) strategy in order to increase the regression’s parsimony.This strategy performs optimization of the model based on minimization of the number of inputs, number of terms, the error of prediction, or any combination of them. Herein, in order to maximize the model’s parsimony, the optimization was made with the target of least possible number of inputs and number of terms, along with the least sum of squared errors (SSE). Therefore, amongst the 19 input parameters, a number of 7 parameters (mi, GSIp, E, s′, a′, γp∗, Ei′) were automatically excluded from the model. This is probably because such parameters are interdependent variables where their effect can be incorporated by other parameters. The program automatically removes the less important or interdependent variables to increase the model’s parsimony via a procedure of trial and error to find the most efficient parameters to be used in the model. Also, as some parameters have similar effects on the GRC, it is believed that the program relies on the variation of Pi with some of them aiming at accurate predictions. Although some of the deleted parameters are still important to the problem, the remaining parameters with similar natures were enough to predict different values of Pi based on the concerning input parameters for the whole range of studied input parameters. As the proposed relationship is efficient for the development of GRCs of rock masses with input parameters, which lie in the range of parameters studied in this paper, some other deleted parameters may be still required for the representation of a new formula for other input ranges (useful for other stress, and material states/strength values). Hence, it is suggested to employ the proposed relationship for cases with adequate consistencies to the problem defined in this paper. Numerous functions in various forms were examined to find a relationship with the highest parsimony. The final relationship to predict the internal support pressure (Pi) of the tunnel is presented by the following equation:(5)PiMPa=1.3746×10−11ν−2σ03D−1m′REDZ3Ui−3−3.431×10−9φp−0.5+0.00056024×φr−0.25REDZ−0.25−3.9274×10−6φr−0.25b0.75Ui−1+1.704×10−13γp−0.5bUi−3−3.2565×10−10biν−1σ03E0′0.5REDZ−1Ui−3+0.16224.Figure5 illustrates a comparison of the predicted values of the support pressure verses the original values existing in the database. The quality of regression can be assessed by statistical measurements such as the coefficient of determination (COD = R2), the sum of squared errors (SSE), and the mean of squared errors (MSE). The COD was obtained as 90.33% which shows the good performance of the obtained mathematical model. Moreover, the SSE and MSE are obtained as 0.63% and 0.58%, respectively, verifying the decent precision of the proposed model. However, as Figure 5 suggests, the error of prediction mostly occurred in the range of support pressure less than 3 MPa. Therefore, it can be inferred that the model gives the best results in the range of 3–8 MPa, which is also the most common range of support pressure in practice.Figure 5
Comparison of the predicted support pressure with those of the database.


## 3.4. Verification and Comparison

In order to verify the validity of the obtained equation (equation (5)), a verification analysis was performed using two sets of available data in the literature. The first case was the results of field measurements at the crown of Hanlingjie tunnel in Hunan, China, reported by Zou et al. [17]. The tunnel is located at the depth of 146 m with the radius of 5.5 meters. It should be noted that some of the required data were back-calculated using proper relationships and correlations found in the literature, due to lack of data, especially regarding the EDZ. Moreover, the effect of intermediate principal stress was also investigated by varying parameter (b) in the range of 0–1. Table 3 presents a summary of the used parameters and Figure 6 depicts the results of calculations in the form of variation of radial displacements (Ui) against the internal support pressure. As can be seen, a meaningful consistency is observed as the required support pressure to restrict the displacement to a defined value can be acceptably predicted by the proposed relationship. It should be noted that some damaged parameters were not specified in this real case and the required values are calculated based on the available relationships. Therefore, the accuracy of predictions can be even higher while accessing to all the required in-field/laboratory measured parameters.Table 3
Values of the parameters used in the verification analysis.
ParameterZou et al. [17]Lee and Pietruszczak [9]bi (m)5.53ν0.350.25mi17.775.88s0 (MPa)4.815GSIp41.8855.25D0.60.6E (MPa)40005700ɸp1215ɸr55γp0.0080.01m′0.4280.6s′3.95E− 050.002a′0.5270.5γ′ (kN/m3)28.652.5Ei′ (MPa)277.185700E0′ (MPa)1257.185700REDZ7.54b0, 0.5, 10.5Figure 6
Verification analysis using field measurements data reported by Zou et al. [17].The second case is the study by Lee and Pietruszczak [9], in which the effect of strain softening was taken into account for the problem of circular tunnel in an elastic-plastic rock mass. The effects of EDZ and the intermediate principal stress were not considered in Lee and Pietruszczak [9]. Therefore, the proposed relationship is more practical for the cases of elastic-plastic-EDZ rock masses with strain softening in the plastic region. Hence, equivalent damaged parameters must be calculated prior to the analysis since the original studied case did not consider the occurred damaged zones. In this regard, damaged zone parameters are assumed to be identical to the residual zone parameters and the radius of the damaged zone is assumed to be 1 m (according to the values proposed in the literature). Other affecting and missing parameters are calculated using the relationships presented in Ghorbani and Hasanzadehshooiili [26]. The parameters used in this analysis are presented via Table 3 while Figure 7 depicts the results of comparison. As can be seen, a good consistency is observed between two studies which shows the strength and validity of the obtained equation by evolutionary polynomial regression (EPR) in this study while keeping simplicity and ease of use. This is important as the obtained equation can be used as a robust tool in the preliminary stages of design and construction in engineering practice by considering the state-of-the-art features of the problem such as strain softening, the intermediate principal stress, and the degradation caused by excavation damaged zone in the problem.Figure 7
Comparison of the results in this study with Lee and Pietruszczak [9].


## 4. Multiple-Variable Parametric Study

In this part, using the obtained relationship (equation (5)), a multivariable parametric study is conducted to investigate the simultaneous effects of various parameters such as the extent of excavation damage zone (REDZ) as well as its properties such as Young’s modulus at the EDZ-plasticity boundary and the constant parameter of m′ on the ground reaction curve. Besides, the effect of the intermediate principal stress parameter (b) is investigated. Such analysis is also helpful in verifying the obtained formula, by illustrating its sensitivity to variation of different parameters. In order to perform this analysis, all the variables were assumed to be constant equal to their default values, except the one being investigated.


### 4.1. Effects ofREDZ on GRC

Figure8 shows the effect of considering different values for the extent of the excavation damaged (EDZ) zone on the GRC by varying the REDZ between 6.5 and 8 meters. The same figure is represented in two different views for better observation. As can be seen, the internal support pressure (Pi) increases with increase of the REDZ. However, a critical boundary can be identified for displacements (obtained as 5–7 cm in this case), above which the REDZ has no effect on the internal support pressure.Figure 8
Effect ofREDZ on GRC: (a) first view, (b) second view.
(a)(b)Conversely, in the case of restricting the tunnel’s convergence to displacements lower than the critical boundary value (5–7 cm), the effect ofREDZ will be more considerable. In fact, this is implying that the effects of the extent of the damaged zone around the tunnel can be neglected provided that the installation time of the support system can have an enough delay, and the associated displacements are higher than the critical value. On the contrary, a special consideration should be paid to the extent of the damaged zone where only low displacements are allowed for the support system.


### 4.2. Simultaneous Effects ofE0′ and REDZ

Figure9 is presented to investigate the effects of both the extent of the damage zone (REDZ) and its strength properties, simultaneously. This figure shows the ground reaction curve with regard to the different radius of the damaged zone while three different values are taken for the elasticity modus at the EDZ-plastic boundary (E0′), namely, 1257, 3800, and 5500 MPa. As depicted, with increasing E0′, GRC drops. As can be seen, (4–6) cm is a critical boundary value for the tunnel’s convergence. For all cases, a higher internal pressure is needed to restrict displacements to values lower than 4 cm. The required internal support pressure increases for higher radiuses of the damaged zone. Such value will be even higher for a rock medium with lower E0′. On the other hand, the effect of E0′ on the GRC decreases as the tunnel’s displacement increases. Indeed, for displacements higher than (4–6) cm, effects of E0′ on the GRC will be negligible.Figure 9
Effect of variation of Young’s modulus of the EDZ on the GRC: (a) first view, (b) second view.
(a)(b)


### 4.3. Simultaneous Effects ofm′ and REDZ

Figure10 presents the simultaneous effect the extent of the excavation damage zone (REDZ) and m′ which represents the material’s constant of the excavation damaged zone. m′ was varied by taking three values: m′ = 0.428, 0.6, and 0.8 regarding the values available in the database. It can be said the GRC is sharply surged by increasing m′ even by small steps, leading to severe increase of the support pressure. This verifies the high contribution of the EDZ properties on the ground reaction curve of tunnel, as already suggested by the CAM analysis. Besides, the boundary value of the convergence is found to be 7 cm, after which the GRC would be constant.Figure 10
Effect of variation ofm′ on the GRC: (a) first view, (b) second view.
(a)(b)


### 4.4. Simultaneous Effects ofb (the Roll of Intermediate Principal Stress) and REDZ

Figure11 is presented to show the simultaneous effect of the intermediate principal stress (b) and REDZ on the GRC. For this purpose, the intermediate principal stress factor was varied between 0 and 1 taking five values in the steps of 0.25. As can be seen, the critical boundary value of convergence is (5–7) cm. Besides, it can be inferred that the value of b has a significant effect on the GRC, such that ignoring this parameter leads to overestimation of displacements and an overdesign supporting system, consequently. Similarly, increasing b-values result in decreasing displacements.Figure 11
Effect of intermediate principal stress andREDZ on the GRC: (a) first view and (b) second view.
(a)(b)


## 4.1. Effects ofREDZ on GRC

Figure8 shows the effect of considering different values for the extent of the excavation damaged (EDZ) zone on the GRC by varying the REDZ between 6.5 and 8 meters. The same figure is represented in two different views for better observation. As can be seen, the internal support pressure (Pi) increases with increase of the REDZ. However, a critical boundary can be identified for displacements (obtained as 5–7 cm in this case), above which the REDZ has no effect on the internal support pressure.Figure 8
Effect ofREDZ on GRC: (a) first view, (b) second view.
(a)(b)Conversely, in the case of restricting the tunnel’s convergence to displacements lower than the critical boundary value (5–7 cm), the effect ofREDZ will be more considerable. In fact, this is implying that the effects of the extent of the damaged zone around the tunnel can be neglected provided that the installation time of the support system can have an enough delay, and the associated displacements are higher than the critical value. On the contrary, a special consideration should be paid to the extent of the damaged zone where only low displacements are allowed for the support system.


## 4.2. Simultaneous Effects ofE0′ and REDZ

Figure9 is presented to investigate the effects of both the extent of the damage zone (REDZ) and its strength properties, simultaneously. This figure shows the ground reaction curve with regard to the different radius of the damaged zone while three different values are taken for the elasticity modus at the EDZ-plastic boundary (E0′), namely, 1257, 3800, and 5500 MPa. As depicted, with increasing E0′, GRC drops. As can be seen, (4–6) cm is a critical boundary value for the tunnel’s convergence. For all cases, a higher internal pressure is needed to restrict displacements to values lower than 4 cm. The required internal support pressure increases for higher radiuses of the damaged zone. Such value will be even higher for a rock medium with lower E0′. On the other hand, the effect of E0′ on the GRC decreases as the tunnel’s displacement increases. Indeed, for displacements higher than (4–6) cm, effects of E0′ on the GRC will be negligible.Figure 9
Effect of variation of Young’s modulus of the EDZ on the GRC: (a) first view, (b) second view.
(a)(b)


## 4.3. Simultaneous Effects ofm′ and REDZ

Figure10 presents the simultaneous effect the extent of the excavation damage zone (REDZ) and m′ which represents the material’s constant of the excavation damaged zone. m′ was varied by taking three values: m′ = 0.428, 0.6, and 0.8 regarding the values available in the database. It can be said the GRC is sharply surged by increasing m′ even by small steps, leading to severe increase of the support pressure. This verifies the high contribution of the EDZ properties on the ground reaction curve of tunnel, as already suggested by the CAM analysis. Besides, the boundary value of the convergence is found to be 7 cm, after which the GRC would be constant.Figure 10
Effect of variation ofm′ on the GRC: (a) first view, (b) second view.
(a)(b)


## 4.4. Simultaneous Effects ofb (the Roll of Intermediate Principal Stress) and REDZ

Figure11 is presented to show the simultaneous effect of the intermediate principal stress (b) and REDZ on the GRC. For this purpose, the intermediate principal stress factor was varied between 0 and 1 taking five values in the steps of 0.25. As can be seen, the critical boundary value of convergence is (5–7) cm. Besides, it can be inferred that the value of b has a significant effect on the GRC, such that ignoring this parameter leads to overestimation of displacements and an overdesign supporting system, consequently. Similarly, increasing b-values result in decreasing displacements.Figure 11
Effect of intermediate principal stress andREDZ on the GRC: (a) first view and (b) second view.
(a)(b)


## 5. Conclusion

In this study, an attempt was made to present a simple, yet accurate and comprehensive method using the genetics algorithms (GA), to construct the ground reaction curve (GRC) of circular tunnels in elastoplastic rock mass considering the state-of-the-art feature of the problem including strain softening, effect of excavation damaged zone (EDZ), the effect of principal stress parameter, and the exponential decaying dilation behavior via the unified strength criterion (USC) framework as the material strength criteria. Incorporating such advanced issues into the problem usually demands advanced numerical analysis in a time-consuming procedure and a high level of knowledge of rock mechanics, constitutive models, etc. The results of this study can be beneficial in the preliminary stages of tunnel design in a simple, yet accurate and comprehensive manner.A database was gathered from the available data sources available in the literature including various tunnels in different conditions with different assumptions. Besides, the database was enriched by adding some new advanced features such as the effect of the intermediate principal stress, weight of EDZ, and the exponential decaying dilation parameter obtained by a newly developed algorithm. A relationship was produced using evolutionary polynomial regression (EPR) technique, which uses the best features of the genetic algorithm (GA) and the least squared method to predict new target functions based on their previous learning of the relationship between input and output parameters. This relationship can be used, via a design spreadsheet, for preliminary design of the supporting system of the tunnel in a timely manner and high accuracy. Validity of the results was confirmed using verification against field measurements data and comparison to other advanced numerical solutions. Sensitivity analysis of the database using the cosine amplitude method showed firstly the importance of taking EDZ into account in the design procedure and secondly choosing proper values for stiffness and strength of this zone. Besides, a multivariable parametric study was performed to show the simultaneous effects of the extent of excavation damage zone (REDZ) and taking proper properties such as m′ (material’s constant for the EDZ) and E′ (Young’s modulus of the EDZ) for this zone on the GRC. In addition, the effect of intermediate principal stress on the GRC was significant, such that ignoring this parameter results in an overdesign supporting system. The boundary values of the tunnel convergence were obtained 5–7 cm for the REDZ, 4–6 cm for Young’s modulus of the EDZ-plasticity boundary (E0′), 7 cm for the material constant of EDZ (m′), and 5 cm for the intermediate principal stress (b). Most importantly, the boundary value of the tunnel convergence showed a relatively strong dependency on the extent of EDZ. Therefore, in case of restricting convergence conditions, i.e., where a small value of tunnel convergence (lesser than the aforementioned boundary values) is desired for the tunnel, the radius of the excavation damaged zone (REDZ) should be determined with cautious.Last but not least, it should be remembered that the method originated from the ground reaction curve in the framework of the confinement convergence method. Therefore, the same assumptions and limitations are applied, as mentioned above (homogenous rock mass, isotropic stress field, circular cross section of the tunnel, plain-strain condition, and neglecting the gravity forces). Besides, the rock quality should be in the same range with the values considered in this study (medium quality rock mass) and the small strain condition is applied.



---

*Source: 2237918-2021-08-29.xml*