#### Abstract

TBM cutterhead driving system is generally an extraordinarily large and complex machine containing lots of parameters; understanding and assessment of its dynamic characteristics are a great challenge as each of these parameters has uncertainty. In this work, a hierarchical modeling method for the dynamic model of the complex gear transmission system is proposed based on the generalized finite element modeling idea. On this basis, the whole machine dynamic model of cutterhead driving system is established; both the characteristics of vibration responses and the meshing force are revealed. Vibration responses under the action of simulated load are estimated and verified by comparing with the data measured on the tunneling field, where the error is about 18%~50%. With the vibration response of the key nodes and the dynamic meshing force for the system dynamic characteristic evaluation index, considering the change of input parameters such as external loads, material parameters, meshing parameters, and coupling parameters, the global parameter sensitivity of system dynamic characteristics is analyzed based on the technology of surrogate model. Finally, variation of dynamic characteristics considering the interaction of polytypic parameters is obtained.

#### 1. Introduction

A tunnel boring machine (TBM) is a typical complex and major equipment, which has been widely used in excavation of tunnels under the subsoil all over the world. TBM based excavation is a very complex and rigorous multidisciplinary and multiobjective engineering problem including mechanical, civil, control, hydraulic, planning, and time-cost problems, whose security has to be primarily guaranteed as the associated risk can be extremely high [1–4]. Despite the external factors such as the human misoperation, the dynamic performance of TBM as well as the geologic condition can be identified as the main causative factor influencing the safe and efficient excavation of tunnels. Therefore, analyzing the dynamic performance of TBM under complex geologic condition is a major concern.

As TBM is an extraordinarily large and complex machine consisting of many subsystems and lots of components, experimental analysis of a full machine or even one of the main subsystems/components under real condition is generally not practicable, which results in a lack of understanding and assessment of its realistic performance and thus the further design optimization. Recently, due to the advances in computer hardware and numerical methods, computational simulation has been becoming a popular and promising method/tool at the early stage of a TBM design to assess its various performances such as structural strength, vibration, and control characteristics. For instance, Swoboda et al. adopted Finite Element Analysis (FEA) to investigate TBM tunneling in saturated porous medium. They found that the realistic modeling of soil behavior can be assessed [5]. Sun et al. carried out dynamics modeling to estimate the dynamics of a TBM; a multidegree-of-freedom coupling dynamic model was presented to obtain the cutterhead joint surface loads [6, 7]. In addition, much work has been done in an effort to understand the complex behavior of TBM and considerable information of value has been obtained [8–12]. However, the aforementioned research falls into the category of conventional deterministic analysis, in which the effect of statistical variation has not been taken into account. Besides, most of them only investigated the effect of a few of parameters (no more than five) on the behavior of TBM, and the conventional analysis significantly relies on “one-factor-at-a-time” (OFAT) method, in which only one parameter is varied while keeping others at their baseline values so that the global and comparative effect among different parameters remain vague. On the contrary, there are usually many parameters involved in even one subsystem or one component; some of them are associated with the variation due to the varying loads and control parameters, and the others are associated with variation due to the uncertainties of material properties, manufacturing errors, and assembling tolerance. Therefore, it is inadequate to investigate a TBM with a deterministic method, and thus a method to deal with these uncertainties involved and to identify and rank the most important varied and/or uncertain parameters is necessary for better understanding and further optimization of a TBM.

