Kriging-based optimization design of deep tunnel in the rheological Burger rock

The principal purpose of this work consists in optimizing the support system of a deep tunnel accounting for the uncertainty of the time-dependent behaviour of the surrounding rock, which is described by the rheological Burger law. The stochastic approach is chosen for this aim. On one hand the Quantile Monte Carlo (QMC) simulation is used to determine the optimal design variables (i.e., the thickness of two liners). On the other hand, the well-known Kriging metamodeling technique is undertaken to approximate the limit state function in the augmented reliability space (i.e., the tensor product between the random variable space and the design variable space). The adopted optimization process allows to derive the optimal tunnel support that verifies two failure modes, namely the support capacity criterion and the maximum tunnel convergence.


Introduction
The optimization in the design process of the deep tunnel with respect to the uncertainty in input mechanical parameters is one of the most critical issues in the rock engineering field. The conservative design of the support element has been usually conducted based on the deterministic approach using the factor of safety of different parameters involved. Although the robustness thanks to its simplicity in the engineering application, such an approach does not permit us to understand the effect of the uncertainty in ground properties, particularly in the case of highly heterogeneous formations and can induce an overestimated result of the support elements.
In the last decade, many scholars attempt to account for the uncertainty in the analysis and the design process of tunnel by using rational statistical treatment. Especially, the performance of some well-known probabilistic methods such as the direct Monte Carlo Simulation (MCS), the Response Surface Method (RSM) in combination with the First or Second-Order Reliability Method (FORM/SORM) or the Kriging metamodeling technique in the optimization of tunnel supports was demonstrated in many contributions. However, this procedure, known as the Reliability-based optimization design (RBOD) was mostly limited in case of time-independent behaviour, while many experimental studies in the field revealed a significant time-dependent response of a large range of rock formations. Due to the timedependent characteristics of the studied problem, the application of Monte Carlo Simulations (MCS) becomes time-consuming and its combination with the metamodeling technique such as the Kriging method is preferred. This latter will be chosen in this work to optimize the thickness of the tunnel support system constructed in the rheological rock which verifies at the same time two failure modes, namely the support capacity criterion and the maximum tunnel convergence.

Deterministic problem
Our considered problem consists of a deep tunnel excavated in a viscoelastic Burger rock. The Burger model is characterized by four parameters (GM, GK, M, K) which present respectively the spring and dashpot of the Maxwell and Kelvin elements. We assume that the tunnel has a circular section (with radius R) and is subjected to a hydrostatic stress ( 0 h p ) at the far field. The tunnel is supported by two liners installed at the same elapsed time (t0) after the instantaneous excavation of tunnel. This tunnel support system composes of a compressible outer liner and a concrete inner one. The behaviour of two liners is linear elastic and characterized by the corresponding elastic moduli E1, 1 and E2, 2. Their thicknesses are respectively d1 and d2. Especially, thanks to its low elastic moduli (E1 << E2, 10), the flexible liner, made from a compressible material, can absorb the convergence at the excavation wall due to the Burger viscoelastic behaviour of the rock and hence reduces the transmission in stress on the final concrete support layer. Furthermore, a perfect adhesive interfaces of rock/outer liner and outer/inner liner (i.e., no-slip condition) is adopted in this study. For the sake of clarity, the coupled flexible/concrete liner of deep tunnel in the Burger rock is illustrated in figure 1.
Analytical solutions have been derived to determine the response of tunnel excavated in the viscoelastic Burger rocks. For example, in the recent contribution of [1], the authors consider the sequential excavation and sequential installation of two liners of the deep tunnel surrounded by a Burger's rock. In comparison with the solution of [1], the deterministic problem in this work can be considered as a particular case. Thus, the closed-form solution of [1] can be straightforwardly applied and the interested reader can refer to this last contribution for the detail of the solution derivation.

