Application of coupled XFEM-BCQO in the structural optimization of a circular tunnel lining subjected to a ground motion

Abstract

A new structural optimization method of coupled extended finite element method and bound constrained quadratic optimization method (XFEM-BCQO) is adopted to quantify the optimum values of four design parameters for a circular tunnel lining when it is subjected to earthquakes. The parameters are: tunnel lining thickness, tunnel diameter, tunnel lining concrete modulus of elasticity and tunnel lining concrete density. Monte-Carlo sampling method is dedicated to construct the meta models so that to be used for the BCQO method using matlab codes. Numerical simulations of the tensile damage in the tunnel lining due to a real earthquake in the literature are created for three design cases. XFEM approach is used to show the cracks for the mentioned design cases. The results of the BCQO method for the maximum design case for the tunnel tensile damage was matching the results obtained from XFEM approach to a fair extent. The new coupled approach manifested a significant capability to predict the cracks and spalling of the tunnel lining concrete under the effects of dynamic earthquakes.

Keywords

ovaling deformation / monte carlo sampling / XFEM-BCQO / maximum principal stress

Introduction

The underground structures have been subjected to less amount of damage compared to overground structures. Probably, an appreciable rate of underground structures manifested significant damages in the last strong earthquakes. Generation of ground strain and the deformation due to wave propagation affect the responses of the tunnels. The wave types are controlling the motion of the soil which always are converted into a longitudinal and transverse components [¹]. The response of tunnels subjected to seismic loadings is complex to some extent until two decades ago which was not believed as a critical issue in the design process. Just latest experimental and numerical researches have made appreciable and important proceedings for the understanding of the mechanisms dominating the ground-tunnel interaction under seismic loadings. Commonly, the seismic response of the tunnel in soft soil is considered to be dominated by the soil around the tunnel and in the same time, the inertial load of the tunnel is neglected [²]. As a result, the alteration of the internal forces in the tunnel lining due to seismic vibration are commonly determined from the ground transient response. When earhquakes occur, deformation exists in the tunnel lining in the transverse and longitudinal sections which are resulted by the soil shear strain in both horizontal and vertical planes because of the non concurrent movement over the tunnel axis. The general methods of the design suggest separate mechanisms for the two deformation, without considering coupling effect [³–⁶].

Very extensive numerical simulations have been performed to determine stresses and deformations of tunnels in recent year. A research study has been performed by Ref. [⁷] using numerical analysis and simulation of a circular tunnel with uncertain three parameters related to the surrounding soil during earthquakes. They quantified the uncertainty of those parameters which are related to the responses of the tunnel through meta-models using an experimental method. Another study was conducted by Ref. [⁸] where they used finite element method to simulate the forces in the tunnel lining and the deformation in the soil for both Steinhaldenfeld and Heinenoord tunnels. Lekhnitskii [⁹] suggested a complicated variable-function method to accurately solve the stresses around a circular hole applying in-plane loading. Furthermore, Lu et al. [¹⁰] used the complex method and displayed the solutions for the stress in a circular tunnel. A two-dimensional elastic solution was founded by Ref. [¹¹] for a deep circular tunnel subjected to the far-field static loading [⁷–¹²].

Numerous computational models have been suggested and advanced since 1926 for the design of tunnels. An elastic continuum was adopted by Ref. [¹³] and they suggested the first analysis method. Another analytical solution was performed by Ref. [¹⁴] considering continuum models, where this solution adopts the elliptical deformation of the tunnel lining. A model for a circular tunnel in an elastic continuum was announced by Ref. [¹⁵] which has a geometrical nonlinearity. It is worthy to mention that the practical tunnel design method was suggested by Ref. [¹⁶]. Another continuum model also a bedded-beam model excluding ground pressure reduction at the crown was introduced for shallow tunnels having a ratio of tunnel depth to its diameter of C/D≤2. The soil-structure interaction was involved in the continuum model automatically, where in the bedded-beam model, this interaction can be secured by using bedding springs with a suitable stiffness. Despite the fact that this model takes into account the ground-tunnel interactions by using Winkler springs, only radial pressures are considered. Furthermore, a sophisticated model was introduces by Ref. [¹⁷] which includes the tangential pressures [¹⁸].

Many research works in the literature can be studied to significantly understand the most efficient methods that are related to the fracture modeling and they are considered similar methods to XFEM approach which is intended to be adopted in this study, for example, methods based on damage and crack detection [¹⁹–²⁹]. To assess the reliability and structural safety of tunnels, advanced computational methods are commonly applied which are capable of capturing the damages. These methods can be classified into continuous and discrete fracture approaches. Continuous approaches to fracture smear the crack over a certain width. They include non-local damage models [³⁰,³¹] and gradient models [³²,³³]. Also the introduction of a viscosity [³⁴–³⁸] smears the crack over a certain width. With the seminal work of Ref. [³⁹], phase field approaches have become another alternative to non-local and gradient models. Instead of relating the non-local damage to internal state variables, the damage is obtained by solving a partial differential equations. Phase field approaches have meanwhile been applied to numerous interesting problems [⁴⁰–⁴⁶] including fracture in thin shells. Discrete crack approaches include meshfree methods [⁴⁷–⁵²], extended finite element methods [⁵³–⁵⁵], smoothed extended finite element methods [⁵⁶–⁵⁸], phantom node methods [⁵⁹–⁶²]. One serious problem in these approaches is tracking the crack path especially in three dimensional and in dynamic fracture.

In this study we will use a combination of two methods as a new optimization technique supporting on a coupled extended finite element method-bound constrained quadratic optimization (XFEM-BCQO) in order to quantify the optimum values of four parameters of a tunnel lining constructed in a soft soil media which is subjected to a ground motion for 10 s duration. The responses of the tunnel system would be optimized regarding minimization of cracks in the tunnel lining.