Sensitivity analysis (SA) is a promising approach to achieve it as it is able to identify the effect of input parameter uncertainty/variation on output response [13]. SA can be categorized into the local sensitivity analysis (LSA) and global sensitivity analysis (GSA) [14–16]. The former emphasizes the effect of small parameter variations on model response and can be used to determine model response changes with individual parameter [17–19]. In contrast, GSA can be used to determine the most contributing input parameter to an output behavior as the noninfluential inputs or ascertain some interaction effects within the model over the whole design space [20]. Specifically, GSA of a complex engineering product/system is able to tell (1) which parameters contribute the most to the system response and particular attention should be given in the next design optimization, (2) which parameters interact with each other and a potential cause could be identified, (3) which parameters are insignificant and can be eliminated to reduce the dimensionality of the design. A number of GSA methods have been developed, like Morris method [21], Iterated Fractional Factorial Design [22], Fourier Amplitude Sensitivity Test (FAST) [23], Sobol’s method [24], and so on. Among these GSA methods, Sobol’s method has been widely used recently in a variety of fields such as geologic, chemical, and mechanical engineering [17, 25–27]. For example, Kala et al. employed Sobol’s method to investigate stability problems of steel plane frames [25]. Queipo et al. used Sobol’s method to rank the input parameters including flow angle, hydrogen area, oxygen area, and oxidizer post tip thickness based on their contribution to the improvement of performance and life of liquid-rocket injector [17]. Song et al. applied Sobol’s method to study the effect of seven parameters on the crashworthiness of thin-walled structures [26]. Fu et al. utilized Sobol’s method to investigate the impact of twenty-one parameters and selected most important three as the design variables for the subsequent optimization of water distribution system [27]. Availability and reliability of Sobol’s GSA have been verified through these researches. Nevertheless, the efficiency of Sobol’s method has also been found inadequate as Sobol’s method necessitates a large number of calculations, which becomes a major barrier for Sobol’s GSA to apply in computationally expensive simulation such as the dynamics modeling of TBM.

There are generally three solutions to overcome this difficulty: (1) to simplify the computationally expensive model or to use low-fidelity model instead of high-fidelity model to reduce the calculation time of each run [28], (2) to use a distributed parallel computing algorithm and/or workstation to accelerate the speed of the computation [29, 30], and (3) to employ surrogate modeling technique to construct a computationally cheap approximation model instead of the original time consuming model [17, 31]. The first one sometimes is not applicable as it is difficult to make a tradeoff between the efficiency and accuracy. In most cases, less simplification will not be able to adequately relieve the computational burden and much more simplification will lead to loss of important characteristics of the model. The second one requires an appropriate computing architecture, sufficient computational resource, and much work on coding and debugging, which makes it unaffordable to researchers from the engineering design community [29, 30]. In this context, surrogate modeling which could considerately alleviate the computational burden with less modeling error becomes a good alternative to be combined with Sobol’s GSA method [17].

Surrogate modeling is an approach to approximate the implicit relationship between input and output with a limited number of computational simulations, namely, sampling data, and thus efficient to infer output values at yet-unseen input points. To date, there has been a wide variety of surrogate models developed and applied in various engineering problems [32–35]. Huang et al. developed a redial basis function model to determine the optimal loading parameters in T-shape tube hydroforming process [32]. Xiang et al. presented a study on applying Support Vector Regression in the design optimization of the protective effect of railway wind barriers [33]. Fang et al. proposed RBF model and compared it with other in the application of multiobjective crashworthiness optimization [34]. Surrogate models such as Kriging model (KRG), Radial Basis Function (RBF) model, and support machine regression (SVR) have been proven reliable and efficient in these examples. However, it also has been clarified that these different surrogate models often provide rather different modeling accuracy, gradient information, and design outcomes. It is yet unclear which surrogate model is most appropriate for a particular engineering case. Researchers have to select one of them in terms of their experience or literatures in their application. The present study is focusing on three typical surrogate models, namely, RBF, KRG, and SVR, for constructing the mapping between the input and output. One of these three surrogate models will be used for the GSA as long as its errors or metrics meet the accuracy requirements.

Therefore, the aim of this paper is to develop a surrogate model-based GSA methodology for uncertainty analysis of TBM, which comprises surrogate model, Sobol’s GSA, and multilevel dynamics model of TBM. It demonstrates how such a large and complicated model can be investigated through the proposed process. The remainder of the paper is organized as follows. Section 2 gives a brief introduction of TBM and describes multilevel modeling of a practical TBM, Section 3 describes the GSA process based on surrogate model, Section 4 depicts the GSA results of the TBM, and Section 5 finally draws some conclusions.

