- Email: [email protected]
A modified optimization algorithm for back analysis of properties for coupled stress-seepage field problems
Tunnelling and Underground Space Technology xxx (xxxx) xxxx
Contents lists available at ScienceDirect
Tunnelling and Underground Space Technology journal homepage: www.elsevier.com/locate/tust
A modified optimization algorithm for back analysis of properties for coupled stress-seepage field problems ⁎
Chuangzhou Wua,b, Yi Honga, Qingsheng Chenc, , Shivakumar Karekald a
College of Civil Engineering and Architecture, Zhejiang University, Hangzhou 310058, China Department of Geophysics, Kangwon National University, Chuncheon, Gangwon-do 200701, Republic of Korea c Department of Civil and Environmental Engineering, National University of Singapore, 21 Lower Kent Ridge Road, 119077, Singapore d School of Civil, Mining & Environmental Engineering, University of Wollongong, Wollongong City, NSW 2522, Australia b
A R T I C LE I N FO
A B S T R A C T
Keywords: Back analysis Coupled stress and seepage fields Levenberg-Marquardt optimization Complex-variable-differentiation method Parameters estimation
The study is aimed at developing a new algorithm for back analysis of rock mass parameters based on the observed displacements of tunnel excavations under coupled stress-seepage field problems which are generally encountered in hydraulic tunnel projects. The back-analysis algorithm is developed by incorporating the Levenberg-Marquardt optimization technique with complex-variable differentiation method. Additional auxiliary technique is also incorporated to enhance the convergence and stability of the proposed algorithm during the back-estimation of multiple rock mass parameters. A hypothetical hydraulic tunnel case was used for testing and validating the proposed algorithm by incorporating the method in a finite element code in which multiple rock mass parameters (such as modulus, permeability, and in situ stress) were treated as target parameters. Results show that the multiple rock mass parameters can be accurately and efficiently estimated by back analysis using a newly developed algorithm for coupled stress and seepage fields encountered in tunneling. The proposed algorithm can be used for predicting excavation behavior, particularly, the stress-induced deformations at subsequent stages of tunnel excavation under coupled multiple fields (e.g. stress and seepage fields).
1. Introduction Stress-seepage coupling analysis is commonly undertaken in geotechnical engineering especially in the design of hydraulic tunnels. Accurate estimation of rock mass properties is necessary for design of geotechnical structures. However, reliable estimation of rock mass properties is difficult in the laboratory. It is not only expensive but also several errors associated with disturbances in sample collection, testing and scale effects of the rock mass which make it difficult to estimate true properties of rock mass. In couple problems involving stress and seepage, it is all the more difficult to assess suitable properties of rock mass (Exadaktylos, 2001; Finno and Calvello, 2005; Gens et al., 1996; Gonzaga et al., 2008; Hakala et al., 2007 and Du et al., 2018; Li et al., 2019; Maleki, 2018; Vardakos et al., 2016; Yang et al., 2017). Furthermore, it is also time-consuming and economically undesirable to determine the mechanical and hydraulic properties of anisotropic rock mass especially in deep and interaction of the multi-field (e.g. stress and seepage coupled field) (Chen et al., 2019; Gonzaga et al., 2008; Hakala et al., 2007). Therefore, the back analysis method based on field observations, as suggested by Peck (1969), becomes important in
⁎
determining the representative and/or dominant geomaterial parameters in practice (Zhang et al., 2015; Goh et al., 2018). The backestimated parameters can then be used in the subsequent predictions of behavior or rock mass under different geometry and boundary conditions. Therefore, it is of great significance to develop an efficient and accurate approach for the parameter determination of fractured rock mass (Exadaktylos, 2001; Gao et al., 2018; Li et al., 2019; Luo et al., 2018; Sun et al., 2018; Tang and Kung, 2009; Zhang and Goh, 2015). The application of the back analysis technique often consists of two procedures: a forward procedure and a backward procedure (Qi and Fourie 2018; Qi et al., 2018). Optimization methods thus utilized during the backward procedure which can be classified into two categories: that base on the artificial intelligence (AI) theory (Qi and Fourie 2018; Qi et al. 2018) the classic optimization method. The optimization methods have been applied to calibrate the geomaterial parameters through back analyses of field observations over the past two decades (Sakurai and Takeuchi, 1983, Arai et al., 1983; Tang and Kung, 2009; Zhang et al., 2015; Wu et al., 2019; Yao et al., 2019). As the optimization method suffers from cancellation errors and slow convergence rates in addition to the likelihood of becoming trapped in local optima.
Corresponding author. E-mail addresses: [email protected] (C. Wu), [email protected] (Y. Hong), [email protected] (Q. Chen), [email protected] (S. Karekal).
https://doi.org/10.1016/j.tust.2019.103040 Received 2 January 2019; Received in revised form 14 March 2019; Accepted 12 July 2019 0886-7798/ © 2019 Elsevier Ltd. All rights reserved.
Please cite this article as: Chuangzhou Wu, et al., Tunnelling and Underground Space Technology, https://doi.org/10.1016/j.tust.2019.103040
Tunnelling and Underground Space Technology xxx (xxxx) xxxx
C. Wu, et al.
iteratively in the optimization process in search of global optimized values. The iterative procedure for modified L-M method used in this study can be expressed as
Thus, the complex-variable differentiation method was developed to evaluate derivatives of real functions in optimization methods (Lyness and Moler, 1967; Gao, 2002). This method is efficient since it has neither cancellation errors nor step-size-dependent problems compared to normal optimization methods. However, in dealing with complex coupled problems (i.e Stressseepage issues or Stress-seepage-thermal issues) involving multiple parameters, the optimization method often fail to converge and results in inaccurate values (Cai et al., 2007). In this study, a modified algorithm is proposed by incorporating the auxiliary techniques in the Levenberg-Marquardt (here mentioned as “L-M”) optimization technique to enhance the convergence and stability of the algorithm in dealing with coupled problems involving multiple fields such as stress-seepage coupling issues and Stress-seepage-thermal coupling issues. The proposed new method is then incorporated in a finite element code. A series of hypothetical excavation cases were used to test and validate the applicability of the proposed method for various scenarios of multiple parameter groups of coupled stress-seepage field problems.
x (k + 1) = x k + α (k + 1) Sk
(2)
Sk = (DkT Dk + μI )−1DkT (um (x k ) − ui*)
(3)
∂u ⎡ ∂x11 ⋯ ⎢ Dk = ⎢ ⋮ ⋱ ⎢ ∂um ⋯ ⎣ ∂x1
(4)
The modified optimization method was developed for the back analysis of excavation to back-figure multiple-parameters simultaneously based on observed deformations (e.g., displacement of the tunnel). The modified algorithm is proposed by incorporating the auxiliary techniques in the L-M optimization technique to enhance the convergence and stability of the algorithm. In mathematics and computing, the Levenberg-Marquardt algorithm (LMA), also known as the damped least-squares (DLS) method, is used to solve non-linear least squares problems. The LMA is more robust than the Gauss–Newton algorithm (GNA), which means that in many cases it finds a solution even if it starts very far off the final minimum. The improved optimization method is capable of conducting multi-parameters back analyses simultaneously of coupled stress and seepage field, which is described in the following sections. Note that such back-figured rock parameters should be regarded as the equivalent rock parameters, not real rock parameters (Tang and Kung, 2009; Zhang, et al., 2015; Espada and Lamas, 2017). Specifically, the deformation behavior of tunnel may be affected by a large number of factors such as the excavation width and depth, surround rock stiffness strength distribution, bolt arrangement, lining spacing and stiffness. The back-figured parameters are generally away from the real parameters and thus regarded as the equivalent parameters.
2.3. Convergence criteria In the proposed study, the objective function need to satisfy the following convergence criterion:
J (x k + ak Sk ) < J (x k )
(5)
The convergence criterion should satisfy the following conditions in order to terminate the process for obtaining reasonable estimates.
J (x k ) ≤ ε
(6)
whereε is tolerance which is a small positive number. 2.4. Auxiliary techniques dealing with complex issues To enhance convergence rate of the proposed optimization method, especially for the coupled cases such as stress and seepage field, the auxiliary technique is employed. In this study, a new relaxation factor was employed by substituting the step size factor, α k (marked as ‘SSF’ method) by step size matrix Y k (marked as SSM method), and Y k has the following expression
Yk
2.1. Objective function An objective function is necessary in a back analysis to determine the optimized rock mass parameters through iteration process. The objective function J (x *) used in this study is defined as
k ⎡ a1 ⋯ 0 ⎤ ⎢ = ⋮ aik ⋮ ⎥ ⎥ ⎢ ⎢ 0 ⋯ amk ⎥ ⎦n ⎣
(7)
aik
is the step size factor for each type of parameters in coupled where (stress-seepage) case scenario. By introducing the step size matrix Y k into the iterative procedure (Eq. (2)), the optimum step size for each type of parameters can be determined and assigned accordingly, which would enhance the convergence rate of the optimization process. In contrast, for normal optimization method only one step size for all types of parameters are used (Tang and Kung, 2009).
m i=1
⎤ ⎥ ⋮ ⎥ ∂um ⎥ ∂xn ⎦ k
where k is the iteration number and S is a search direction. The step size factor in terms of the scalar quantity α k defines the moving distance in the search direction, S . Dk , denotes the sensitivity matrix. n is the number of the parameters. In this test, the complex-variable-differentiation method was employed to calculate the sensitivity matrix Dk rather than direct differentiation used in normal optimization algorithms which will give rise to the cancellation errors or step-size-dependent problems. The algorithm of complex-variable-differentiation method (CVDM) is presented in detail in following sub-sections.
2. Proposed hybrid optimization algorithm
J (x *) = minJ (x ) = min ∑ (ui (x ) − ui*)2
∂u1 ∂xn
(1)
where, m is the number of measurement points used for back-predicting the intended multiple parameters in an optimization analysis. ui* is the measured displacement and ui is the computed displacement using FEM. x is a parameter vector composed of the intended multiple parameters. Note that the multiple types of parameters such as modulus, poison’s ratio, permeability, in situ stress and water pressure can be simultaneously calculated based on observed displacement of the tunnel.
2.5. Complex-variable-differentiation method In order to improve the convergence rate of the optimization process and to improve the accuracy of back-estimated multiple parameters for couple multi-field, the complex-variable-differentiation method (CVDM) developed by Lyness and Moler (1967) was incorporated in L-M method. This would eliminate the cancellation error and heavily dependent step-size effect. In this CVDM, the variable x of a real function u (x ) is replaced by a complex one of x + ih , with the imaginary part h which is very small (usuallyh = 10−20 y). The function u (x + ih) can be expanded into
2.2. Search direction and line search An initial set of target variables, x0 (i.e., initial condition) is required in optimization algorithms, and then these variables are updated 2
Tunnelling and Underground Space Technology xxx (xxxx) xxxx
C. Wu, et al.
Fig. 1. 2D FEM model (a) Tunnel geometry, (b) Mesh around the tunnel and (c) Layout of convergence measurement.
matrix, see Eqs. (4)) can be obtained by dividing the Im [u (x + ih)] with h as
Table 1 Mechanical and hydraulic property of the rock materials. Rock mass
Self-weight (kN/m3)
Eh (GPa)
Ev (GPa)
kh (10−8 m/s)
kv (10−8 m/s)
Type 1 Type 2 Type 3
27.2 24.0 26.1
15 1.1 16
8.0 1.8 22
5.14 1.92 0.63
4.2 1.92 0.63
du Im [u (x + ih)] = dx h
where Im denotes the imaginary part. It is shown that CVDM is a very promising method since the derivatives require to be evaluated. It can be used as easy as the direct difference method and has neither cancellation errors nor step-size-dependent problems. In CVDM, the variables should be complex types, and the derivatives of an inverted parameter can be obtained by adding a small imaginary part of ih . The numerical method, such as the direct differentiation method, is the easiest way to compute the sensitivity coefficients of Eqs. (4).
Note: E – Elastic modulus of rock mass, k – Hydraulic conductivity, Subscripts of h and v – indicate the relative parameters in the horizontal and vertical directions representatively
Taylor’s series as
u (x + ih) = u (x ) + ih
du h2 d2u − + o (h3) dx 2 dx 2
(9)
2.6. Solution procedure
(8)
The computational steps based on the displacement can be
Since h is very small, derivatives of a real function (i.e. sensitivity Table 2 Case studies information. Target parameters
Case 1 Case 2 Case 3
Auxiliary technique
Modulus
Hydraulic Conductivity
Field stress
Eh1, Ev1, Eh2, Ev2, Eh3, Ev3 Eh1, Ev1 Eh1, Ev1, Eh2, Ev2
– k1h, k1v k2
ph, pv ph, pv
Note: Initial values setting: Moduli (E) as 1 GPa, Hydraulic conductivity (k) as 10 × 10−8 m/s and field stress (p) as 10 MPa. 3
SSF SSF SSM (compared with SSF)
Tunnelling and Underground Space Technology xxx (xxxx) xxxx
C. Wu, et al.
25
20
J16=2.6*10-13
Eh3
Eh1 Ev1
10
J5=4*10-4
Ev2 Eh2
J0=0.246
0
Initial values Values after 5 steps of iteration Values after 13 steps of iteration True values
20
J13=4.7*10-15
E1h
15
10
E1v J5=3.2*10-7
5
J0=1.6*10-3
0 0
-10 -10
0
10
20
30
5
10
Calculated permeability, k (10-8 m/s)
Initial values Values after 5 steps of iteration Values after 16 steps of iteration
pv ph
J16=2.6*10-13
20
Iteration process
J5=4*10-4 J0=0.246
10
J0=1.6*10-3
10
Iteration process
8
kh1
6
J5=3.2*10-7
kv1
J13=4.7*10-15
4
Initial values Values after 5 steps of iteration Values after 13 steps of iteration True values
2 0 0
2
4
6
8
10
12
Measured permeability, k (10-8 m/s)
(b)
0 0
10
20
30
40
50
50
Calculated loading pressure, p (MPa)
Measured loading pressure, p (MPa)
(b) Fig. 2. The iteration process of inverted parameters for case 1 (Eh1, Ev1, Eh2, Ev2, Eh3, Ev3, ph, and pv) with objective function value after different iterations.
summarized as follows: Step 1. Input measured displacements, u* and initialize target parameters xi. Step 2. Call the subroutines (except the Step 3 in the direct problem) to obtain back-estimated values of parameters xi. Step 3. Terminate the iterative procedure if the stopping criterion is achieved; otherwise, employing SSM method (Eqs. (7)) and then go to Step 2.
Initial values Values after 5 steps of iteration
40
Values after 13 steps of iteration True values
J13=4.7*10-15 J5=3.2*10-7
pv
30
ph 20
10
J0=1.6*10-3
Iteration process
Calculated loading pressure, p (MPa)
25
(a)
50
30
20
12
(a)
True values
15
Measured modulus, Em (GPa)
40
Measured Modulus, Em (GPa)
40
Iteration process
30
Calculated modulus, Ec (GPa)
Initial values Values after 5 steps of iteration Values after 16 steps of iteration Ev3 True values
Iteration process
Calculated Modulus, Ec (GPa)
40
0 0
10
20
30
40
50
Measured loading pressure, p (MPa) 3. FEM code and Stress-seepage coupling analysis
(c)
The proposed optimization method was compiled and incorporated into the FEM code on MATLAB platform. A hypothetical case is considered for testing the proposed algorithm. A tunnel is excavated in deep rock mass with three different types of stratum bounded by high in situ stress and seepage pressure (He et al., 2011; Yang et al., 2014; Ma et al., 2017; Grenon et al., 2017). In terms of the seepage-stress coupling analysis, the following assumptions should be confirmed. (i) initially, the rock mass is treated as an elastic, homogeneous and anisotropic; and (ii) seepage tensor is defined in the Representative Element
Fig. 3. The iteration process of inverted parameters for case 2: (a) Eh1 and Ev1 (b) ph, and pv (c) kh1 and kv1 with objective function value at different iteration number.
Volume (REV) of the rock mass (Yang et al., 2014). For orthotropic rock mass in plane problems, the elastic constitutive material model with no pore water pressure is given as by:
4
Tunnelling and Underground Space Technology xxx (xxxx) xxxx
C. Wu, et al.
10 E2v
E1v
E1h
8
6
10
Using SSM code after 10 iteration of SSF
Real value
Modulus, E (GPa)
k2
4
5 2
0
Initial value
0
5
10
15
20
25
( V = − (K y
u12 u23 u13 u14 u45 u15 u16 u67 u17
19.847 21.428 30.531 38.805 47.537 52.181 37.676 39.485 31.388
Relative errors(%)
Identified convergence
Relative errors(%)
19.885 21.473 30.622 38.917 47.716 52.358 37.723 39.586 31.472
0.19 0.21 0.30 0.29 0.38 0.34 0.12 0.26 0.27
19.856 21.445 30.559 38.834 47.575 52.227 37.697 39.506 31.406
0.09 0.08 0.09 0.07 0.08 0.09 0.06 0.05 0.06
2
∂H ∂y
(12)
(13)
4.1. Model settings
Identified convergence
2
+ Kyy
4. Validation of the proposed optimization algorithm
0
A two-dimensional plane strain finite element model was considered for back predicting the properties of rock mass surrounding a tunnel excavation subjected to stress and seepage fields (coupled fields). A horse-shoe shaped tunnel with a height of 8 m and a width of 6 m is considered. The boundary is set to about 4 times of width of the tunnel. The mechanical and hydraulic properties of the rock mass were collected from a practical project completed in China (Wu et al., 2018). The rock mass was subjected to in-situ stress with major in-plane principal stress (pv) of 30 MPa in vertical direction and the minor inplane principal stress ph of 24 MPa in horizontal direction. The seepage pressure of the rock mass at the boundary is considered as 10 MPa and no water head inside of the tunnel. The seepage direction is from outside to inside of the tunnel. The constitutive model of the rock mass in this study is the linear elastic model. The model setting and meshes are shown in Fig. 1a and b, respectively. Several points surrounding the tunnel surface are selected as measuring points as shown in Fig. 1c, thus the convergence of the tunnel face could be represented by the relative displacements between the measuring points uij. Here, nine convergence measurements are adopted for the following studies. The rock mass mechanical parameters are listed in Table 1, which are further adopted as target values in the back analysis. The tunnel is driven through stratified rock mass consists of three different rock types, characterizing anisotropic mechanical behavior. In the following section, these mechanical parameters and the field stress are compared with those calculated from the proposed back analysis approach with nine convergence measurements. To illustrate the application of the proposed method, three case studies are carried out as shown in Table 2. In general, Case 1 and Case 2 are adopted to assess the efficiency and accuracy of the proposed optimization method through the determination of different types of parameters, such as modulus, permeability, in situ stress, etc. individually (separately) and combinedly (coupled manner). On the other hand, it should be noted that the 2nd type of rock mass has much lower modulus and larger hydraulic conductivity than the other two types which would be essential for the tunnel stability analysis. Therefore, Case 3 is used to present a more crucial case which is dealing with nonconvergence situations due to the high variation of sensitivity of different type of parameters which need be back-estimated
SSM
v v ⎡ 1 − v31 − ⎛ E12 + E31 ⎞ E3 3⎠ ⎢ E1 ⎝ 1 ε11 ⎡ ⎤ ⎢ 2 2 v v v 1 − E31 ⎢ ε22 ⎥ = ⎢− ⎛ E12 + E31 ⎞ E1 3⎠ 3 ⎣ ε12 ⎦ ⎢ ⎝ 1 ⎢ 0 0 ⎢ ⎣
∂H xy ∂x
)⎫⎪ )⎬⎪⎭
where Kf is the current groundwater hydraulic conductivity; K 0 is the initial hydraulic conductivity; σ is the normal stress on the joint plane, and α is the coupling parameter that reflect the influence of stress.
Table 3 Comparisons of the tunnel convergence in case 3. SSF
∂H ∂y
Kf = K 0 e−ασ
Fig. 4. Improved convergence of inverting parameters for case 3 by using SSM after termination of iteration using SSF code at iteration n = 10.
Measured values
+ Kxy
The rock mass permeability is controlled by the normal stress. The stress is coupled to seepage according to Eq. (11) by Louis (1972).
Iteration number, n
Convergence
∂H ∂x
Vx = − Kxx Permeability, k (10-8 m/s)
15
E2h
⎤ ⎥ ⎥ ⎡ σ11 ⎤ 0 ⎥ ⎢ σ22 ⎥ ⎥ ⎣ σ12 ⎦ 2(1 + v12 ) ⎥ ⎥ E1 ⎦ 0
(10)
where E1 and E3 are the Young’s moduli and v31 and v12 are Poisson’s ratio The relationship between the hydraulic pressure and the hydraulic conductivity, based on the Biot’s theory (Biot, 1941), is expressed as
kij ∇2 P = 0
(11)
where kij is the hydraulic conductivity, P is the hydraulic pressure. The equation of hydraulic conductivity for a plane strain model can be expressed as:
Table 4 Analysis results of the hypothetical case 3 using SSF and SSM code. Parameter type
Eh1 (GPa) Ev1 (GPa) Eh2 (GPa) Ev2 (GPa) k2 (10−8 m/s)
True values
15 8 1.1 1.8 1.92
Initial values
1 1 1 1 1
SSF
SSM
Convergence values
Relative error (%)
Convergence values
Relative error (%)
14.6 7.7 1.12 1.86 2.41
2.7 3.9 1.8 5.4 25.7
15 8.0 1.1 1.8 1.93
< 0.01 < 0.01 < 0.01 < 0.01 0.42
5
Tunnelling and Underground Space Technology xxx (xxxx) xxxx
C. Wu, et al.
such a case, the residuals of the computed objective function for this case are as low as 8 × 10−12, meaning that the solution is converged for the back-estimated two kinds of parameters. The results of case 3 suggest that the proposed SSM method improved the optimization process and it seems to be more suitable for the back analysis when multiple parameters are encountered. The case 3 is a theoretical case study in which the exact rock mass properties are known, which allowed to realize that the expected solution in case 3 did not converge with the SSF method, and so it was necessary to apply the SSM. As far as field application was concerned, the displacements obtained in case 3 with the SSF method shows that the SSF method is quite satisfactory for back-analysis of rock mass properties including coupling issues the relative errors of which is lower than 1%. However the proposed method using SSM method shows higher precision for complex multiple parameters for multiple field coupling issues (i.e. seepage field, stress field and thermal field).” Based on the results of the case study, the proposed optimization model appears to be sound and suitable for back analysis for the coupled seepage-stress problems. However, it should be noted that, further cross validation will be needed on the accuracy of the estimation of rock mass parameters for the field case with measured in-situ data in order to use the proposed optimization algorithm for engineering applications.
simultaneously. 4.2. Analysis of case 1 and case 2 In case 1, the back-calculated parameters were set as the moduli (Eh1, Ev1, Eh2, Ev2, Eh3, Ev3) and the field stresses as (ph, pv). The initial values of moduli are set as 1.0 GPa, while those of field stresses are set as 10.0 MPa. Fig. 2 presents the convergence process of moduli and filed stresses after the 0rd, 5th and 16th iteration with their objective functions. In general, the objective function reduces with the increase in the number of iterations. It can be seen that the back-calculated parameters approaches the target values after the 5th iteration as the residual of the objective function becomes J5 = 4 × 10−4. The analysis stopped at the 16th iteration with J16 = 2.6 × 10−13. In case 2, the back-calculated parameters were set as the moduli, the hydraulic conductivity and the field stresses. Here, the parameters of rock type 1 are studied in this case, such as the moduli (Eh1, Ev1), the hydraulic conductivity (k1h, k1v) and the field stresses (ph, pv). As shown in Fig. 3a and Fig. 3b, the back-calculated hydraulic conductivity and field stresses are close to the targets values after 5 iterations with J5 = 3.2 × 10−7. However, the back-calculated moduli needed several iterations more to achieve the target value. Finally, all the back-calculated parameters could be close to the target values after 13th iteration (J13 = 4.7 × 10−15). The results from case 1 and case 2 indicate that the proposed backanalysis scheme is useful in identifying the multiple types of parameters such as modulus, field stress and hydraulic conductivity in dealing with coupling issues. More iterations for case 1 are needed in search of global solution than that of case 2, which is perhaps due to initial conditions of number of moduli parameters which would have allowed the solution to get caught into local minimum, thus increasing the iteration number until it reaches the final solution or final minimum. In case 2, the proposed optimization model seems effective in dealing with three types of parameters simultaneously by considering coupling effect of the hydraulic parameters and in-situ stresses.”
5. Concluding remarks A modified optimization technique is proposed, tested and validated to back predict the multiple properties simultaneously for the coupled stress-seepage related problems based on the observed displacements. The proposed optimization algorithm can also be used for coupled problems in anisotropic rock mass material models. Some of the main derived conclusions are: (i) The proposed optimization algorithm shows better precision in estimating parameters such as modulus, permeability and in situ stress, simultaneously, in coupled problems. In case 2, for different types of parameters, the proposed optimization model seems more sensitive to the hydraulic parameters and in-situ stresses than modulus. (ii) In dealing with the complex problem such as the target parameters which are associated with coupling issues of anisotropic materials (see case 3), an auxiliary technique by substituting the step size factor (SSF method) with factor matrix (SSM method) can be employed to improve the convergence, which is useful in solving the large convergence problems with low relative errors.
4.3. Analysis of case 3 In case 3, the back-calculated parameters were set as the moduli of both rock types 1 and 2 and the hydraulic conductivity of rock mass type 2. The SSF method is firstly adopted to back-analyze the target parameters. Totally 10 iterations are carried out using SSF method. As shown in Fig. 4, the moduli (Eh1, Ev1, Eh2, Ev2) are approaching to the target values as number of iterations increases. However, the hydraulic conductivity k2 could not be well identified at the 10th iteration and the relative error between the identified value and measured value is as large as 25.7%. Looking at the tunnel convergences shown in Table 3, they are also not well identified. This might be caused by the high deformability and high hydraulic conductivity of rock type 2. In the model, the convergence measurement u15 becomes as high as 52.2 mm which is much larger compared to other locations. The other two measurement u14 and u45 are also relatively large which could affect the accuracy of the calculated k2 value. To improve the identification of hydraulic conductivity k2, the auxiliary technique (SSM method) is further adopted by introducing a matrix factor to substitute the step size factor. As shown in Fig. 4, k2 achieved the target value progressively using SSM method. The error between the identified value and the measured one is as low as 0.42% after another ten iterations. At the same time, there is no influences on the identified modulus (Eh1, Ev1, Eh2, Ev2) which are approaching the target values as well. The increasing accuracy to identify the rock properties are listed in the Table 4 that the errors of moduli become less than 0.1% after SSM method and hydraulic conductivity becomes as low as 0.42%. As shown in Table 3, the accuracies of the tunnel convergence are also improved after using SSM method and the errors between the identified and measured values become less than 0.1%. In
Acknowledgement This work was financially supported by State Key Laboratory for GeoMechanics and Deep Underground Engineering, China University of Mining & Technology, China (Grant Nos: SKLGDUEK1521 and SKLGDUEK1819), National Natural Science Foundation of China (Grant No: 51509219), and Basic Science Research Program through the National Research Foundation of Korea (NRF) funded by the Ministry of Education (2019R1A6A1A03033167). References Arai, K., Ohta, H., Yasui, T., 1983. Simple optimization techniques for evaluating deformation moduli from field observations. Soils Found. 23 (1), 107–113. Biot, M.A., 1941. General theory of three-dimensional consolidation. J. Appl. Phys. 12 (2), 155–164. Cai, M., Morioka, H., Kaiser, P.K., Tasaka, Y., Kurose, H., Minami, M., Maejima, T., 2007. Back-analysis of rock mass strength parameters using AE monitoring data. Int. J. Rock Mech. Min. Sci. 44, 538–549. Chen, F., Wang, L., Zhang, W., 2019. Reliability assessment on stability of tunneling perpendicularly beneath an existing tunnel considering spatial variabilities of rock mass properties. Tunn. Undergr. Space Technol. 88 (2), 276–289. https://doi.org/10.
6
Tunnelling and Underground Space Technology xxx (xxxx) xxxx
C. Wu, et al.
Lyness, J.N., Moler, C.B., 1967. Numerical differentiation of analytic functions. SIAM J. Numer. Anal. 4 (2), 202–210. Ma, G.W., Wang, H.D., Fan, L.F., Wang, B., 2017. Simulation of two-phase flow in horizontal fracture networks with numerical manifold method. Adv. Water Resour. 108, 293–309. Maleki, M.R., 2018. Groundwater Seepage Rate (GSR); a new method for prediction of groundwater inflow into jointed rock tunnels. Tunn. Undergr. Space Technol. 71, 505–517. Peck, R.B., 1969. Advantages and limitations of the observational method in applied soil mechanics. Geotechnique 19 (2), 171–187. Qi, C., Fourie, A., 2018. A real-time back-analysis technique to infer rheological parameters from field monitoring. Rock Mech. Rock Eng. 51, 3029–3043. Qi, C., Fourie, A., Zhao, X., 2018. Back-analysis method for stope displacements using gradient-boosted regression tree and firefly algorithm. J. Comput. Civil Eng. 32 (5), 04018031. Sakurai, S., Takeuchi, K., 1983. Back analysis of measured displacements of tunnels. Rock Mech. Rock Eng. 16 (3), 173–180. Sun, Y., Jiang, Q., Yin, T., Zhou, C., 2018. A back-analysis method using an intelligent multi-objective optimization for predicting slope deformation induced by excavation. Eng. Geol. 239, 214–228. Tang, Y.G., Kung, G.T.C., 2009. Application of nonlinear optimization technique to back analyses of deep excavation. Comput. Geotech. 36 (1–2), 276–290. Vardakos, S., Gutierrez, M., Xia, C., 2016. Back-analysis of tunnel response from field monitoring using simulated annealing. Rock Mech. Rock Eng. 49 (12), 4833–4852. Wu, C., Chen, Q., Basack, S., Karekal, S., 2018. Laboratory investigation on rheological properties of greenschist considering anisotropy under multi-stage compressive creep condition. J. Struct. Geol. 114, 111–120. Wu, C., Chu, J., Wu, S., Hong, Y., 2019. 3D characterization of the microbially induced carbonate precipitation within rock fracture and the resulted permeability reduction. Eng. Geol. 249, 23–30. https://doi.org/10.1016/j.enggeo.2018.12.017. Yang, T.H., Jia, P., Shi, W.H., Wang, P.T., Liu, H.L., Yu, Q.L., 2014. Seepage–stress coupled analysis on anisotropic characteristics of the fractured rock mass around roadway. Tunn. Undergr. Space Technol. 43, 11–19. Yang, F., Zhang, C., Zhou, H., Liu, N., Zhang, Y., Azhar, M.U., Dai, F., 2017. The long-term safety of a deeply buried soft rock tunnel lining under inside-to-outside seepage conditions. Tunn. Undergr. Space Technol. 67, 132–146. Yao, K., Xiao, H., Chen, D., Liu, Y., 2019. A direct assessment for the stiffness development of artificially cemented clay. Géotechnique 69 (8), 741–747. https://doi.org/ 10.1680/jgeot.18.T.010. Zhang, W.G., Goh, A.T.C., 2015. Regression models for estimating ultimate and serviceability limit states of underground rock caverns. Eng. Geol. 188, 68–76.
1016/j.tust.2019.03.013. Du, H., Chen, W., Zhu, Q., Liu, S., Zhou, J., 2018. Identification of weak peaks in X-ray fluorescence spectrum analysis based on the hybrid algorithm combining genetic and Levenberg Marquardt algorithm. Appl. Radiat. Isot. 141, 149–155. Espada, M., Lamas, L., 2017. Back analysis procedure for identification of anisotropic elastic parameters of overcored rock specimens. Rock Mech. Rock Eng. 50 (3), 513–527. Exadaktylos, G.E., 2001. On the constraints and relations of elastic constants of transversely isotropic geomaterials. Int. J. Rock Mech. Min. Sci. 38 (7), 941–956. Finno, R.J., Calvello, M., 2005. Supported excavations: observational method and inverse modeling. J. Geotech. Geoenviron. Eng. 131 (7), 826–836. Gao, X.W., Liu, D.D., Chen, P.C., 2002. Internal stresses in inelastic BEM using complexvariable differentiation. Comput. Mech. 28 (1), 40–46. Gao, X., Yan, E.C., Yeh, T.C.J., Wang, Y.L., Cai, J.S., Hao, Y.H., 2018. Sequential back analysis of spatial distribution of geomechanical properties around an unlined rock cavern. Comput. Geotech. 99, 177–190. Gens, A., Ledesma, A., Alonso, E.E., 1996. Estimation of parameters in geotechnical back analysis-II. application to a tunnel excavation problem. Comput. Geotech. 18 (1), 29–46. Goh, A.T.C., Zhang, W., Zhang, Y., Xiao, Y., Xiang, Y., 2018. Determination of earth pressure balance tunnel-related maximum surface settlement: a multivariate adaptive regression splines approach. Bull. Eng. Geol. Environ. 77 (2), 489–500. Gonzaga, G.G., Leite, M.H., Corthésy, R., 2008. Determination of anisotropic deformability parameters from a single standard rock specimen. Int. J. Rock Mech. Min. Sci. 45 (8), 1420–1438. Grenon, M., Caudal, P., Amoushahi, S., Turmel, D., Locat, J., 2017. Analysis of a large rock slope failure on the east wall of the LAB chrysotile mine in Canada: back analysis, impact of water infilling and mining activity. Rock Mech. Rock Eng. 50 (2), 403–418. Hakala, M., Kuula, H., Hudson, J.A., 2007. Estimating the transversely isotropic elastic intact rock properties for in situ stress measurement data reduction: a case study of the Olkiluoto mica gneiss, Finland. Int. J. Rock Mech. Min. Sci. 44 (1), 14–46. He, J., Chen, S.H., Shahrour, I., 2011. Back analysis of equivalent permeability tensor for fractured rock masses from packer tests. Rock Mech. Rock Eng. 44 (4), 491. Li, Y., Sun, S., Tang, C., 2019. Analytical prediction of the shear behaviour of rock joints with quantified waviness and unevenness through wavelet analysis. Rock Mech. Rock Eng. 1–13. https://doi.org/10.1007/s00603-019-01817-5. Louis, C., 1972. Rock hydraulics. In: Rock mechanics. Springer, Vienna, pp. 299–387. Luo, Y., Chen, J., Chen, Y., Diao, P., Qiao, X., 2018. Longitudinal deformation profile of a tunnel in weak rock mass by using the back analysis method. Tunn. Undergr. Space Technol. 71, 478–493.
7
Contents lists available at ScienceDirect
Tunnelling and Underground Space Technology journal homepage: www.elsevier.com/locate/tust
A modified optimization algorithm for back analysis of properties for coupled stress-seepage field problems ⁎
Chuangzhou Wua,b, Yi Honga, Qingsheng Chenc, , Shivakumar Karekald a
College of Civil Engineering and Architecture, Zhejiang University, Hangzhou 310058, China Department of Geophysics, Kangwon National University, Chuncheon, Gangwon-do 200701, Republic of Korea c Department of Civil and Environmental Engineering, National University of Singapore, 21 Lower Kent Ridge Road, 119077, Singapore d School of Civil, Mining & Environmental Engineering, University of Wollongong, Wollongong City, NSW 2522, Australia b
A R T I C LE I N FO
A B S T R A C T
Keywords: Back analysis Coupled stress and seepage fields Levenberg-Marquardt optimization Complex-variable-differentiation method Parameters estimation
The study is aimed at developing a new algorithm for back analysis of rock mass parameters based on the observed displacements of tunnel excavations under coupled stress-seepage field problems which are generally encountered in hydraulic tunnel projects. The back-analysis algorithm is developed by incorporating the Levenberg-Marquardt optimization technique with complex-variable differentiation method. Additional auxiliary technique is also incorporated to enhance the convergence and stability of the proposed algorithm during the back-estimation of multiple rock mass parameters. A hypothetical hydraulic tunnel case was used for testing and validating the proposed algorithm by incorporating the method in a finite element code in which multiple rock mass parameters (such as modulus, permeability, and in situ stress) were treated as target parameters. Results show that the multiple rock mass parameters can be accurately and efficiently estimated by back analysis using a newly developed algorithm for coupled stress and seepage fields encountered in tunneling. The proposed algorithm can be used for predicting excavation behavior, particularly, the stress-induced deformations at subsequent stages of tunnel excavation under coupled multiple fields (e.g. stress and seepage fields).
1. Introduction Stress-seepage coupling analysis is commonly undertaken in geotechnical engineering especially in the design of hydraulic tunnels. Accurate estimation of rock mass properties is necessary for design of geotechnical structures. However, reliable estimation of rock mass properties is difficult in the laboratory. It is not only expensive but also several errors associated with disturbances in sample collection, testing and scale effects of the rock mass which make it difficult to estimate true properties of rock mass. In couple problems involving stress and seepage, it is all the more difficult to assess suitable properties of rock mass (Exadaktylos, 2001; Finno and Calvello, 2005; Gens et al., 1996; Gonzaga et al., 2008; Hakala et al., 2007 and Du et al., 2018; Li et al., 2019; Maleki, 2018; Vardakos et al., 2016; Yang et al., 2017). Furthermore, it is also time-consuming and economically undesirable to determine the mechanical and hydraulic properties of anisotropic rock mass especially in deep and interaction of the multi-field (e.g. stress and seepage coupled field) (Chen et al., 2019; Gonzaga et al., 2008; Hakala et al., 2007). Therefore, the back analysis method based on field observations, as suggested by Peck (1969), becomes important in
⁎
determining the representative and/or dominant geomaterial parameters in practice (Zhang et al., 2015; Goh et al., 2018). The backestimated parameters can then be used in the subsequent predictions of behavior or rock mass under different geometry and boundary conditions. Therefore, it is of great significance to develop an efficient and accurate approach for the parameter determination of fractured rock mass (Exadaktylos, 2001; Gao et al., 2018; Li et al., 2019; Luo et al., 2018; Sun et al., 2018; Tang and Kung, 2009; Zhang and Goh, 2015). The application of the back analysis technique often consists of two procedures: a forward procedure and a backward procedure (Qi and Fourie 2018; Qi et al., 2018). Optimization methods thus utilized during the backward procedure which can be classified into two categories: that base on the artificial intelligence (AI) theory (Qi and Fourie 2018; Qi et al. 2018) the classic optimization method. The optimization methods have been applied to calibrate the geomaterial parameters through back analyses of field observations over the past two decades (Sakurai and Takeuchi, 1983, Arai et al., 1983; Tang and Kung, 2009; Zhang et al., 2015; Wu et al., 2019; Yao et al., 2019). As the optimization method suffers from cancellation errors and slow convergence rates in addition to the likelihood of becoming trapped in local optima.
Corresponding author. E-mail addresses: [email protected] (C. Wu), [email protected] (Y. Hong), [email protected] (Q. Chen), [email protected] (S. Karekal).
https://doi.org/10.1016/j.tust.2019.103040 Received 2 January 2019; Received in revised form 14 March 2019; Accepted 12 July 2019 0886-7798/ © 2019 Elsevier Ltd. All rights reserved.
Please cite this article as: Chuangzhou Wu, et al., Tunnelling and Underground Space Technology, https://doi.org/10.1016/j.tust.2019.103040
Tunnelling and Underground Space Technology xxx (xxxx) xxxx
C. Wu, et al.
iteratively in the optimization process in search of global optimized values. The iterative procedure for modified L-M method used in this study can be expressed as
Thus, the complex-variable differentiation method was developed to evaluate derivatives of real functions in optimization methods (Lyness and Moler, 1967; Gao, 2002). This method is efficient since it has neither cancellation errors nor step-size-dependent problems compared to normal optimization methods. However, in dealing with complex coupled problems (i.e Stressseepage issues or Stress-seepage-thermal issues) involving multiple parameters, the optimization method often fail to converge and results in inaccurate values (Cai et al., 2007). In this study, a modified algorithm is proposed by incorporating the auxiliary techniques in the Levenberg-Marquardt (here mentioned as “L-M”) optimization technique to enhance the convergence and stability of the algorithm in dealing with coupled problems involving multiple fields such as stress-seepage coupling issues and Stress-seepage-thermal coupling issues. The proposed new method is then incorporated in a finite element code. A series of hypothetical excavation cases were used to test and validate the applicability of the proposed method for various scenarios of multiple parameter groups of coupled stress-seepage field problems.
x (k + 1) = x k + α (k + 1) Sk
(2)
Sk = (DkT Dk + μI )−1DkT (um (x k ) − ui*)
(3)
∂u ⎡ ∂x11 ⋯ ⎢ Dk = ⎢ ⋮ ⋱ ⎢ ∂um ⋯ ⎣ ∂x1
(4)
The modified optimization method was developed for the back analysis of excavation to back-figure multiple-parameters simultaneously based on observed deformations (e.g., displacement of the tunnel). The modified algorithm is proposed by incorporating the auxiliary techniques in the L-M optimization technique to enhance the convergence and stability of the algorithm. In mathematics and computing, the Levenberg-Marquardt algorithm (LMA), also known as the damped least-squares (DLS) method, is used to solve non-linear least squares problems. The LMA is more robust than the Gauss–Newton algorithm (GNA), which means that in many cases it finds a solution even if it starts very far off the final minimum. The improved optimization method is capable of conducting multi-parameters back analyses simultaneously of coupled stress and seepage field, which is described in the following sections. Note that such back-figured rock parameters should be regarded as the equivalent rock parameters, not real rock parameters (Tang and Kung, 2009; Zhang, et al., 2015; Espada and Lamas, 2017). Specifically, the deformation behavior of tunnel may be affected by a large number of factors such as the excavation width and depth, surround rock stiffness strength distribution, bolt arrangement, lining spacing and stiffness. The back-figured parameters are generally away from the real parameters and thus regarded as the equivalent parameters.
2.3. Convergence criteria In the proposed study, the objective function need to satisfy the following convergence criterion:
J (x k + ak Sk ) < J (x k )
(5)
The convergence criterion should satisfy the following conditions in order to terminate the process for obtaining reasonable estimates.
J (x k ) ≤ ε
(6)
whereε is tolerance which is a small positive number. 2.4. Auxiliary techniques dealing with complex issues To enhance convergence rate of the proposed optimization method, especially for the coupled cases such as stress and seepage field, the auxiliary technique is employed. In this study, a new relaxation factor was employed by substituting the step size factor, α k (marked as ‘SSF’ method) by step size matrix Y k (marked as SSM method), and Y k has the following expression
Yk
2.1. Objective function An objective function is necessary in a back analysis to determine the optimized rock mass parameters through iteration process. The objective function J (x *) used in this study is defined as
k ⎡ a1 ⋯ 0 ⎤ ⎢ = ⋮ aik ⋮ ⎥ ⎥ ⎢ ⎢ 0 ⋯ amk ⎥ ⎦n ⎣
(7)
aik
is the step size factor for each type of parameters in coupled where (stress-seepage) case scenario. By introducing the step size matrix Y k into the iterative procedure (Eq. (2)), the optimum step size for each type of parameters can be determined and assigned accordingly, which would enhance the convergence rate of the optimization process. In contrast, for normal optimization method only one step size for all types of parameters are used (Tang and Kung, 2009).
m i=1
⎤ ⎥ ⋮ ⎥ ∂um ⎥ ∂xn ⎦ k
where k is the iteration number and S is a search direction. The step size factor in terms of the scalar quantity α k defines the moving distance in the search direction, S . Dk , denotes the sensitivity matrix. n is the number of the parameters. In this test, the complex-variable-differentiation method was employed to calculate the sensitivity matrix Dk rather than direct differentiation used in normal optimization algorithms which will give rise to the cancellation errors or step-size-dependent problems. The algorithm of complex-variable-differentiation method (CVDM) is presented in detail in following sub-sections.
2. Proposed hybrid optimization algorithm
J (x *) = minJ (x ) = min ∑ (ui (x ) − ui*)2
∂u1 ∂xn
(1)
where, m is the number of measurement points used for back-predicting the intended multiple parameters in an optimization analysis. ui* is the measured displacement and ui is the computed displacement using FEM. x is a parameter vector composed of the intended multiple parameters. Note that the multiple types of parameters such as modulus, poison’s ratio, permeability, in situ stress and water pressure can be simultaneously calculated based on observed displacement of the tunnel.
2.5. Complex-variable-differentiation method In order to improve the convergence rate of the optimization process and to improve the accuracy of back-estimated multiple parameters for couple multi-field, the complex-variable-differentiation method (CVDM) developed by Lyness and Moler (1967) was incorporated in L-M method. This would eliminate the cancellation error and heavily dependent step-size effect. In this CVDM, the variable x of a real function u (x ) is replaced by a complex one of x + ih , with the imaginary part h which is very small (usuallyh = 10−20 y). The function u (x + ih) can be expanded into
2.2. Search direction and line search An initial set of target variables, x0 (i.e., initial condition) is required in optimization algorithms, and then these variables are updated 2
Tunnelling and Underground Space Technology xxx (xxxx) xxxx
C. Wu, et al.
Fig. 1. 2D FEM model (a) Tunnel geometry, (b) Mesh around the tunnel and (c) Layout of convergence measurement.
matrix, see Eqs. (4)) can be obtained by dividing the Im [u (x + ih)] with h as
Table 1 Mechanical and hydraulic property of the rock materials. Rock mass
Self-weight (kN/m3)
Eh (GPa)
Ev (GPa)
kh (10−8 m/s)
kv (10−8 m/s)
Type 1 Type 2 Type 3
27.2 24.0 26.1
15 1.1 16
8.0 1.8 22
5.14 1.92 0.63
4.2 1.92 0.63
du Im [u (x + ih)] = dx h
where Im denotes the imaginary part. It is shown that CVDM is a very promising method since the derivatives require to be evaluated. It can be used as easy as the direct difference method and has neither cancellation errors nor step-size-dependent problems. In CVDM, the variables should be complex types, and the derivatives of an inverted parameter can be obtained by adding a small imaginary part of ih . The numerical method, such as the direct differentiation method, is the easiest way to compute the sensitivity coefficients of Eqs. (4).
Note: E – Elastic modulus of rock mass, k – Hydraulic conductivity, Subscripts of h and v – indicate the relative parameters in the horizontal and vertical directions representatively
Taylor’s series as
u (x + ih) = u (x ) + ih
du h2 d2u − + o (h3) dx 2 dx 2
(9)
2.6. Solution procedure
(8)
The computational steps based on the displacement can be
Since h is very small, derivatives of a real function (i.e. sensitivity Table 2 Case studies information. Target parameters
Case 1 Case 2 Case 3
Auxiliary technique
Modulus
Hydraulic Conductivity
Field stress
Eh1, Ev1, Eh2, Ev2, Eh3, Ev3 Eh1, Ev1 Eh1, Ev1, Eh2, Ev2
– k1h, k1v k2
ph, pv ph, pv
Note: Initial values setting: Moduli (E) as 1 GPa, Hydraulic conductivity (k) as 10 × 10−8 m/s and field stress (p) as 10 MPa. 3
SSF SSF SSM (compared with SSF)
Tunnelling and Underground Space Technology xxx (xxxx) xxxx
C. Wu, et al.
25
20
J16=2.6*10-13
Eh3
Eh1 Ev1
10
J5=4*10-4
Ev2 Eh2
J0=0.246
0
Initial values Values after 5 steps of iteration Values after 13 steps of iteration True values
20
J13=4.7*10-15
E1h
15
10
E1v J5=3.2*10-7
5
J0=1.6*10-3
0 0
-10 -10
0
10
20
30
5
10
Calculated permeability, k (10-8 m/s)
Initial values Values after 5 steps of iteration Values after 16 steps of iteration
pv ph
J16=2.6*10-13
20
Iteration process
J5=4*10-4 J0=0.246
10
J0=1.6*10-3
10
Iteration process
8
kh1
6
J5=3.2*10-7
kv1
J13=4.7*10-15
4
Initial values Values after 5 steps of iteration Values after 13 steps of iteration True values
2 0 0
2
4
6
8
10
12
Measured permeability, k (10-8 m/s)
(b)
0 0
10
20
30
40
50
50
Calculated loading pressure, p (MPa)
Measured loading pressure, p (MPa)
(b) Fig. 2. The iteration process of inverted parameters for case 1 (Eh1, Ev1, Eh2, Ev2, Eh3, Ev3, ph, and pv) with objective function value after different iterations.
summarized as follows: Step 1. Input measured displacements, u* and initialize target parameters xi. Step 2. Call the subroutines (except the Step 3 in the direct problem) to obtain back-estimated values of parameters xi. Step 3. Terminate the iterative procedure if the stopping criterion is achieved; otherwise, employing SSM method (Eqs. (7)) and then go to Step 2.
Initial values Values after 5 steps of iteration
40
Values after 13 steps of iteration True values
J13=4.7*10-15 J5=3.2*10-7
pv
30
ph 20
10
J0=1.6*10-3
Iteration process
Calculated loading pressure, p (MPa)
25
(a)
50
30
20
12
(a)
True values
15
Measured modulus, Em (GPa)
40
Measured Modulus, Em (GPa)
40
Iteration process
30
Calculated modulus, Ec (GPa)
Initial values Values after 5 steps of iteration Values after 16 steps of iteration Ev3 True values
Iteration process
Calculated Modulus, Ec (GPa)
40
0 0
10
20
30
40
50
Measured loading pressure, p (MPa) 3. FEM code and Stress-seepage coupling analysis
(c)
The proposed optimization method was compiled and incorporated into the FEM code on MATLAB platform. A hypothetical case is considered for testing the proposed algorithm. A tunnel is excavated in deep rock mass with three different types of stratum bounded by high in situ stress and seepage pressure (He et al., 2011; Yang et al., 2014; Ma et al., 2017; Grenon et al., 2017). In terms of the seepage-stress coupling analysis, the following assumptions should be confirmed. (i) initially, the rock mass is treated as an elastic, homogeneous and anisotropic; and (ii) seepage tensor is defined in the Representative Element
Fig. 3. The iteration process of inverted parameters for case 2: (a) Eh1 and Ev1 (b) ph, and pv (c) kh1 and kv1 with objective function value at different iteration number.
Volume (REV) of the rock mass (Yang et al., 2014). For orthotropic rock mass in plane problems, the elastic constitutive material model with no pore water pressure is given as by:
4
Tunnelling and Underground Space Technology xxx (xxxx) xxxx
C. Wu, et al.
10 E2v
E1v
E1h
8
6
10
Using SSM code after 10 iteration of SSF
Real value
Modulus, E (GPa)
k2
4
5 2
0
Initial value
0
5
10
15
20
25
( V = − (K y
u12 u23 u13 u14 u45 u15 u16 u67 u17
19.847 21.428 30.531 38.805 47.537 52.181 37.676 39.485 31.388
Relative errors(%)
Identified convergence
Relative errors(%)
19.885 21.473 30.622 38.917 47.716 52.358 37.723 39.586 31.472
0.19 0.21 0.30 0.29 0.38 0.34 0.12 0.26 0.27
19.856 21.445 30.559 38.834 47.575 52.227 37.697 39.506 31.406
0.09 0.08 0.09 0.07 0.08 0.09 0.06 0.05 0.06
2
∂H ∂y
(12)
(13)
4.1. Model settings
Identified convergence
2
+ Kyy
4. Validation of the proposed optimization algorithm
0
A two-dimensional plane strain finite element model was considered for back predicting the properties of rock mass surrounding a tunnel excavation subjected to stress and seepage fields (coupled fields). A horse-shoe shaped tunnel with a height of 8 m and a width of 6 m is considered. The boundary is set to about 4 times of width of the tunnel. The mechanical and hydraulic properties of the rock mass were collected from a practical project completed in China (Wu et al., 2018). The rock mass was subjected to in-situ stress with major in-plane principal stress (pv) of 30 MPa in vertical direction and the minor inplane principal stress ph of 24 MPa in horizontal direction. The seepage pressure of the rock mass at the boundary is considered as 10 MPa and no water head inside of the tunnel. The seepage direction is from outside to inside of the tunnel. The constitutive model of the rock mass in this study is the linear elastic model. The model setting and meshes are shown in Fig. 1a and b, respectively. Several points surrounding the tunnel surface are selected as measuring points as shown in Fig. 1c, thus the convergence of the tunnel face could be represented by the relative displacements between the measuring points uij. Here, nine convergence measurements are adopted for the following studies. The rock mass mechanical parameters are listed in Table 1, which are further adopted as target values in the back analysis. The tunnel is driven through stratified rock mass consists of three different rock types, characterizing anisotropic mechanical behavior. In the following section, these mechanical parameters and the field stress are compared with those calculated from the proposed back analysis approach with nine convergence measurements. To illustrate the application of the proposed method, three case studies are carried out as shown in Table 2. In general, Case 1 and Case 2 are adopted to assess the efficiency and accuracy of the proposed optimization method through the determination of different types of parameters, such as modulus, permeability, in situ stress, etc. individually (separately) and combinedly (coupled manner). On the other hand, it should be noted that the 2nd type of rock mass has much lower modulus and larger hydraulic conductivity than the other two types which would be essential for the tunnel stability analysis. Therefore, Case 3 is used to present a more crucial case which is dealing with nonconvergence situations due to the high variation of sensitivity of different type of parameters which need be back-estimated
SSM
v v ⎡ 1 − v31 − ⎛ E12 + E31 ⎞ E3 3⎠ ⎢ E1 ⎝ 1 ε11 ⎡ ⎤ ⎢ 2 2 v v v 1 − E31 ⎢ ε22 ⎥ = ⎢− ⎛ E12 + E31 ⎞ E1 3⎠ 3 ⎣ ε12 ⎦ ⎢ ⎝ 1 ⎢ 0 0 ⎢ ⎣
∂H xy ∂x
)⎫⎪ )⎬⎪⎭
where Kf is the current groundwater hydraulic conductivity; K 0 is the initial hydraulic conductivity; σ is the normal stress on the joint plane, and α is the coupling parameter that reflect the influence of stress.
Table 3 Comparisons of the tunnel convergence in case 3. SSF
∂H ∂y
Kf = K 0 e−ασ
Fig. 4. Improved convergence of inverting parameters for case 3 by using SSM after termination of iteration using SSF code at iteration n = 10.
Measured values
+ Kxy
The rock mass permeability is controlled by the normal stress. The stress is coupled to seepage according to Eq. (11) by Louis (1972).
Iteration number, n
Convergence
∂H ∂x
Vx = − Kxx Permeability, k (10-8 m/s)
15
E2h
⎤ ⎥ ⎥ ⎡ σ11 ⎤ 0 ⎥ ⎢ σ22 ⎥ ⎥ ⎣ σ12 ⎦ 2(1 + v12 ) ⎥ ⎥ E1 ⎦ 0
(10)
where E1 and E3 are the Young’s moduli and v31 and v12 are Poisson’s ratio The relationship between the hydraulic pressure and the hydraulic conductivity, based on the Biot’s theory (Biot, 1941), is expressed as
kij ∇2 P = 0
(11)
where kij is the hydraulic conductivity, P is the hydraulic pressure. The equation of hydraulic conductivity for a plane strain model can be expressed as:
Table 4 Analysis results of the hypothetical case 3 using SSF and SSM code. Parameter type
Eh1 (GPa) Ev1 (GPa) Eh2 (GPa) Ev2 (GPa) k2 (10−8 m/s)
True values
15 8 1.1 1.8 1.92
Initial values
1 1 1 1 1
SSF
SSM
Convergence values
Relative error (%)
Convergence values
Relative error (%)
14.6 7.7 1.12 1.86 2.41
2.7 3.9 1.8 5.4 25.7
15 8.0 1.1 1.8 1.93
< 0.01 < 0.01 < 0.01 < 0.01 0.42
5
Tunnelling and Underground Space Technology xxx (xxxx) xxxx
C. Wu, et al.
such a case, the residuals of the computed objective function for this case are as low as 8 × 10−12, meaning that the solution is converged for the back-estimated two kinds of parameters. The results of case 3 suggest that the proposed SSM method improved the optimization process and it seems to be more suitable for the back analysis when multiple parameters are encountered. The case 3 is a theoretical case study in which the exact rock mass properties are known, which allowed to realize that the expected solution in case 3 did not converge with the SSF method, and so it was necessary to apply the SSM. As far as field application was concerned, the displacements obtained in case 3 with the SSF method shows that the SSF method is quite satisfactory for back-analysis of rock mass properties including coupling issues the relative errors of which is lower than 1%. However the proposed method using SSM method shows higher precision for complex multiple parameters for multiple field coupling issues (i.e. seepage field, stress field and thermal field).” Based on the results of the case study, the proposed optimization model appears to be sound and suitable for back analysis for the coupled seepage-stress problems. However, it should be noted that, further cross validation will be needed on the accuracy of the estimation of rock mass parameters for the field case with measured in-situ data in order to use the proposed optimization algorithm for engineering applications.
simultaneously. 4.2. Analysis of case 1 and case 2 In case 1, the back-calculated parameters were set as the moduli (Eh1, Ev1, Eh2, Ev2, Eh3, Ev3) and the field stresses as (ph, pv). The initial values of moduli are set as 1.0 GPa, while those of field stresses are set as 10.0 MPa. Fig. 2 presents the convergence process of moduli and filed stresses after the 0rd, 5th and 16th iteration with their objective functions. In general, the objective function reduces with the increase in the number of iterations. It can be seen that the back-calculated parameters approaches the target values after the 5th iteration as the residual of the objective function becomes J5 = 4 × 10−4. The analysis stopped at the 16th iteration with J16 = 2.6 × 10−13. In case 2, the back-calculated parameters were set as the moduli, the hydraulic conductivity and the field stresses. Here, the parameters of rock type 1 are studied in this case, such as the moduli (Eh1, Ev1), the hydraulic conductivity (k1h, k1v) and the field stresses (ph, pv). As shown in Fig. 3a and Fig. 3b, the back-calculated hydraulic conductivity and field stresses are close to the targets values after 5 iterations with J5 = 3.2 × 10−7. However, the back-calculated moduli needed several iterations more to achieve the target value. Finally, all the back-calculated parameters could be close to the target values after 13th iteration (J13 = 4.7 × 10−15). The results from case 1 and case 2 indicate that the proposed backanalysis scheme is useful in identifying the multiple types of parameters such as modulus, field stress and hydraulic conductivity in dealing with coupling issues. More iterations for case 1 are needed in search of global solution than that of case 2, which is perhaps due to initial conditions of number of moduli parameters which would have allowed the solution to get caught into local minimum, thus increasing the iteration number until it reaches the final solution or final minimum. In case 2, the proposed optimization model seems effective in dealing with three types of parameters simultaneously by considering coupling effect of the hydraulic parameters and in-situ stresses.”
5. Concluding remarks A modified optimization technique is proposed, tested and validated to back predict the multiple properties simultaneously for the coupled stress-seepage related problems based on the observed displacements. The proposed optimization algorithm can also be used for coupled problems in anisotropic rock mass material models. Some of the main derived conclusions are: (i) The proposed optimization algorithm shows better precision in estimating parameters such as modulus, permeability and in situ stress, simultaneously, in coupled problems. In case 2, for different types of parameters, the proposed optimization model seems more sensitive to the hydraulic parameters and in-situ stresses than modulus. (ii) In dealing with the complex problem such as the target parameters which are associated with coupling issues of anisotropic materials (see case 3), an auxiliary technique by substituting the step size factor (SSF method) with factor matrix (SSM method) can be employed to improve the convergence, which is useful in solving the large convergence problems with low relative errors.
4.3. Analysis of case 3 In case 3, the back-calculated parameters were set as the moduli of both rock types 1 and 2 and the hydraulic conductivity of rock mass type 2. The SSF method is firstly adopted to back-analyze the target parameters. Totally 10 iterations are carried out using SSF method. As shown in Fig. 4, the moduli (Eh1, Ev1, Eh2, Ev2) are approaching to the target values as number of iterations increases. However, the hydraulic conductivity k2 could not be well identified at the 10th iteration and the relative error between the identified value and measured value is as large as 25.7%. Looking at the tunnel convergences shown in Table 3, they are also not well identified. This might be caused by the high deformability and high hydraulic conductivity of rock type 2. In the model, the convergence measurement u15 becomes as high as 52.2 mm which is much larger compared to other locations. The other two measurement u14 and u45 are also relatively large which could affect the accuracy of the calculated k2 value. To improve the identification of hydraulic conductivity k2, the auxiliary technique (SSM method) is further adopted by introducing a matrix factor to substitute the step size factor. As shown in Fig. 4, k2 achieved the target value progressively using SSM method. The error between the identified value and the measured one is as low as 0.42% after another ten iterations. At the same time, there is no influences on the identified modulus (Eh1, Ev1, Eh2, Ev2) which are approaching the target values as well. The increasing accuracy to identify the rock properties are listed in the Table 4 that the errors of moduli become less than 0.1% after SSM method and hydraulic conductivity becomes as low as 0.42%. As shown in Table 3, the accuracies of the tunnel convergence are also improved after using SSM method and the errors between the identified and measured values become less than 0.1%. In
Acknowledgement This work was financially supported by State Key Laboratory for GeoMechanics and Deep Underground Engineering, China University of Mining & Technology, China (Grant Nos: SKLGDUEK1521 and SKLGDUEK1819), National Natural Science Foundation of China (Grant No: 51509219), and Basic Science Research Program through the National Research Foundation of Korea (NRF) funded by the Ministry of Education (2019R1A6A1A03033167). References Arai, K., Ohta, H., Yasui, T., 1983. Simple optimization techniques for evaluating deformation moduli from field observations. Soils Found. 23 (1), 107–113. Biot, M.A., 1941. General theory of three-dimensional consolidation. J. Appl. Phys. 12 (2), 155–164. Cai, M., Morioka, H., Kaiser, P.K., Tasaka, Y., Kurose, H., Minami, M., Maejima, T., 2007. Back-analysis of rock mass strength parameters using AE monitoring data. Int. J. Rock Mech. Min. Sci. 44, 538–549. Chen, F., Wang, L., Zhang, W., 2019. Reliability assessment on stability of tunneling perpendicularly beneath an existing tunnel considering spatial variabilities of rock mass properties. Tunn. Undergr. Space Technol. 88 (2), 276–289. https://doi.org/10.
6
Tunnelling and Underground Space Technology xxx (xxxx) xxxx
C. Wu, et al.
Lyness, J.N., Moler, C.B., 1967. Numerical differentiation of analytic functions. SIAM J. Numer. Anal. 4 (2), 202–210. Ma, G.W., Wang, H.D., Fan, L.F., Wang, B., 2017. Simulation of two-phase flow in horizontal fracture networks with numerical manifold method. Adv. Water Resour. 108, 293–309. Maleki, M.R., 2018. Groundwater Seepage Rate (GSR); a new method for prediction of groundwater inflow into jointed rock tunnels. Tunn. Undergr. Space Technol. 71, 505–517. Peck, R.B., 1969. Advantages and limitations of the observational method in applied soil mechanics. Geotechnique 19 (2), 171–187. Qi, C., Fourie, A., 2018. A real-time back-analysis technique to infer rheological parameters from field monitoring. Rock Mech. Rock Eng. 51, 3029–3043. Qi, C., Fourie, A., Zhao, X., 2018. Back-analysis method for stope displacements using gradient-boosted regression tree and firefly algorithm. J. Comput. Civil Eng. 32 (5), 04018031. Sakurai, S., Takeuchi, K., 1983. Back analysis of measured displacements of tunnels. Rock Mech. Rock Eng. 16 (3), 173–180. Sun, Y., Jiang, Q., Yin, T., Zhou, C., 2018. A back-analysis method using an intelligent multi-objective optimization for predicting slope deformation induced by excavation. Eng. Geol. 239, 214–228. Tang, Y.G., Kung, G.T.C., 2009. Application of nonlinear optimization technique to back analyses of deep excavation. Comput. Geotech. 36 (1–2), 276–290. Vardakos, S., Gutierrez, M., Xia, C., 2016. Back-analysis of tunnel response from field monitoring using simulated annealing. Rock Mech. Rock Eng. 49 (12), 4833–4852. Wu, C., Chen, Q., Basack, S., Karekal, S., 2018. Laboratory investigation on rheological properties of greenschist considering anisotropy under multi-stage compressive creep condition. J. Struct. Geol. 114, 111–120. Wu, C., Chu, J., Wu, S., Hong, Y., 2019. 3D characterization of the microbially induced carbonate precipitation within rock fracture and the resulted permeability reduction. Eng. Geol. 249, 23–30. https://doi.org/10.1016/j.enggeo.2018.12.017. Yang, T.H., Jia, P., Shi, W.H., Wang, P.T., Liu, H.L., Yu, Q.L., 2014. Seepage–stress coupled analysis on anisotropic characteristics of the fractured rock mass around roadway. Tunn. Undergr. Space Technol. 43, 11–19. Yang, F., Zhang, C., Zhou, H., Liu, N., Zhang, Y., Azhar, M.U., Dai, F., 2017. The long-term safety of a deeply buried soft rock tunnel lining under inside-to-outside seepage conditions. Tunn. Undergr. Space Technol. 67, 132–146. Yao, K., Xiao, H., Chen, D., Liu, Y., 2019. A direct assessment for the stiffness development of artificially cemented clay. Géotechnique 69 (8), 741–747. https://doi.org/ 10.1680/jgeot.18.T.010. Zhang, W.G., Goh, A.T.C., 2015. Regression models for estimating ultimate and serviceability limit states of underground rock caverns. Eng. Geol. 188, 68–76.
1016/j.tust.2019.03.013. Du, H., Chen, W., Zhu, Q., Liu, S., Zhou, J., 2018. Identification of weak peaks in X-ray fluorescence spectrum analysis based on the hybrid algorithm combining genetic and Levenberg Marquardt algorithm. Appl. Radiat. Isot. 141, 149–155. Espada, M., Lamas, L., 2017. Back analysis procedure for identification of anisotropic elastic parameters of overcored rock specimens. Rock Mech. Rock Eng. 50 (3), 513–527. Exadaktylos, G.E., 2001. On the constraints and relations of elastic constants of transversely isotropic geomaterials. Int. J. Rock Mech. Min. Sci. 38 (7), 941–956. Finno, R.J., Calvello, M., 2005. Supported excavations: observational method and inverse modeling. J. Geotech. Geoenviron. Eng. 131 (7), 826–836. Gao, X.W., Liu, D.D., Chen, P.C., 2002. Internal stresses in inelastic BEM using complexvariable differentiation. Comput. Mech. 28 (1), 40–46. Gao, X., Yan, E.C., Yeh, T.C.J., Wang, Y.L., Cai, J.S., Hao, Y.H., 2018. Sequential back analysis of spatial distribution of geomechanical properties around an unlined rock cavern. Comput. Geotech. 99, 177–190. Gens, A., Ledesma, A., Alonso, E.E., 1996. Estimation of parameters in geotechnical back analysis-II. application to a tunnel excavation problem. Comput. Geotech. 18 (1), 29–46. Goh, A.T.C., Zhang, W., Zhang, Y., Xiao, Y., Xiang, Y., 2018. Determination of earth pressure balance tunnel-related maximum surface settlement: a multivariate adaptive regression splines approach. Bull. Eng. Geol. Environ. 77 (2), 489–500. Gonzaga, G.G., Leite, M.H., Corthésy, R., 2008. Determination of anisotropic deformability parameters from a single standard rock specimen. Int. J. Rock Mech. Min. Sci. 45 (8), 1420–1438. Grenon, M., Caudal, P., Amoushahi, S., Turmel, D., Locat, J., 2017. Analysis of a large rock slope failure on the east wall of the LAB chrysotile mine in Canada: back analysis, impact of water infilling and mining activity. Rock Mech. Rock Eng. 50 (2), 403–418. Hakala, M., Kuula, H., Hudson, J.A., 2007. Estimating the transversely isotropic elastic intact rock properties for in situ stress measurement data reduction: a case study of the Olkiluoto mica gneiss, Finland. Int. J. Rock Mech. Min. Sci. 44 (1), 14–46. He, J., Chen, S.H., Shahrour, I., 2011. Back analysis of equivalent permeability tensor for fractured rock masses from packer tests. Rock Mech. Rock Eng. 44 (4), 491. Li, Y., Sun, S., Tang, C., 2019. Analytical prediction of the shear behaviour of rock joints with quantified waviness and unevenness through wavelet analysis. Rock Mech. Rock Eng. 1–13. https://doi.org/10.1007/s00603-019-01817-5. Louis, C., 1972. Rock hydraulics. In: Rock mechanics. Springer, Vienna, pp. 299–387. Luo, Y., Chen, J., Chen, Y., Diao, P., Qiao, X., 2018. Longitudinal deformation profile of a tunnel in weak rock mass by using the back analysis method. Tunn. Undergr. Space Technol. 71, 478–493.
7