Analytical solution of a circular tunnel

The easiest way to guess ovaling deformation in a circular tunnel is to suppose the deformations to be similar to free-field, as a result disregarding the tunnel-ground interaction. This hypothesis is suitable when the ovaling stiffness of the tunnel lining is equal to the stiffness of the surrounding soil. A continuous medium is considered for the circular tunnel-ground shearing (referring to as non-perforated ground) considering the absence of the tunnel (Fig. 1) where the diametric strain for a circular section of the tunnel would be calculated as Eq. (1) in Ref. [²]:

(1){\displaystyle \frac{{{\displaystyle \mathrm{\Delta }d}}_{\mathrm{f}\mathrm{r}\mathrm{e}\mathrm{e}-\mathrm{f}\mathrm{i}\mathrm{e}\mathrm{l}\mathrm{d}}}{d}}=\pm {\displaystyle \frac{{{\displaystyle \gamma }}_{\max }}{\mathrm{2}}}.

When the ovaling stiffness of the tunnel lining is too small compared to the surrounding soil, the tunnel diametric strain is calculated proposing an unlined tunnel as Eq. (2) in Ref. [²]. referred to as perforated ground:

(2){\displaystyle \frac{{{\displaystyle \mathrm{\Delta }d}}_{\mathrm{f}\mathrm{r}\mathrm{e}\mathrm{e}-\mathrm{f}\mathrm{i}\mathrm{e}\mathrm{l}\mathrm{d}}}{d}}=\pm {{\displaystyle \mathrm{2}\gamma }}_{\max }\left(\mathrm{1}-{{\displaystyle v}}_{\mathrm{m}}\right).

That {{\displaystyle \mathrm{\Delta }d}}_{\mathrm{f}\mathrm{r}\mathrm{e}\mathrm{e}-\mathrm{f}\mathrm{i}\mathrm{e}\mathrm{l}\mathrm{d}} is free-field diametric deflection in non-perforated ground, d is the diameter of the tunnel, {{\displaystyle \gamma }}_{\max } is maximum free-field shear strain of soil or rock medium, and {{\displaystyle v}}_{\mathrm{m}} is Poisson’s ratio of soil and this deformation is much greater in the case where the presence of the tunnel is included compared to the case where only the continuous ground deformation is assumed. Also, the maximum shear strain in the soil for both constant and variable shear strain may be calculated as Ref. [³]:

(3){{\displaystyle G}}_{\mathrm{m}}={\displaystyle \frac{{{\displaystyle E}}_{\mathrm{m}}}{\mathrm{2}(\mathrm{1}+{{\displaystyle v}}_{\mathrm{m}})}},

(4){{\displaystyle C}}_{\mathrm{m}}=\sqrt{{\displaystyle \frac{{{\displaystyle G}}_{\mathrm{m}}}{{{\displaystyle \rho }}_{\mathrm{m}}}}},

(5){{\displaystyle \gamma }}_{\max }={\displaystyle \frac{{{\displaystyle V}}_{\mathrm{s}}}{{{\displaystyle C}}_{\mathrm{m}}}},

where {{\displaystyle \rho }}_{\mathrm{m}}={{\displaystyle \gamma }}_{\mathrm{s}}/g, {{\displaystyle E}}_{\mathrm{m}} is Young’s modulus, {{\displaystyle G}}_{\mathrm{m}} is shear modulus, {{\displaystyle V}}_{\mathrm{s}} is peak particle velocity associated with S-waves, and {{\displaystyle C}}_{\mathrm{m}} is apparent velocity of S-wave propagation in soil. Full-slip assumption under simple shear strain, may cause important under speculation of the maximum thrust [²,³,⁶³,⁶⁴].

Earthquake amplitudes

A real time history from the literature with adjustment up to the requirements of the parameric analysis has been used for the duration of 10 s in the simulation of the tunnel responses when subjected to the ground motion due to an earthquake, see Figs. 2 and 3.

Finite element model

Figure 4 illustrates the 2D model of both the tunnel lining and the soil media generated in ABAQUS with dimensions to be prepared for simulating the earthquake events for the duration of 10 s both in vertical and horizontal directions.

The tunnel lining thickness, diameter, concrete modulus of elasticity, and concrete density are the selected parameters to predict the responses of the tunnel regarding the damage in the tunnel lining. The dimensions of the soil media are 100 m × 50 m. The tunnel lining material is made of reinforced concrete where the reinforcement effect is neglected for simplicity. Concrete damaged plasticity model used for damage detection analysis supporting on real data in Ref. [⁶⁵]. The XFEM approach has been used to detect the cracks and damages in the tunnel lining due to the 10 s duration earthquake supporting on the maximum principal stress output.

The density and elastic properties of the soil were assigned with different values for meta-modeling process. Mohr Coulomb plasticity model was used for the soil media to define its shear strengths at different effective stresses. In the initial step, two predefined fields were created for the soil, geostatic stress of 20601 Pa and void ratio of 1. A geostatic step of 1 s duration was generated, where the gravity loads of the model system were assigned. Two different boundary conditions were assigned for both the base and the sides of the soil media, by preventing the displacement in x-axis direction for both of them and providing two degrees of freedom for the sides in y-direction and rotation about z-axis for both of them. A static step with duration of 1e–10 s was created for the model system. A dynamic implicit step with 10 s duration was created to simulate the ground motion due to the earthquake by assigning acceleration boundary conditions in both x and y directions with magnitude of 4.905 m/s² . A tie constraint was assigned to model the surface contact between the soil media and the tunnel lining. The soil media and the tunnel lining were modeled with standard plain strain linear with reduced integration hourglass control elements CPE4R.

Monte Carlo sampling