#### 2. Dynamic Modeling and Vibration Prediction of TBM Cutterhead Driving System

##### 2.1. Dynamic Modeling

###### 2.1.1. Hierarchical Modeling Method and the Dynamic Model

An open-type TBM used for the current study is shown in Figure 1 (left); it can be decomposed into two systems including the main engine system and the back auxiliary system, which are connected by the bridge shelf. The main engine system mainly consists of some subsystems, like the cutterhead and cutter system, the cutterhead driving system, the supporting and the thrust system, etc. Among these subsystems, the cutterhead driving system (see Figure 1 (right)) can be recognized as the most important one as it not only provides the platform to support the other subsystems, but also undergoes and/or absorbs the impact load arising at the interface between the cutterhead and the front rock face. Therefore, the performance of the whole machine significantly depends on the performance of the cutterhead driving system.

To evaluate the dynamics characteristics of the cutterhead driving system, a hierarchical dynamics model is constructed as shown in Figure 2 (see our previous work [1, 6, 7]). In the model, the cutterhead driving system is divided into three subsystems including the cutterhead-bearing-gear subsystem, the driveshaft subsystem, and the carrier-girder subsystem. The basic components in each subsystem are equivalent to virtual shaft, which are modeled as Timoshenko beam elements. The gear meshing element is modeled in the space coordinate, considering the whole freedoms. The bearing/bolt/spline connections are equivalent to the spring-damper element for supporting or connection. The interaction effect between the cutterhead and the rock is equivalent to the coupling element acting on one supporting node.

Therefore, the dynamic equation of the driving system can be defined aswhere is the mass matrix; is the damping matrix; is the stiffness matrix; is the force vector including the external additional load; is the displacement vector of the all nodes.

###### 2.1.2. Definition of Parameters in the Dynamic Model

The definition of parameters including time-varying meshing stiffness, error, backlash, and the supporting stiffness and damping of bearing that has been introduced in [1]. In this work, the gear meshing stiffness is calculated using Ishikawa formula, and the tangential composite tolerance is considered as the gear meshing error. The mesh phases of the multipinions meshing system and the planetary gear system are considered in both the time-varying gear meshing stiffness and error, which are shown in Figure 3. The gear backlash leads to significant nonlinearity into the dynamic model. The damping of gear mesh is calculated by the gear mass and meshing stiffness according to the damping ratio, which is between 0.03 and 0.17.

###### 2.1.3. Numerical Method of the Nonlinear Dynamic Equation

In this paper, the nonlinear dynamic equation (1) is solved by using the precise time-integration method [3], which adopts Hamilton method to translate the ordinary equation to the state space equation aswhere is a -dimensional vector; is the -dimensional constant matrix which can be calculated by , , and in (1); is the column vector constituted by in (1). The common solution of (2) iswhere in the first term can be calculated by using the approximate value of exponential matrix; the integral term is calculated by using the Runge-Kutta method.

##### 2.2. Load Transmission Characteristics and Vibration Prediction

###### 2.2.1. Load Transmission Characteristics

The dynamic rock breaking loads are simulated by the three-dimensional rock breaking model established in LS-DYNA with three types of hob groups including the central cutterhead, the direct cutterhead, and the side cutterhead. The cutterhead loads are integrated by each type of cutter load according to the cutter distributions. The cutterhead of DZ101 TBM which is used as an example in this paper is shown in Figure 4(left), where the cutterhead diameter is 7.93m, and the quantities of three types of cutter are 4, 32, and 11, respectively. Figure 4 (right) shows the time-domain curves and frequency spectrum of the simulated TBM cutterhead loads. Table 1 lists the mean values and main frequency components of the forces and torques at the cutter head in three directions, which have been validated by the field test data recorded in the water diversion project of Jilin province from June to July in 2015.

###### 2.2.2. Vibration Analysis