Kriging-based reliability analysis
Recently, [1] performed the Monte Carlo Simulation (MCS) study to investigate the influence of uncertainty of the Burger parameters on the failure probability of the deep tunnel in this linear viscoelastic rock. Although the direct evaluations of the structure response in MCS provide an accurate estimation of the failure probability, an enormous number of trials is required in the MCS reduces its efficiency. Thus, in this study, the so-called Kriging metamodeling technique whose performance has been tested in [4,5] is preferential.
The Kriging metamodel aims at approximating the limit state function (LSF) G(X), which separates the safety and failure domains of structure response in the space of random variables X, by a Gaussian process ( ) [4,5,7] written in the form: The vectors β and k(X) in equation (1) designate correspondingly the vector of the regression coefficient and the vector of basis-functions of n elements while the first term k(X)β represents the mean value of the Gaussian process (i.e., the trend of the process). The second term Z(X) with the zero-mean IOP Conf. Series: Earth and Environmental Science 833 (2021) 012155 IOP Publishing doi:10.1088/1755-1315/833/1/012155 3 and variance ( , ′) = 2 ( , , ') ( 2 constant process variance) characterizes the stationary Gaussian process. The kernel function R(, X, X') presents the prescribed auto-correlation function with respect to the hyperparameter vector .
The unknown parameters 2 , , of the Kriging metamodel in equation (1) are calibrated from an optimization process by using the exact results (gathered in a vector y) of the LSF ( ) that are evaluated at different observation points (called the training points) of the Design of Experiment (DoE). For the iterative reconstruction of the Kriging surrogate, the initial (gathered in a matrix S) can be generated from the Latin Hypercube Sampling (LHS) technique. During the iterative process, the DoE is enriched by adding the new training points thanks to a so-called learning function. The final Kriging metamodel ( ) is obtained after the convergence criterion is satisfied.
Using the constructed metamodel ( ), the value of the LSF function can be predicted/interpolated for any realization of the random input vector X, which follows the Normal distribution. Following that, Kriging predictor provides the following mean and variance values ( ( ), 2 ( )) [4, 5, 6]: with: Then the well-known MCS can be applied by using Kriging predictor ( ) to interpolate the results of NMCS random samples based on which the failure probability can be estimated: In this work, we adopt the two failure modes as presented in [1] to investigate the stability at 100 years of the deep tunnel. Following that, we investigate the tunnel stability with respect to the allowable convergence and the allowable stress in the concrete liner. For the first failure mode, the convergence at the surface of tunnel is limited to a maximum value umax: For the second failure mode, we consider also that the tunnel instability occurs when the maximum equivalent stress qL2 in the concrete liner exceeds the allowable stress cl of the constitutive material: The random variable vector X in equations (5 and 6) presents the four parameters of Burger model (i.e., X = [GM, GK, M, K]). In addition, the time parameter t is omitted in equations (5 and 6) due to the monotone response of tunnel over time. In this case, the time-dependent failure probability is a priori known when it reaches the maximum at the chosen design lifetime of the tunnel (i.e., t=100 years). For the numerical applications, in   Figure 2 illustrates the long-term failure probability of tunnel with respect to the thickness of the coupled compressible/concrete liners. The reliability analysis is performed by assuming a coefficient of variation COV=25% of the random variables X (with X = [GM, GK, M,K]). The results in figure 2 highlight the strong dependence of the concrete liner stability on the thickness of both liners. However, the failure probability relating to the tunnel convergence seems only to be controlled by the thickness of the compressible layer d1. For the compressible layer with a thickness lower than 10(cm) we observe a very low probability of exceedance (largely smaller than 1%) of tunnel convergence with respect to the allowable value umax=45(mm). Correspondingly to this thickness range of the compressible liner (i.e., d1 < 0.1m), the concrete liner must support much more the transmitted loading due to the convergence of the surrounding rheological rock. In figure 3 are performed the iso-value contours of the probability of two failure modes on the plane of the thickness of liners which confirm the observations in figure 2.
Based on this iso-value contours, the thickness range of two liners can be potentially chosen corresponding to each pair of target probability of failure -the necessary parameters of the design process of the tunnel. For example, by fixing a target probability at 5% with the allowable stress cl=36 (MPa) in the concrete liner, a pretty large range we can be chosen for each pair of thickness for the tunnel support system. However, if the same target probability of 5% is adopted for the other failure mode (i.e., the stability of tunnel related to its convergence), the possibility to choose the thickness of both liners that verify two conditions of failure probability seems much smaller when the potential thickness range of two liners becomes much narrower ( figure 3). The procedure to find out the most appropriated thickness of liners based on the reliability analysis (also known as the reliability-based optimization design RBOD) will be shown and discussed in the next section.