Monte Carlo techniques are commonly used to perform uncertainty and sensitivity analysis. A key element of Monte Carlo method is the sampling of input parameters for the simulation, where the goal is to explore the entire input space with a reasonable sample size N. The sample size determines the computational cost of the analysis since N is equal to the required number of simulation runs. The unbiased mean and variance of the model output can be calculated by the following equations.

(6)\overline{y}={\displaystyle \frac{\mathrm{1}}{\mathrm{n}}}\mathrm{\sum }_{i=\mathrm{1}}^N{{\displaystyle y}}_i,

(7)Var(y)={\displaystyle \frac{\mathrm{1}}{N-\mathrm{1}}}\mathrm{\sum }_{i=\mathrm{1}}^N{{\displaystyle \left({{\displaystyle y}}_i-\overline{y}\right)}}^{\mathrm{2}}.

The mean and the variance resulting from the sample and calculated with these Eqs. (6) and (7) are uncertain [⁶⁶,⁶⁷].

Response surface model

A response surface model (RSM) is a collection of statistical and mathematical techniques that are useful for developing, improving, and optimizing processes.The choice of RSM for a given computational model depends on the knowledge of the computational model itself [⁶⁸,⁶⁹]. It is used in the development of an adequate functional relationship between a response of interest y, and a number of associated input parameters denoted by (x₁, x₂, . . ., x_k). In general, such a relationship is unknown but can be approximated by a low-degree polynomial model of the form:

(8)\mathit{y}=\mathit{f}\text{'}(\mathit{x})\mathit{\beta }+\mathit{ϵ},

where \mathit{x} = (x₁, x₂, . . ., x_k), f(x) is a vector function of p elements that consists of powers and cross-products of powers of x₁, x₂, . . ., x_k up to a certain degree denoted by d (≥1), b is a vector of p unknown constant coefficients referred to as parameters, and ε is a random experimental error assumed to have a zero mean. This is conditioned on considering the model provides an adequate representation of the response. In this case, the quantity \mathit{f}\text{'}(\mathit{x})\mathit{\beta } represents the mean response, that is, the expected value of y, and is denoted by m(x). Two important models are commonly used in RSM. These are special cases of model in Eq. (8) and include the first-degree model (d = 1)

(9)\mathit{y}={{\displaystyle \mathit{\beta }}}_0+\mathrm{\sum }_{i=\mathrm{1}}^k{{\displaystyle \mathit{\beta }}}_i{{\displaystyle \mathit{x}}}_i+\mathit{ϵ},

and the second-degree model (d = 2)

(10)

where x₁, x₂, . . . , x_k are the input factors which influence the responsey;β₀, β_ii (i = 1,2,…,k), β_ij (i = 1,2, . . . ,k ; j = 1,2,. . . ,k) are unknown parameters and \mathit{ϵ} is a random error. The β coefficients are obtained by the least squares method. A series of n experiments should first be carried out, in each of which the response y is measured (or observed) for specified settings of the control parameters. The totality of these settings constitutes the so-called response surface design, or just design, which can be represented by a matrix, denoted by D, of order n× k called the design matrix,

(11)D=\left(\begin{array}{l}{{\displaystyle x}}_{\mathrm{11}}{{\displaystyle x}}_{\mathrm{12}}\mathrm{...}{{\displaystyle x}}_{\mathrm{1}k}\\ {{\displaystyle x}}_{\mathrm{21}}{{\displaystyle x}}_{\mathrm{22}}\mathrm{...}{{\displaystyle x}}_{\mathrm{2}k}\\ \mathrm{...}\\ {{\displaystyle x}}_{n\mathrm{1}}{{\displaystyle x}}_{n\mathrm{2}}\mathrm{...}{{\displaystyle x}}_{nk}\end{array}\right),

where {{\displaystyle x}}_{ui} denotes the uth design setting of {{\displaystyle x}}_i (i = 1,2 . . .,k; u = 1,2,. . .,n). Each row of D represents a point, referred to as a design point, in a k-dimensional Euclidean space. Let {{\displaystyle y}}_u denote the response value obtained as a result of applying the u th setting of x, namely {{\displaystyle x}}_u = ({{\displaystyle x}}_{u\mathrm{1}},{{\displaystyle x}}_{u\mathrm{2}},\mathrm{...},{{\displaystyle x}}_{uk}), (u = 1, 2, . . ., n). From Eq. (8), we then have

(12){{\displaystyle \mathit{y}}}_u=\mathit{f}\text{'}({{\displaystyle \mathit{x}}}_u)\mathit{\beta }+{\mathit{ϵ}}_u,u=\mathrm{1},\mathrm{2},\mathrm{...},n,

where {\mathit{ϵ}}_u denotes the error term at the uth experimental run. Equation (12) can be expressed in matrix form as:

(13)\mathit{y}=\mathit{X}\mathit{\beta }+\mathit{ϵ},

where \mathit{y} = ({{\displaystyle y}}_{\mathrm{1}},{{\displaystyle y}}_{\mathrm{2}},\mathrm{...},{{\displaystyle y}}_n), \mathit{X} is a matrix of order n× p whose uth row is \mathit{f}\text{'}({{\displaystyle \mathit{x}}}_u), and \mathit{ϵ}=({ϵ}_{\mathrm{1}},{ϵ}_{\mathrm{2}},\mathrm{...},{ϵ}_n). Note that the first column of X is the column of ones 1_n.

Assuming that \mathit{ϵ} has a zero mean, the so called ordinary least-squares estimator of b is [⁷⁰]:

(14)\widehat{\mathit{\beta }}={(\mathit{X}\text{'}\mathit{X})}^{-\mathrm{1}}\mathit{X}\text{'}\mathit{y}.

Range of parameters

The parameters under study to construct the surrogate model for the damage detection are illustrated in Table 1 where the name and symbol of each parameter would be used in the optimization process. It is worthy to mention that quadratic and interaction terms would be used for the surrogate model.

Results of response surface model

To efficiently construct the surrogate model for the maximum principal stress, we are comparing the actual maximum principal stress from the numerical simulations for the 40 samples of the ground motion simulation with the predicted maximum principal stress quantified from the surrogate model and determine the coefficient of regression so that to construct an efficient representation of the structural system responses as much as we could. The coefficient of regression was R² = 0.9736 (see Fig. 5) which is a very good representation that can be support on to predict the responses of the tunnel during the 10 s of ground motion.

The surrogate model equation for the maximaum principal stress at the top point inside the tunnel lining has been quantified as follows:

(15)\begin{array}{l}Y=-\mathrm{1220729}+\mathrm{538971}{{\displaystyle X}}_{\mathrm{1}}-\mathrm{125488}{{\displaystyle X}}_{\mathrm{2}}+\mathrm{0.01968}{{\displaystyle X}}_{\mathrm{3}}+\mathrm{1410.92}{{\displaystyle X}}_{\mathrm{4}}+\mathrm{796483}{{\displaystyle X}}_{\mathrm{1}}^{\mathrm{2}}\\ +\mathrm{7603.66}{{\displaystyle X}}_{\mathrm{2}}^{\mathrm{2}}-\mathrm{0.00000000025}{{\displaystyle X}}_{\mathrm{3}}^{\mathrm{2}}-\mathrm{0.19799}{{\displaystyle X}}_{\mathrm{4}}^{\mathrm{2}}+\mathrm{16126.3}{{\displaystyle X}}_{\mathrm{1}}{{\displaystyle X}}_{\mathrm{2}}\\ -\mathrm{0.00323}{{\displaystyle X}}_{\mathrm{1}}{{\displaystyle X}}_{\mathrm{3}}-\mathrm{576.641}{{\displaystyle X}}_{\mathrm{1}}{{\displaystyle X}}_{\mathrm{4}}+\mathrm{0.00143}{{\displaystyle X}}_{\mathrm{2}}{{\displaystyle X}}_{\mathrm{3}}-\mathrm{15.5855}{{\displaystyle X}}_{\mathrm{2}}{{\displaystyle X}}_{\mathrm{4}}\\ -\mathrm{0.000004645}{{\displaystyle X}}_{\mathrm{3}}{{\displaystyle X}}_{\mathrm{4}}.\end{array}

BCQO method

The purpose of optimization is to find the best solution or optimal solution from a set of solutions called feasible region. Mathematically, most of the interesting optimization problems can be formulated to optimize some objective function, subject to some equality and/or inequality constraints. The general form of optimization problems is: Min. or Max. f(x₁,x_2,…,x_n) (objective function) Subject to g_i(x) {≤, = ,≥} b_i(i = 1,2,…,m) which are the constraints and l_j≤x_j≤u_j (j= 1,2,…,n) are the bounded constrains.

The goal is to find x₁,x₂,…,x_n that satisfy the constraints and achieve min(max) objective function value.

Regarding the objective function, the classification based on the number of objective functions, objective functions can be classified as single and multiobjective programming problems, linear or nonlinear and the structure of the formulation of the variables.

The formulation problem leads to a quadratic optimization problem which is called Bound-Constrained Quadratic Optimization. It is an objective function subject to bound constraints on the values of the variables. We have a special case of the bound constrained which we have only simple constrained lb\le x\le ub for each variable.

Bound constrained optimization problems has important role in engineering and industries applications. Also they arise where the parameters that describe physical quantities are constrained to be in a given range. The mathematical form is

(16)\mathrm{M}\mathrm{i}\mathrm{n}\hspace{0.17em}\mathit{Z}=\mathrm{0.5}{\mathit{X}}^{\mathrm{T}}\mathit{Q}\mathit{X}+{\mathit{C}}^{\mathrm{T}}\mathit{X},

where \mathit{Q} ∈ R^{n× n} is symmetric and positive definite matrix [⁷¹].

Solving this type of the problem depends on type of the constraints. A solution was introduced by Ref. [⁷²]. through extending the simplex method to solve constrained minimization problem where the constraints were in simple and inequalities form.

In our objective function, X is a vector of variables (x₁, x_2,x₃, x₄)^T, Q is a (4×4) not positive definite matrix and C^T₌ (c₁, c₂, c_3,c₄) [⁷³].

A set of lower and upper bounds on the design variables in x are defined and Matlab program is used to calculate the maximum value of the parameters with regard to the tensile damage in the tunnel lining during the earthquake.

Optimum results

Optimization completed because the objective function is non-decreasing in feasible directions within the default value of the function tolerance, and constraints are satisfied within the default value of the constraint tolerance. The maximum values of the parameters are as follows {{\displaystyle X}}_{\mathrm{1}}=\mathrm{0.4}\mathrm{m},{{\displaystyle X}}_{\mathrm{2}}=\mathrm{10}\mathrm{m},{{\displaystyle X}}_{\mathrm{3}}=\mathrm{20.2}\mathrm{G}\mathrm{P}\mathrm{a},\mathrm{a}\mathrm{n}\mathrm{d}\hspace{0.17em}{{\displaystyle X}}_{\mathrm{4}}=\mathrm{2458.6}\mathrm{k}\mathrm{g}/{\mathrm{m}}^3.

The maximum principal stress in the lining= 5.3894e+ 05 Pa. This value would be compared to the maximum principal stress calculated using XFEM approach for the maximum case to double check the accuracy of the coupling method.

Optimization regarding tunnel damages

The optimization has been studied regarding the responses of the circular tunnel during an earthquake. The parameters are used in three situations as shown in Table 2 which are the maximum, minimum and medium cases by taking the stiffness of the tunnel lining as a criteria, for example for the maximum case the highest value of the concrete density was chosen, the highest lining thickness of the lining was considered, lowest value of the lining diameter was selected and the lowest value of the modulus of elasticity of the concrete was adopted. While for the minimum case, the same mentioned parameters values are the lowest, lowest, highest and highest simultaneously. While for the medium case, the medium values of the parameters were used. Two outputs are considered for the damage responses of the circular tunnel which are compressive damage and tensile damage of the lining.

Compressive damage–design cases

The minimum case for the design (worst case) and when the peak ground acceleration is 4.905 m/s², where the values of the four parameters for the tunnel lining are (density= 2300 kg/m³, thickness= 0.3 m, diameter= 10 m, and Young’s modulus of concrete= 31 GPa) . The results of the simulation at the last 10 s of the earthquake loading are as follows:

The maximum compressive damage in the tunnel lining is 1.44e–02 (see Fig. 6). The damage appears in the right and left sides inside the tunnel lining.

While for the medium case where the values of the four parameters are (density= 2400 kg/m³, thickness= 0.4 m, diameter= 9 m, and Young’s modulus of concrete= 24 GPa), the maximum compressive damage appears in the same region but it decreases to the value 8.015e–03 (see Fig. 7).

But when we consider the maximum design case with the four parameters values (density= 2500 kg/m³, thickness= 0.5 m, diameter= 8 m, and Young’s modulus of concrete= 17 GPa), the maximum compressive damage in the tunnel lining reaches 2.284e–03 (see Fig. 8) which is seen in the left and right sides with a small portion in the top part of the tunnel lining, which means that it has been decreased.

The results of the compressive damage show a semi-linear behavior of the tunnel lining against the earthquake when changing the four parameters, as a result we recognize that the maximum case is the better design case for the tunnel structure to withstand the compressive stresses in the tunnel lining because the maximum compressive damage in this case is 2.284e–03 which is smaller than its values both in minimum and medium cases which are 1.44e–02 and 8.015e–03, respectively.

Tensile damage – design cases

In the same way, depending on the results of the optimization process, the minimum case for the design where the values of the four parameters for the tunnel lining are (density= 2300 kg/m³, thickness= 0.3 m, diameter= 10 m, and Young’s modulus of concrete= 31 GPa) at the last 10 s of the earthquake loading the maximum tensile damage in the tunnel lining is 6.501e–01 (see Fig. 9). The damage appears in the top and bottom sides inside the tunnel lining.

While for the medium case where the values of the four parameters are (density= 2400 kg/m³, thickness= 0.4 m, diameter= 9 m, and Young’s modulus of concrete= 24 GPa), the maximum tensile damage appears in the same region but it increases to the value 7.883e–01 (see Fig. 10).

But when we consider the maximum design case with the four parameters values (density= 2500 kg/m³, thickness= 0.5 m, diameter= 8 m, and Young’s modulus of concrete= 17 GPa), the maximum tensile damage in the tunnel lining reaches 7.479e–01 (see Fig. 11) which is seen in the top and bottom sides of the tunnel lining, which means that it has been decreased.

The results of the tensile damage show a nonlinear behavior of the tunnel lining against the earthquake when changing the four parameters, we conclude that the minimum case is the better design case for the tunnel structure to withstand the tensile stresses in the tunnel lining because the maximum tensile damage in this case is 6.51e–01 which is smaller than its values both in medium and maximum cases which are 7.883e–01 and 7.479e–01, respectively.

Application of XFEM for tensile damage

We are dedicating the extended finite element XFEM approach in detecting the cracks in addition to RSM approach for doubling the accuracy of the structural design optimization of the tunnel lining. In the following section we will discuss the results in details.

Minimum case

When the minimum case is considered for the tensile damage, and after applying the XFEM approach we see the start of the crack at the top part inside the tunnel lining a little to the left side as seen in Fig. 12 starting at 3.607 s after earthquake loading.The crack continues to propagate starting from 6.138 s toward the top of the tunnel lining and the crack stops without further propagation, see Fig. 13. The maximum principal stress in the minimum case is 3.234e+ 06 Pa where the output is very crucial in detecting the propafation of cracks in the tunnel lining. The principal stress depends on the normal stresses which are normal to the cross-section at each point of consideration.

Medium case

While for the medium case, we see the start of the crack at the top part inside the tunnel lining as seen in Fig. 14 starting at 2.454 s. The crack continues to propagate starting from 5.565 s toward the top of the tunnel lining and the crack stops from propagation, see Fig. 15. It is worthy to mention that the maximum principal stress in this case is 6.253e+ 06 Pa. The magnitude of this output has been increased to double of its value approximately compared to the previous minimum case, which is an indication that the cracks appears earlier in the medium case regardless of increasing the thickness, decreasing the diameter, decreasing the Young’s modulus and increasing the density of the concrete for the tunnel lining. The head of the overburden of the soil on the top part of the tunnel becomes higher with 1 m height which generates further static and impact stresses in the tunnel lining at this region which is beyond the generation of cracks earlier in this case.

Maximum case

But for the maximum case, the crack starts to appear at the top part inside the tunnel lining a little to the right side as seen in Fig. 16 starting at 1.981 s. The crack continues to propagate starting from 4.937 s toward the left side of the tunnel lining and the crack propagates further in the same direction at 4.941 s and finally stops from propagation, see Figs. 17 and 18 and it is worthy to mention that the damage in this case is near to spalling than the normal cracks compared to the previous two cases. For this case, the maximum principal stress is 3.314e+ 06 Pa. The magnitude of this output has been decreased to half of its value approximately compared to the medium case, also the cracks appears earlier in the maximum case regardless of increasing the thickness, decreasing the diameter, decreasing the Young’s modulus and increasing the density of the concrete for the tunnel lining. The head of the overburden of the soil on the top part of the tunnel becomes higher with 2 m height in this case which generates further static and impact stresses in the tunnel lining at this region, but the decrease in the circumference of the tunnel lining decreases the static and impact stresses which is the main reason of decreasing the maximum principal stress.

The responses of the tunnel lining with respect to the tensile damage and the crack during the earthquake for 10 s showed a nonlinear behavior regarding optimization of the four parameters of the tunnel lining, where the normal stresses are not dependent on stronger concrete only but the overburden height of the soil also plays a vital role in determining the optimum design case of the tunnel lining.

The maximum principal stress is 3.314e+ 06 Pa compared to the same output value 5.3894e+ 05 Pa calculated previously in the BCQO method are similar to the extent of 83.74%, which means that we can use this coupling method for the purpose of structural optimization for the tensile damage in the tunnel lining due to ground motion.

Conclusions

The results of the numerical simulations using XFEM-BCQO detected many significant points for the structural design of tunnel lining as follows:

1) The Monte Carlo sampling method showed efficient approach in detecting the responses of the tunnel lining during earthquakes supporting on the numerical simulations through quantifying optimum design case considering the surrogate model.

2) The XFEM approach manifested a great role in predicting the time and positions of the tensile cracks in the tunnel lining by standing on the results of the surrogate model from Monte Carlo sampling mehod in addition to using BCQO method to verify the results of the structural optimization of the four parameters which showed significant efficiency.

3) The concrete tensile strength of the tunnel lining isn’t the only controlling factor to avoid cracks and failure during earthquakes, but the surrounding soil type properties, the head of the overburden weight of the soil above the tunnel lining and the external circumferencial area of the tunnel are effective factors, too.

4) The most sensitive part of the tunnel lining is the top region of it where the combination of the static and the dynamic impact load of the overburden soil above the tunnel during the earthquake is critical in generating the early cracks and spalling as seen in the maximum design case of the tunnel lining.

References

[1] Pescara M, Gaspari G M, Repetto L. Design of underground structures under seismic conditions: A long deep tunnel and a metro tunnel. In: Colloquium on Seismic Design of Tunnels. Torino: Geodata Engineering SpA, 2011

[2] Hashash Y M A, Hook J J, Schmidt B, I-Chiang Yao J. Seismic design and analysis of underground structure. Journal of Tunneling and Underground Space Technology, 2001, 16(4): 247–293

[3] Hashash Y M A, Park D, Yao J I. Ovaling deformations of circular tunnels under seismic loading: An update on seismic design and analysis of underground structures. Journal of Tunneling and Underground Space Technology, 2005, 20(5): 435–441

[4] St John C M, Zahrah T F. Aseismic design of underground structures. Tunnelling and Underground Space Technology, 1987, 2(2): 165–197

[5] Kawashima K. Seismic design of underground structures in soft ground: A review. In: Kusakabe, Fujita, Miyazaki, eds. Geotechnical Aspects of Underground Construction in Soft Ground. Rotterdam, 1999

[6] Fabozzi S. Behaviour of segmental tunnel lining under static and dynamic loads. Dissertation for the Doctoral Degree. Naples: University of Naples Federico II, 2017

[7] Nariman N A, Hussain R R, Msekh M A, Karampour P. Prediction meta-models for the responses of a circular tunnel during earthquakes. Underground Space, 2019, 4(1): 31–47

[8] Moller S C, Vermeer P A. On numerical simulation of tunnel installation. Tunnelling and Underground Space Technology (Oxford, England), 2008, 23(4): 461–475

[9] Lekhnitskii S G. Anisotropic plates. London: Foreign Technology Div Wright-Patterson Afb Oh, 1968

[10] Lu A Z, Zhang L Q, Zhang N. Analytic stress solutions for a circular pressure tunnel at pressure and great depth including support delay. International Journal of Rock Mechanics and Mining Sciences, 2011, 48(3): 514–519

[11] Yasuda N, Tsukada K, Asakura T. Elastic solutions for circular tunnel with void behind lining. Tunnelling and Underground Space Technology (Oxford, England), 2017, 70: 274–285

[12] Han X, Xia Y. Analytic solutions of the forces and displacements for multicentre circular arc tunnels. Hindawi Mathematical Problems in Engineering, 2018, 2018: 8409129

[13] Schmid H. Static problems of tunnels and pressure tunnels construction and their mutual relationships. Berlin: Springer, 1926

[14] Morgan H. A contribution to the analysis of stress in a circular tunnel. Geotechnique, 1961, 11(1): 37–46

[15] Windels R. Kreisring im elastischen continuum. Bauingenieur, 1967, 42: 429–439

[16] Duddeck H, Erdmann J. On structural design models for tunnels in soft soil. Underground Space (United States), 1985, 9(5–6): 246–259

[17] Do N A, Dias D, Oreste P, Djeran-Maigre I. A new numerical approach to the hyperstatic reaction method for segmental tunnel linings. International Journal for Numerical and Analytical Methods in Geomechanics, 2014, 38(15): 1617–1632

[18] Vu Minh N, Broere W, Bosch J W. Structural analysis for shallow tunnels in soft soils. International Journal of Geomechanics, 2017, 17(8): 04017038

[19] Ghasemi H, Park H S, Rabczuk T. A level-set based IGA formulation for topology optimization of flexoelectric materials. Computer Methods in Applied Mechanics and Engineering, 2017, 313: 239–258

[20] Ghasemi H, Brighenti R, Zhuang X, Muthu J, Rabczuk T. Optimum fiber content and distribution in fiber-reinforced solids using a reliability and NURBS based sequential optimization approach. Structural and Multidisciplinary Optimization, 2015, 51(1): 99–112

[21] Zhang C, Nanthakumar S S, Lahmer T, Rabczuk T. Multiple cracks identification for piezoelectric structures. International Journal of Fracture, 2017, 206(2): 151–169

[22] Nanthakumar S, Zhuang X, Park H, Rabczuk T. Topology optimization of flexoelectric structures. Journal of the Mechanics and Physics of Solids, 2017, 105: 217–234

[23] Nanthakumar S, Lahmer T, Zhuang X, Park H S, Rabczuk T. Topology optimization of piezoelectric nanostructures. Journal of the Mechanics and Physics of Solids, 2016, 94: 316–335

[24] Nanthakumar S, Lahmer T, Zhuang X, Zi G, Rabczuk T. Detection of material interfaces using a regularized level set method in piezoelectric structures. Inverse Problems in Science and Engineering, 2016, 24(1): 153–176

[25] Nanthakumar S, Valizadeh N, Park H, Rabczuk T. Surface effects on shape and topology optimization of nanostructures. Computational Mechanics, 2015, 56(1): 97–112

[26] Nanthakumar S S, Lahmer T, Rabczuk T. Detection of flaws in piezoelectric structures using extended FEM. International Journal for Numerical Methods in Engineering, 2013, 96(6): 373–389

[27] Vu-Bac N, Lahmer T, Zhuang X, Nguyen-Thoi T, Rabczuk T. A software framework for probabilistic sensitivity analysis for computationally expensive models. Advances in Engineering Software, 2016, 100: 19–31

[28] Vu-Bac N, Silani M, Lahmer T, Zhuang X, Rabczuk T. A unified framework for stochastic predictions of mechanical properties of polymeric nanocomposites. Computational Materials Science, 2015, 96: 520–535

[29] Vu-Bac N, Rafiee R, Zhuang X, Lahmer T, Rabczuk T. Uncertainty quantification for multiscale modeling of polymer nanocomposites with correlated parameters. Composites. Part B, Engineering, 2015, 68: 446–464

[30] Rabczuk T, Akkermann J, Eibl J. A numerical model for reinforced concrete structures. International Journal of Solids and Structures, 2005, 42(5–6): 1327–1354

[31] Bažant Z P. Why continuum damage is nonlocal: Micromechanics arguments. Journal of Engineering Mechanics, 1991, 117(5): 1070–1087

[32] Thai T Q, Rabczuk T, Bazilevs Y, Meschke G. A higher-order stress-based gradient-enhanced damage model based on isogeometric analysis. Computer Methods in Applied Mechanics and Engineering, 2016, 304: 584–604

[33] Fleck N A, Hutchinson J W. A phenomenological theory for strain gradient effects in plasticity. Journal of the Mechanics and Physics of Solids, 1993, 41(12): 1825–1857

[34] Rabczuk T, Eibl J. Simulation of high velocity concrete fragmentation using SPH/MLSPH. International Journal for Numerical Methods in Engineering, 2003, 56(10): 1421–1444

[35] Rabczuk T, Eibl J, Stempniewski L. Numerical analysis of high speed concrete fragmentation using a meshfree Lagrangian method. Engineering Fracture Mechanics, 2004, 71(4–6): 547–556

[36] Rabczuk T, Xiao S P, Sauer M. Coupling of meshfree methods with nite elements: Basic concepts and test results. Communications in Numerical Methods in Engineering, 2006, 22(10): 1031–1065

[37] Rabczuk T, Eibl J. Modelling dynamic failure of concrete with meshfree methods. International Journal of Impact Engineering, 2006, 32(11): 1878–1897

[38] Etse G, Willam K. Failure analysis of elastoviscoplastic material models. Journal of Engineering Mechanics, 1999, 125(1): 60–69

[39] Miehe C, Hofacker M, Welschinger F. A phase field model for rate-independent crack propagation: Robust algorithmic implementation based on operator splits. Computer Methods in Applied Mechanics and Engineering, 2010, 199(45–48): 2765–2778

[40] Amiri F, Millan D, Arroyo M, Silani M, Rabczuk T. Fourth order phase-field model for local max-ent approximants applied to crack propagation. Computer Methods in Applied Mechanics and Engineering, 2016, 312(C): 254–275

[41] Areias P, Rabczuk T, Msekh M. Phase-field analysis of finite-strain plates and shells including element subdivision. Computer Methods in Applied Mechanics and Engineering, 2016, 312(C): 322–350

[42] Msekh M A, Silani M, Jamshidian M, Areias P, Zhuang X, Zi G, He P, Rabczuk T. Predictions of J integral and tensile strength of clay/epoxy nanocomposites material using phase-eld model. Composites. Part B, Engineering, 2016, 93: 97–114

[43] Hamdia K, Msekh M A, Silani M, Vu-Bac N, Zhuang X, Nguyen-Thoi T, Rabczuk T. Uncertainty quantication of the fracture properties of polymeric nanocomposites based on phase field modeling. Composite Structures, 2015, 133: 1177–1190

[44] Msekh M A, Sargado M, Jamshidian M, Areias P, Rabczuk T. ABAQUS implementation of phase-field model for brittle fracture. Computational Materials Science, 2015, 96: 472–484

[45] Amiri F, Millán D, Shen Y, Rabczuk T, Arroyo M. Phase-field modeling of fracture in linear thin shells. Theoretical and Applied Fracture Mechanics, 2014, 69: 102–109

[46] Hamdia K M, Zhuang X, He P, Rabczuk T. Fracture toughness of polymeric particle nanocomposites: Evaluation of Models performance using Bayesian method. Composites Science and Technology, 2016, 126: 122–129

[47] Rabczuk T, Belytschko T, Xiao S P. Stable particle methods based on Lagrangian kernels. Computer Methods in Applied Mechanics and Engineering, 2004, 193(12–14): 1035–1063

[48] Rabczuk T, Belytschko T. Adaptivity for structured meshfree particle methods in 2D and 3D. International Journal for Numerical Methods in Engineering, 2005, 63(11): 1559–1582

[49] Nguyen V P, Rabczuk T, Bordas S, Duflot M. Meshless methods: A review and computer implementation aspects. Mathematics and Computers in Simulation, 2008, 79(3): 763–813

[50] Zhuang X, Cai Y, Augarde C. A meshless sub-region radial point interpolation method for accurate calculation of crack tip elds. Theoretical and Applied Fracture Mechanics, 2014, 69: 118–125

[51] Zhuang X, Zhu H, Augarde C. An improved meshless Shepard and least squares method possessing the delta property and requiring no singular weight function. Computational Mechanics, 2014, 53(2): 343–357

[52] Zhuang X, Augarde C, Mathisen K. Fracture modelling using meshless methods and level sets in 3D: Framework and modelling. International Journal for Numerical Methods in Engineering, 2012, 92(11): 969–998

[53] Chen L, Rabczuk T, Bordas S, Liu G R, Zeng K Y, Kerfriden P.Extended finite element method with edge-based strain smoothing (Esm-XFEM) for linear elastic crack growth. Computer Methods in Applied Mechanics and Engineering, 2012, 209–212(4): 250–265

[54] Belytschko T, Black T. Elastic crack growth in finite elements with minimal remeshing. International Journal for Numerical Methods in Engineering, 1999, 45(5): 601–620

[55] Moës N, Dolbow J, Belytschko T. A finite element method for crack growth without remeshing. International Journal for Numerical Methods in Engineering, 1999, 46(1): 131–150

[56] Vu-Bac N, Nguyen-Xuan H, Chen L, Lee C K, Zi G, Zhuang X, Liu G R, Rabczuk T. A phantom-node method with edge-based strain smoothing for linear elastic fracture mechanics. Journal of Applied Mathematics, 2013, 2013: 978026

[57] Bordas S P A, Natarajan S, Kerfriden P, Augarde C E, Mahapatra D R, Rabczuk T, Pont S D. On the performance of strain smoothing for quadratic and enriched finite element approximations (XFEM/GFEM/PUFEM). International Journal for Numerical Methods in Engineering, 2011, 86(4–5): 637–666

[58] Bordas S P A, Rabczuk T, Hung N X, Nguyen V P, Natarajan S, Bog T, Quan D M, Hiep N V. Strain smoothing in FEM and XFEM. Computers & Structures, 2010, 88(23–24): 1419–1443

[59] Rabczuk T, Zi G, Gerstenberger A, Wall W A. A new crack tip element for the phantom node method with arbitrary cohesive cracks. International Journal for Numerical Methods in Engineering, 2008, 75(5): 577–599

[60] Chau-Dinh T, Zi G, Lee P S, Rabczuk T, Song J H. Phantom-node method for shell models with arbitrary cracks. Computers & Structures, 2012, 92–93: 242–256

[61] Song J H, Areias P M A, Belytschko T. A method for dynamic crack and shear band propagation with phantom nodes. International Journal for Numerical Methods in Engineering, 2006, 67(6): 868–893

[62] Areias P M A, Song J H, Belytschko T. Analysis of fracture in thin shells by overlapping paired elements. Computer Methods in Applied Mechanics and Engineering, 2006, 195(41–43): 5343–5360

[63] Zamani R, Motahari M R. The effect of soil stiffness variations on Tunnel Lining Internal Forces under seismic loading and Case comparison with existing analytical methods. Ciência e Natura, Santa Maria, 2015, 37(1): 476–487

[64] Möller S C. Tunnel induced settlements and structural forces in linings. Dissertation for the Doctoral Degree. Stuttgart: University of Stuttgart, 2006

[65] Lu Q, Chen S, Chan Y, He C. Comparison between numerical and analytical analysis of the dynamic behavior of circular tunnels. Earth Sciences Research Journal, 2018, 22(2): 119–128

[66] Saltelli A, Ratto M, Andres T, Campolongo F, Cariboni J, Gatelli D. Global Sensitivity Analysis: The Primer. Hoboken: John Wiley & Sons Ltd., 2008

[67] Burhenne S, Jacob D, Henze G P. Sampling based on Sobol sequences for monte carlo techniques applied to building simulations. In: The 12th Conference of International Building Performance Simulation Association. Sydney, 2011

[68] Myers R H, Montgomery D C. Response Surface Methodology: Product and Process Op-timization Using Designed Experiments. 2nd ed. New York: John Wiley & Sons, 2002

[69] Zhao J, Tiede C. Using a variance-based sensitivity analysis for analyzing the relation between measurements and unknown parameters of a physical model. Nonlinear Processes in Geophysics, 2011, 18(3): 269–276

[70] Khuri1 A I, Mukhopadhyay S. Response surface methodology. WIREs Computational Statistics, 2010, 2(2): 128–149

[71] Luenberger D G, Ye Y. Linear and Non-linear programming. In: International Series in Operations Research & Management Science. Palo Alto, CA: Stanford University, 2015

[72] Box M J. A new method of constrained optimization and a comparison with other methods. Computer Journal, 1965, 8(1): 42–52

[73] Hunt B R, Lipsman R L, Rosenberg J M. A Guide to MATLAB for Beginner and Experienced Users. Cambridge: Cambridge University Press, 2006

RIGHTS & PERMISSIONS

Higher Education Press and Springer-Verlag GmbH Germany, part of Springer Nature