By using the precise time-integration method of nonlinear dynamic equation, the vibration responses of TBM cutterhead driving system under the simulated loading can be predicted according to the dynamic model established in Section 2.1. Figure 5 shows the time-domain and frequency-domain curves of the cutterhead vibration acceleration and the dynamic pinion-ring gear meshing force. In Figure 5(b), the main vibration frequencies of the axial vibration acceleration on the cutterhead include the axial natural frequency (A1) and the adjacent loading frequency (130.3Hz); Figure 5(d) shows that the main vibration frequencies of the dynamic pinion-ring gear meshing force are the meshing frequency () and its frequency multiplication, as well as the loading frequency (0.89Hz, 130.3Hz).

##### 2.3. Result Comparison of Simulation and Field Test

In order to validate the dynamic model, the vibration results of simulation have been compared with the field data, which is recorded by the wireless acceleration sensors arranged on the surface of the key components such as the cutterhead, shield, walkways, girder, and gripper of the DZ101 TBM. Figure 6 shows the scheme of the detecting system, which is composed of wireless accelerometer, data receiving system, and data processing system.

Figure 7 shows the vibration acceleration of the cutterhead point and the RMS value of each test point. Figure 7(a) shows that the vibration acceleration on the tested point varies in the range of ±4g. In Figure 7(b), the RMS values of vibration acceleration on the cutterhead point in three directions are greater than the values on the girder and gripper successively, which indicates that the vibration decreases progressively along the transfer direction of the power flow.

Table 2 compares the measured and numerical results of the axial vibration acceleration on the three typical test points, where “Node” means the node number marked in the hierarchical dynamic model of TBM cutterhead driving system as shown in Figure 2. “Relative Error” is calculated as where and mean the measured value and predicted value, respectively.

In Table 2, the errors between the calculation value and measured value are among the range of 18%~50%, where the maximum error occurs on the gripper point. However, considering the complexity of both the structure characteristics of TBM cutterhead driving system and the construction environment, the relative error between predicted and measured values is acceptable. That is, the numerical model established in this paper is able to obtain the dynamic characteristics of TBM cutterhead driving system and can be used for the subsequent analysis.

#### 3. Global Sensitivity Analysis Based on Surrogate Models

##### 3.1. Definition of GSA Problems

###### 3.1.1. Outputs for Evaluating the Dynamic Response

Vibration displacements, velocities, and acceleration of the key components are usually used as the measurable indicators of the dynamic characteristics for the oscillatory system. In this paper, the vibration displacement (RMS value) of the meshing node on ring gear (node 8), the vibration acceleration (maximum value) of the meshing node on the pinion (node 14), and the meshing force of the pinion-ring gear mesh are chosen as the indicators for the dynamic response analysis. Both the vibration displacement and acceleration in three directions are considered. And the outputs of dynamic meshing force contain the maximum value and the dynamic deviation factor, which is calculated aswhere means the dynamic deviation factor;* N* is the number of pinions;* T* is the sampling amount; is the transient meshing force of the* j*th pinion at the* i*th sampling point; is the calculated static meshing force. The three types of outputs for evaluating the dynamic responses are listed in Table 3.

###### 3.1.2. Input Parameters and Their Variation Ranges

The uncertain input parameters considered in this paper include the external conditions, the material properties, the meshing parameters, and the connection parameters. The 25 input parameters as well as their variation are listed in Table 4.

*(1) External Conditions*. External conditions contain the rotating speed of cutterhead, the cutterhead torque, and the thrust force. Their initial values and upper and lower bounds are defined according to the statistical data shown in Figure 5. The variation of the portrait and crosswise unbalance force and overturning moment are considered in this paper as well. The ratio between the initial and instant value of torque and thrust force is applied to unbalance forces and moments, respectively.

*(2) Material Properties*. The material properties contain the damping ratio, the elasticity modulus, and the Poisson ratio. The initial values of these parameters and their variations are defined based on experimental test.

*(3) Meshing Parameters*. The meshing parameters contain the backlash, the meshing error, the average meshing stiffness, and meshing damping ratio. Calculation methods of these meshing parameters have been introduced in our previous work. The bounds of gear backlash are calculated according to the center distance and thickness of the engagement pair.