Optimisation design
The RBOD of a tunnel consists in optimizing the thickness of the liners to ensure the stability of the underground structure with respect to the target failure probability. Mathematically, the RBOD method aims at minimizing the objective function f(d) while the probability of each failure mode must be smaller the chosen target failure probability [7]: In this equation, the design variables consist of the thickness of two liners (d = [d1, d2]) and the minimum thickness of the whole support system is chosen as the principal objective of the tunnel design. The methodology of RBOD has been largely developed in the last two decades and an overview of different approaches can be found in the interesting contributions of [8]. In general, the classical RBOD techniques can be classified on the two-level, mono-level and decoupled approaches [8]. The two-level approach uses two loops: the outer loop explores the design space by using a suitable optimization scheme whilst the reliability analysis in the random variable space is conducted in the inner loop. The time consuming of this approach is high when the reliability analysis must be repeated in each iteration of the outer loop (i.e., each value of design variables). The mono-level (or single loop) approach is proposed to reduce the cost of the previous approach. In this approach, the RBOD problem is solved by avoiding the reliability analysis and enforcing optimality conditions. It has shown however that this approach can fail to convergence when the starting point of the optimization problem is chosen far from the optimal solution [8]. In the decoupled approach, the RBOD is conducted in a sequence manner when an approximate deterministic optimization problem is solved using the information from the previous reliability analysis.
In their recent contribution, the authors of [9] proposed a quantile-based optimization procedure that uses the idea of augmented reliability space Z where all combination of deterministic or design/random variables can be considered (Z = [d, X], i.e., the tensor product of two spaces). Following this approach, the optimization problem defined in equation (8) (8) where is the allowable values of each failure mode (i.e., umax or cl in this work) and the corresponding quantile of each mode is defined as: In equation (9), Rj(Z) is the structure response (i.e., tunnel convergence or maximum equivalent stress in the concrete liner in our considered problem). By considering crude Monte Carlo sampling technique to estimate the quantile in each iteration for the current design variables, these authors proposed to approximate each LSF by a Kringing metamodel in the augmented reliability space Z. For the whole detail of the QMC approach, the interest reader can refer to [9]. In this work, the QMC calculation process are programmed in Matlab in which the toolbox DACE [4,10] is taken to construct the Kriging metamodel as previously mentioned. Following this QMC approach, two Kriging surrogates is iteratively built up in the augmented reliability space corresponding two chosen failure modes.
For the sake of clarity, in table 2 we highlight different steps of this Kriging-based Quantile Monte Carlo approach which was adapted for the design optimization of the deep tunnel defined in equation (8). Following this summarized procedure of Kriging-based QMC, the number of random samples of design variables N d =10 2 and random variables NMCS=10 4 of X are regenerated in each iteration of the optimization process which seem to be sufficient for the target probability equal or higher than 1% in practical design of tunnel. Thus, the total number of 10 6 sample in the augmented space Z is used in the interpolation to calculate each quantile in each iteration of the optimization process. Furthers, the new training points Z* to enrich each Kriging metamodel is chosen based on the U learning function [6] which is now adapted in the QMC approach (see also [9]) as follow: Knowing the quantile of two failure modes, solving the optimization problem with constraints in Eq. (23) is classical when different algorithms can be applied. In this work, the well-known SQP algorithm in the fmincom function that is available in Matlab is chosen for this purpose. Concerning the convergence condition, the stopping criterion written in equation (11) can be adopted (see [5])  Generate NMCS random samples of X to generate the augmented space Z 3 Corresponding to each value d (k) (k=1,2,..,N d ): interpolate the value of tunnel response for NMCS random samples in the augmented space Z (l) =( d (k) , X (l) ) (l=1,2,.., NMCS) by using the constructed Kriging predictor 4 Calculate the quantile of the failure mode 5 Identify the new training points Z* owing the minimum U value (equation 10) and update the DoE and hence the matrix S 6 Evaluate the tunnel response of the new training points (i.e., 1 ( * ) or 2 ( * )) and update y. 7 Update the Kriging metamodel (i.e., 1 ( ) or 2 ( )) from S and y Solving the optimization problem: Determine the optimal result of d that minimizes the cost function f(d) and verifies the constraints defined by two quantiles (equation 8) Check the convergence of the optimization problem (i.e., verify the stopping-criterion) defined in equation (11) If the convergence is verified: obtain the optimal value of d; otherwise, go to step 1 of the iteration process As an example, in figure 4 is captured the DoE of each Kriging surrogate in the thickness liners plan (d1 and d2) in which the iso-value contour of the failure probability calculated in the previous section is also highlighted for the illustration purpose. The corresponding target probability of each mode in this example is 1% for the first mode and 15% for the second mode. In addition, for this iterative optimization design, an initial DoE with 12 samples is generated by the LHS (Latin Hypercube Sampling) for each Kriging metamodel. Concerning the augmented reliability space, it is constructed as the product of 100 samples of design variables d and 10 4 samples of random variables X. For each iteration, the augmented reliability space is regenerated.   The results illustrated in figure 4 show that at each iteration, the adding new training points of each Kriging metamodel, which is selected from the U learning function, lie near the limit state function representing by the contour of the chosen target probability. With respect to the target failure probability of 15%, all the enriched points chosen in the second Kriging metamodel (i.e., the failure mode associated with the concrete liner stability) follow quite accurately the contour of this probability value. For the other failure mode, the enriched points seem more dispersed which can be explained by the small value of the chosen target probability of 1%.
In figure 5 are presented the evolution of the optimal design variables (i.e., the thickness of two liners). All the results of the thickness of each liner and hence the cost function (i.e., the total thickness of tunnel support system) attain the convergence criterion after 21 iterations. The results of the optimal thickness calculated by the Kriging-based QMC d1=21.14(cm) and d2=58.88(cm) are similar to the ones calculated with the two-level approach with d1=21.12(cm) and d2=58.41(cm). More precisely, this latter approach bases on the two loops: the result of the Kriging-based reliability analysis of the previous section is considered as the solution of the inner loop whilst the SQP algorithm in the fmincom function of Matlab is also chosen to solve the optimization problem in the outer loop.

Conclusions
In this work, the optimization design of a deep tunnel excavated in the viscoelastic Burger rock is conducted. The procedure allows accounting for the uncertainty of the time-dependent behavior of the rheological rock and consider at the same time two failure modes concerning the liner stability and convergence at the tunnel wall. To attain this purpose, the stochastic approach based on the Quantile Monte Carlo simulation that is combined with the Kriging metamodeling technique is chosen. The efficiency of the Kriging metamodel is investigated when we consider both problems: reliability analysis of the tunnel stability at 100 years and the optimization design of thickness of the tunnel support system. The study shows the robustness of the proposed RBOD for the deep tunnel surrounded by a rock with time-dependent behavior.