*(4) Connection Parameters*. The connection parameters contain the connecting stiffness and damping ratios of the rock contact, the bolt connection, the bearing supporting, the spline connection, and the hydraulic supporting in the dynamic model. The initial values of these connection parameters have been verified in Section 2.3, which are listed in Table 5. Different fluctuated ratios are considered in this paper, where the variation range of damping ratios is [, 5], and the stiffness is [, 2].

##### 3.2. Design of Experiments (DoE)

DoE are required to generate the initial sampling points before constructing the surrogate models for the sensitivity analysis. In this paper, the Latin Hypercube Sampling (LHS) method is used as the method of DoE due to its good space exploring capacity [36, 37]. 200 sampling points are originally generated according to the bound of 25 input parameters as shown in Table 4.

##### 3.3. Surrogate Models Construction and Accuracy Check

The commonly used surrogate models contain PRS, RBF, KRG, and SVR models [32, 35]. Considering the large number of input parameters and potential high nonlinearity, the surrogate models are constructed using RBF, KRG, and SVR model in this paper [38–41]. Additional 20 sampling points generated by LHS are also used to check the accuracy of the constructed approximation models (Figure 8), where the R-square value is used as the error metric:where* m* is the number of test points; is the response of the test points; is the calculation result of surrogate model; is the average value of responses. According to (6), the closer the value approximates to 1, the better the fitting accuracy is.

#### 4. Results and Discussions

##### 4.1. Construction of Surrogate Models

Figures 9 and 10 show the accuracy comparison results of the three different surrogate models. In Figure 10, s of the vibration displacement (RDef) are the greatest of all, and the smallest is corresponding to the dynamic deviation factor (MF_Se), but it is also greater than 0.8, which means that the accuracies of the constructed surrogate models are all sufficiently high.

To further check the accuracy of the constructed approximation models, a cross validation using 10-fold strategy is conducted [42]. The results are as listed in Table 6. Compared with the results in terms of testing points (Figure 10), the smallest* R*^{2} is greater than 0.7, which means that the accuracy of the constructed surrogate models is acceptable. Thus, the SVR model is eventually chosen as the surrogate model for GSA due to its highest accuracy.

##### 4.2. Global Sensitivity Analysis

A number of GSA methods have been developed and applied in practical application during the past decades. Among them, Sobol’s method is a variance-based Monte Carlo method that allows the computation of both the sensitivity indices of individual parameters,* S*_{i}, and the interactions between these parameters through the ratio of each sensitivity index to the corresponding total sensitivity index [21, 43]. GSA generally requires a huge number of Monte Carlo (MC) calculations to ensure its accuracy. In this work, 20,000, 40,000, 60,000, 80,000, and 100,000 MC simulations are used for the calculation of the Sobol indices, respectively. It is found that 100,000 MC simulations can ensure the convergence of the Sobol indices. Thus, 100,000 MC calculations are carried out for each objective of each metamodel in this case [44, 45].

###### 4.2.1. RMS Values of Vibration Displacement (RDef)

Figure 11 provides a guideline to the effects of 25 input parameters on the RMS values of vibration displacements (RDef) in three directions at the ring gear node by varying them over the uncertainty space, respectively. It can be found in Figures 11(a)–11(c) that the connected stiffness and the cutterhead load are the main sensitive parameters for RDef. Figure 11(a) shows the variation of RDef_A with respect to the five importer sensitive parameters, such as (*x*3), (*x*16), (*x*18), (*x*19), and (*x*22), which are chosen from Figure 11(a). In Figure 12(a), the RMS value of the axial vibration displacement increases with the increasing of and decreases as the connected stiffness increases. This illustrates that the vibration displacement can be decreased effectively by increasing the connected stiffness at the key points.

###### 4.2.2. Peak Vibration Acceleration (PAcc)

The sensitivities to the peak vibration acceleration at the node of the pinion in three directions are shown in Figure 13. In Figure 13(a), the damping ratio of steel (*x*4) is the most significant parameter with an individual effect of about 80%. While the cutterhead torque (*x*2) is the most significant parameter, Figure 13(c), the individual effect is about 60%. By varying the five greater sensitive parameters (*x*1), (*x*2), (*x*11), (*x*23), and (*x*24) over the uncertainty space, the variation of PAcc_L is obtained and shown in Figure 12(b). In Figure 12(b), the peak value of lateral vibration displacement increases with the increasing of and , and decreases with increasing , , and . This illustrates that the vibration acceleration can be decreased efficiently by increasing the connected stiffness at the key points.

###### 4.2.3. Dynamic Meshing Force (MF)

Figure 14 shows the GSA results of the dynamic meshing force. In Figure 14(a), the maximum dynamic meshing force (MF_M) of the pinion-ring gear mesh is only sensitive to the torque of cutterhead (*x*2), while the damping ratio of rock (*x*11) is the most significant parameter with an individual effect of about 80% for the dynamic deviation factor (MF_Se). This is because the meshing force of the pinion-ring gear mesh is directly caused by the cutterhead torque, and the dynamic deviation factor not only reflects the load distribution among the multiple driving elements, but also gives expression to the vibration amplitude of meshing force. Figure 15 shows the variation of MF_M and MF_Se with respect to the four importer sensitive parameters, such as (*x*2), (*x*8), (*x*11), and (*x*20), which are chosen from Figure 14. In Figures 15(a) and 15(b), by increasing the torque of cutterhead, both the maximum value and the dynamic deviation factor of the meshing force increase, and the variation amplitude of the latter is smaller than the former. This can be illustrated by (5) where both the dynamic force and the static force increase with the increasing of . Meanwhile, the dynamic deviation factor increases by increasing the damping ratio () and the connected stiffness (), where the variation rule is oppositely compared with the vibration displacements and acceleration.

From the results, it can be found the value of the individual effect is close to that of the total effect for each output, respectively, and the sum of the individual effects is almost 1.0. It means that the interactive effect of parameters on the output is relatively small, and the negligible parameters can be detected by the low values of the individual effects [44]. Thus, the ranking results only based on the individual effect in this section are convincing.

#### 5. Conclusions

In this work, a hierarchical dynamic model of the complex gear transmission system is proposed based on the generalized finite element modeling idea. Then, a procedure based on the combination of metamodels with global sensitivity analysis is proposed to study the importance ranking of the different sources of uncertainty during operation. The main contribution is concluded as follows:

(1) To analyze the complex dynamic characteristics of the TBM cutterhead driving system, a simulation-based GSA is carried out, where a hierarchical dynamic model of the complex gear transmission system is proposed based on the generalized finite element modeling idea. It has been shown that the proposed simulation-based GSA is feasible and efficient to rank the sensitive parameters of the driving system.

(2) Vibration responses under the action of simulated load are estimated and verified by comparing with the data measured on site. It can be concluded that the theoretical analysis model established in this paper is able to obtain the dynamic characteristics of the TBM cutterhead driving system as the error is sufficiently small.

(3) The vibration displacements can be decreased effectively by increasing the connected stiffness, especially by increasing the supporting stiffness of cutterhead-rock and the main bearing. The vibration acceleration can be decreased by increasing the connected damping ratio, especially the damping ratio of rock supporting. However, increasing the connected stiffness and the damping ratio will increase the dynamic deviation factor of the meshing force. Therefore, both of the influence of connected parameters on the dynamic responses and load uniformity among the multiple driving elements should be considered in the design process of TBM driving system.

It needs to be noted that the metamodel uncertainty has an important effect on sensitivity measures. Although the metamodel in this work is with high accuracy, the metamodel uncertainty still cannot be avoided. In the future, our research will focus on the effect of metamodel uncertainty on sensitivity measures and apply sensitivity analysis to other engineering systems in TBM.

#### Conflicts of Interest

There are no conflicts of interest related to this paper.

#### Acknowledgments

The research is supported by National Natural Science Foundation of China (Grant nos. 51505061 and U1608256).