Numerical Simulation of EPB Shield Tunnelling with TBM Operational Condition Control Using Coupled DEM–FDM

Abstract

This study demonstrates a three-dimensional numerical simulation of earth pressure balance (EPB) shield tunnelling using a coupled discrete element method (DEM) and a finite difference method (FDM). The analysis adopted the actual size of a spoke-type EPB shield tunnel boring machine (TBM) consisting of a cutter head with cutting tools, working chamber, screw conveyor, and shield. For the coupled model to reproduce the in situ ground condition, the ground formation was generated partially using the DEM (for the limited domain influenced by excavation), with the rest of the domain being composed of FDM grids. In the DEM domain, contact parameters of particles were calibrated via a series of large-scale triaxial test analyses. The model simulated tunnelling as the TBM operational conditions were controlled. The penetration rate and the rotational speed of the screw conveyor were automatically adjusted as the TBM advanced to prevent the generation of excessive or insufficient torque, thrust force, or chamber pressure. Accordingly, these parameters were maintained consistently around their set operational ranges during excavation. The simulation results show that the proposed numerical model based on DEM–FDM coupling could reasonably simulate EPB driving while considering the TBM operational conditions.

1. Introduction

The versatility of an earth pressure balance (EPB) shield has led it to become the most popular type of tunnel boring machine (TBM) for use in various ground conditions around the world [1,2]. During tunnelling, an EPB shield uses the excavated soil to fill the working chamber and stabilize the tunnel face. The support pressure of the tunnel face is attained by applying pressure to the soil that fills the chamber and screw conveyor, thus preventing a sudden inflow of soil and groundwater into the machine. This distinctive tunnelling mechanism directly uses muck as a supporting material, and so the dynamic interaction between the machine and the excavated soil is more meaningful in EPB shield tunnelling than in other tunnelling methods. The mechanical relationship between the machine and surrounding soil is generally observed indirectly in real time during tunnelling via data from the TBM. These include the cutter head torque, thrust force, chamber pressure, quantity and quality of discharged muck, and abrasion of cutting tools. Analyzing these machine data is considered essential for assessing the performance of EPB shields.

Collecting and interpreting machine data recorded during tunnelling is especially important to the machine operators who must carefully maintain a TBM within certain operational ranges. The appropriate operational ranges for a TBM are commonly determined by a machine superintendent and provided to the operator before each day’s work. Maintaining the optimal operational conditions during tunnelling is known to enhance significantly the overall performance of a TBM [3]. Accordingly, efficient tunnelling relies on there being sufficient machine data to allow engineers and operators to analyze carefully an interaction of the machine with the soil. However, collecting enough machinery data of sufficient quality in the field for analysis is difficult, because the number of sites around the world where data can be obtained is limited and the data are generally not published. A further consideration is that the different conditions of each tunnelling site complicate the direct comparison of data. Nevertheless, many studies regarding TBM operation have been conducted using machine data accumulated from tunnelling sites, and there always has been demand for large quantities of machine data for TBM operators, engineers, and researchers [4,5].

The limitations of collecting field data make numerical analysis a cost- and time-effective alternative for gathering and reproducing TBM tunnelling data. Many numerical approaches have been suggested to simulate tunnel construction by both conventional and machine methods. The most representative numerical methods for tunnelling simulation are finite difference methods (FDMs) and finite element methods (FEMs). The former have been used predominantly to evaluate ground displacement induced by tunnelling [6,7,8,9]. A recent convergence study conducted inverse analysis of TBM parameters with an artificial neural network [10]. Tunnelling studies have also used FEMs in various ways: for example, simulation of shield tunnelling [11,12], suggestion of critical face pressure and backfill pressure [13], dynamic load analysis on a cutterhead [14], and muck flow simulation in a EPB shield [15]. However, both FDMs and FEMs face challenges regarding the simulation of the interaction between a TBM and the soil in real time. That is, it is hard to obtain the TBM data that are significant for analysis of the dynamic interaction between machinery and the soil.

To overcome the drawbacks of these methods, a discrete element method (DEM) is now widely adopted to employ the recent rapid advances in computing power in simulating tunnelling [16,17,18,19]. With a number of discrete elements (e.g., particles, clumps, and blocks) whose motion follows Newton’s second law and is calculated by an explicit time-stepping scheme, DEM can simulate the behavior of a granular and discrete medium like soil and rock in geotechnical engineering [20,21]. DEM represents a straightforward approach to modeling ground behavior and, therefore, could be applied to collecting real-time TBM data concerning the dynamic interaction of the soil and structure during shield tunnelling [22,23,24,25]. However, the computational cost increases considerably as the number of discretized elements increases. Moreover, to apply macroscopic parameters (e.g., Druker–Prager, Mohr–Coulomb, and modified Cam–Clay parameters), microscopic parameters (i.e., contact parameters) must be calibrated before DEM modeling. Namely, the contact model and contact parameters must be determined via preliminary numerical analyses. Geotechnical engineering commonly uses direct shear [26], uniaxial [27,28], and triaxial [29] tests to calibrate the contact parameters. However, these preliminary analyses might magnify numerical errors and make the analysis time consuming.

The computational load caused by numerous elements and the additional calibration process are critical barriers to applying DEM in practical engineering problems. Therefore, to mitigate the limitation, DEM is frequently coupled with other numerical methods such as a boundary element method (BEM) [30,31], computational fluid dynamics (CFD) [32], FDM [33], and FEM [34,35]. The coupling of DEM and FDM has received notable attention for simulating shield tunnelling, because both methods adopt the same explicit time-stepping, which facilitates data exchange at each time step. Using DEM–FDM coupling, Yin et al. [36] examined the progressive failure of a tunnel face during shield tunnelling. To assess the effect of particle shape on the tunnel face failure region, settlement and face supporting pressure were studied. Qu et al. [37] also adopted coupled DEM–FDM to simulate scaled EPB shield tunnelling in sandy ground. The analysis applied a conditioning technique with particles having different friction coefficients. The thrust force, cutter head torque, and conditioned muck distribution in the machine were investigated with varying soil conditioning.

In accordance with the current demand for, and progress in, the real-time numerical simulating of EPB shield tunnelling, this study couples DEM and FDM to establish a numerical model for EPB shield tunnelling. The numerical code is established using commercial numerical programs: Particle Flow Code in three dimensions (PFC3D, Itasca Consulting Group, Inc.) for DEM [38] and Fast Lagrangian Analysis of Continua in three dimensions (FLAC3D, Itasca Consulting Group, Inc.) for FDM [39]. Virtual ground for excavation during tunnelling is formulated using both DEM and FDM. Generation of the excavation area by DEM uses numerical calibration to determine the contact parameters of particles. The rest of the ground domain (outside the excavation area) consists of FDM grids to alleviate the computation load. A spoke-type TBM (diameter, 6.14 m) comprising a cutter head with cutting tools, working chamber, screw conveyor, and shield is then imported into the model. Within the model, data for the TBM (e.g., torque, thrust force, chamber pressure, soil discharge, and settlement) are monitored as it advances. To verify the applicability of the proposed model, the analysis is conducted by controlling the TBM operational conditions; the simulated TBM data are presented here and discussed.

2. Numerical Model for Earth Pressure Balance (EPB) Shield Tunnelling Simulation

The numerical modeling of EPB shield tunnelling first needs a virtual ground to be excavated. This is 23 m-deep, homogeneous, weathered granite soil formed above bedrock. The center of the TBM is driven around 18.9 m below the ground level.

Table 1. Ground conditions of numerical model for earth pressure balance (EPB) shield tunnelling simulation.

Property	Value
Type of soil	Weathered granite soil
Dry unit weight [kN/m3]	19.5
Internal friction angle, $\phi$ [°]	27.3
Cohesion, $c$ [kPa]	26.3
Lateral earth pressure coefficient at rest, $K_0$ [-]	0.5
Thickness of soil layer [m]	23
Cover depth [m]	15.8
Diameter of excavation area [m]	6.2

summarizes the ground conditions for the simulation, including the soil type and properties, ground dimensions, and diameter of the excavation area. Figure 1 depicts the model set-up. To apply the DEM–FDM coupling scheme, the ground formation is divided into two domains: one for each method. The DEM domain is modeled with ball and wall elements using PFC3D. Through this domain, the EPB shield model advances with an excavation area 6.2 m in diameter; tunnelling can be conducted by the EPB shield model with a diameter of 6.14 m. Meanwhile, the FDM domain is composed of zone elements (grids) using FLAC3D to surround the DEM domain. This domain reduces the number of ball elements (particles) by substituting them with zone elements. This reduction in the number of ball elements can significantly decrease computational cost. Furthermore, it is possible to apply horizontal in situ stress effectively through a wall–zone coupling method.

As mentioned, DEM modeling requires the contact model governing translational and rotational motion of particles to be determined for the input macroscopic strength parameters. In this study, the governing contact model has to contain an adhesion component as the soil type in the excavation area (weathered granite) has weak cohesion. Among the contact models built into the PFC3D program, the bonding-based contact model (linear contact and parallel bond model) is the most widely applied, especially for simulating the attractive force between particles for rock mechanics [40,41,42,43]. However, it cannot reproduce the cohesion of soil, because once adhesion is destroyed by the separation of particles, the attractive interparticle force is not regenerated. Therefore, this study simulates the behavior of cohesive soil using the adhesive rolling resistance linear contact model, which provides the linear elastic, frictional (including rolling), and adhesive interparticle behaviors [44,45].

Figure 2 outlines this contact model. Figure 2a shows activation and deactivation of contact and adhesion depending on the surface gap ($g_s$) and specified attraction range ($D_0$). Figure 2b depicts the components of the contact model, specifically the contact force ($F_c$) and contact moment ($M_c$) between two ball elements. The contact force is the summation of linear force ($F^l$), dashpot force ($F^d$), and attractive force ($F^a$) expressed with normal and shear components (labeled with subscripts $n$ and $s$, respectively): ($F_c=F^l+F^d+F^a$); the contact moment is equal to the rolling resistance moment ($M_r$): ($M_c=M^r$). The contact force and moment depend on the following governing contact parameters: normal stiffness ($k_n$), shear stiffness ($k_s$), friction coefficient ($\mu$), reference gap ($g_r$), effective modulus ($E^{*}$), rolling friction coefficient (${\mu }_r$), normal critical damping ratio (${\beta }_n$), shear critical damping ratio (${\beta }_s$), maximum attraction force ($F_0$), and attraction range ($D_0$). These parameters have to be calibrated before application in the model. In this study they are determined in several large-scale triaxial tests. Lastly, Figure 2c gives details of the adhesion component, showing the linear approximation of the relationship between the attractive force (like the van der Waals interaction) and the surface gap (i.e., separation).

2.1. EPB Shield Model

Modeling EPB shield tunnelling uses a spoke-type TBM of 6.14 m diameter. All its components, including the cutter head with cutting tools, working chamber, screw conveyor, and shield, are imported from 3D model data (STL; Standard Triangle Language file format). The imported data are transformed into wall elements in PFC3D. Unlike a ball element, a wall element is suitable for simulating TBM operation, because the programmer can specify the machine’s translational and rotational motion during time-stepping.

Table 2. Dimensions of EPB shield model for numerical simulation.

Property	Value
Type of cutter head	Spoke-type
Diameter of cutter head [m]	6.14
Number of spokes	6
Number of cutting tools	91
Opening ratio [%]	75.5
Type of screw conveyor	Shaft-type
Diameter of screw [m]	0.9
Length of screw conveyor [m]	12
Pitch of screw [m]	0.5
Inclination of screw conveyor [°]	28
Length of shield [m]	8

lists the dimensions of the imported EPB shield model. The 91 cutting tools on the spoke-type cutter head consist of 60 tooth bits, 24 precut bits, six shell bits, and one fishtail bit, as shown in Figure 3a,b, which also gives the dimensions of the EPB shield model in side and sectional views.

2.2. Calibration of Contact Parameters

As a series of calibration analyses have to be completed before generating particles for the designated numerical model, large-scale triaxial tests were conducted to determine the contact parameters [46,47].

Table 3. Details of large-scale triaxial test model for calibration of contact parameters.

Property	Value
Dimension of specimen, $D$ × $H$ [m]	2 × 4
Packing of ball element	Hexagonal close packing
Diameter of ball element [cm]	20
Number of ball element	2101
Density of ball element [kg/m3]	2700
Contact model	Adhesive rolling resistance linear
Number of shell element	576
Thickness of shell element [mm]	6
Density of shell element [kg/m3]	90
Elastic modulus of shell element [kPa]	1400
Confining pressure [kPa]	125, 250, 375, 500
Strain rate [%/min]	0.03

lists details of the established triaxial test model. To apply the calibration results to the actual-size EPB shield tunnelling model, a relatively large specimen (diameter, 2 m; height, 4 m) and particle size are considered. Therefore, the soil specimen is modeled with ball elements with a single diameter of 20 cm. When DEM modeling with a single size of particle, their packing is a significant factor, because it determines the behavior of particle assemblies. Accordingly, this study generates a specimen with hexagonal-close-packed balls. After that, shell elements of FLAC3D are generated around the ball elements to act as membranes, which transmit the confining pressure to the specimen before and during the test.

Figure 4 shows the configuration of the simulated triaxial test model with the generated ball and shell elements. Analyses use different confining pressures of 125, 250, 375, and 500 kPa and a constant strain rate of 0.03%/min. With the fixed bottom plate, the load is applied by the top plate (consisting of wall elements) until 80% of the maximum deviatoric stress is observed. The analysis is repeated by trial and error with different combinations of contact parameters until the designated Mohr–Coulomb strength parameters, internal friction angle ($\phi$), and cohesion ($c$) are obtained from a Mohr–Coulomb failure envelope plot.

Figure 5 presents the results of the triaxial tests. Figure 5a shows four stress–strain curves at different confining pressures: the maximum deviatoric stress increases with increasing confining pressure. Low confining pressure induces plastic behavior with softening, while brittle behavior occurs at high confining pressure. The plotted Mohr–Coulomb failure envelope (Figure 5b) from the stress–strain relationship allows $\phi$ and $c$ to be attained by trial and error. The absolute errors between the target and calibrated Mohr–Coulomb parameters are under 0.5%.

Table 4. Calibrated parameters after large-scale triaxial tests.

Calibrated Contact Parameters of Adhesive Rolling Resistance Linear Contact Model
Parameter		Value
Normal stiffness, $k_n$ [MN/m]		4.5
Shear stiffness, $k_s$ [MN/m]		3
Friction coefficient, $\mu$ [-]		0.07
Reference gap, $g_r$ [mm]		0
Effective modulus, $E^{*}$ [MN/m2]		10
Rolling friction coefficient, ${\mu }_r$ [-]		0.03
Normal critical damping ratio, ${\beta }_n$ [-]		0.05
Shear critical damping ratio, ${\beta }_s$ [-]		0.05
Maximum attractive force, $F_0$ [N]		0.01
Attraction range, $D_0$ [mm]		2.5
Calibrated Mohr–Coulomb Failure Criterion Parameters
Parameter	Calibrated Value		Target Value	Error [%]
Internal friction angle, $\phi$ [°]	27.22		27.3	0.29
Cohesion, $c$ [kPa]	26.17		26.3	0.49

summarizes the determined contact parameters of the adhesive rolling resistance linear contact model and Mohr–Coulomb failure criterion parameters. The obtained calibrated values are consequently applied to the DEM domain of ground formation.

2.3. Ground Formation with Coupled Discrete Element Method-Finite Difference Method (DEM-FDM)

The DEM–FDM coupling program can generate a ground formation for tunnelling using the following codes for PFC3D and FLAC3D: (a) one-dimensional structural element coupling, (b) wall–zone coupling, and (c) ball–zone coupling. This study mainly adopts the wall–zone coupling scheme. This method allocates the wall elements of PFC3D along the zone faces of FLAC3D. The created walls allow the contact force and moment applied by discrete elements on the walls to be transmitted to gridpoints of the zones. Conversely, information about the zone gridpoints (e.g., force, displacement, and velocity) can be delivered not only to the wall elements but also to ball elements. The program guide gives further details of the coupling methods [38,39]. To reproduce the initial stress condition in the ground model based on wall–zone coupling, the model has to be gravitationally stabilized before installation of the EPB model. The ground is formulated and stabilized as follows.

(1)

Wall compaction. Within the specified DEM domain, after the ball elements are generated using the conditions determined by preliminary calibration (

Table 5. Summary of input parameters for ground formation.

DEM Domain
Property	Value
Packing of ball element	Hexagonal close packing
Diameter of ball element [m]	0.1
Number of ball element	59,024
Density of ball element [kg/m3]	2700
Contact model of ball element	Adhesive rolling resistance linear
Contact parameters	Refer to Table 4
FDM Domain
Property	Value
Size of zone element, $W$ × $H$ [m]	0.5 × 0.5
Number of zone element	50,432
Density of zone element [kg/m3]	1987.77
Constitutive model of zone element	Mohr–Coulomb failure criterion
Elastic modulus, $E$ [MPa]	24
Poisson’s ratio, $v$ [-]	0.3
Internal friction angle, $\phi$ [°]	27.3
Cohesion, $c$ [kPa]	26.3
Lateral earth pressure coefficient at rest, $K_0$ [-]	0.5

), 6 walls (12 facets) are installed around the ball assembly. Excluding the side where tunnelling is initiated and the bottom wall, the designated loads are applied through five walls (one top and four side walls) with a servo-mechanism. On the top wall, vertical earth pressure is loaded considering the cover depth. On the four side walls, the average lateral earth pressure is applied considering the lateral earth pressure coefficient at rest ($K_0$ = 0.5; Figure 6a,b).

(2)

Wall–zone compaction. Zone elements are generated around wall-compacted balls to simulate the rest domain of the ground formation. Along zone faces of the FDM and DEM domain boundaries (i.e., wall–zone boundaries), walls are created by wrapping zone faces. After that, input parameters (Table 5) are applied to zones and time-stepping is conducted with gravitational acceleration until the convergence criterion is satisfied. During wall–zone compaction, a quarter of the target elastic modulus ($E$ = 6 MPa) is applied as zones should be sufficiently deformed to guarantee efficient transmission of stress into the DEM domain. Upon completion of ground stabilization, the elastic modulus is modified to the designated value ($E$ = 24 MPa; Figure 6c,d).

After the ground model is formulated and stabilized, the coupling between the DEM and FDM domains is verified by measuring the in-situ stress with depth. The stress of the FDM domain is obtained by reading the values of the zone gridpoints. The DEM domain stress is measured by installing measurement balls with a diameter of 0.8 m, four times that of the ball elements. Figure 7 illustrates the measured vertical and horizontal in-situ stresses (with designated $\pm$5% allowable error) along with the exact (analytical) solutions of in situ stress. The stress distribution in the FDM domain appears almost identical to the exact solution. After wall compaction, the horizontal stress distribution in the DEM domain exceeds the allowable error with a relatively steep gradient of stress with depth. However, after wall–zone compaction, the horizontal stress is within the allowable error, suggesting that ground modeling with DEM and FDM coupling is reasonably complete.

2.4. Numerical Model and Analysis

The EPB shield tunnelling model is completed by installing the EPB shield model within the stabilized ground formation. Before importing the machine model, the initiation of the tunnelling side wall is replaced with a wall with a hole for tunnelling. After that, the shield model is inserted up to the depth of the bulkhead wall, because unless the chamber is filled with prestressed balls, abrupt settlement at the DEM-FDM boundary will occur, and the coupling will be broken. In addition, when importing the EPB shield model, balls overlapping the TBM are eliminated to reduce numerical errors caused by the overlapping. Figure 8 presents the final configuration of the EPB shield tunnelling model with coupled DEM–FDM.

The model can simulate EPB shield tunnelling under the designated TBM operational conditions. The penetration rate and rotational speed of the cutter head and screw conveyor can be manually assigned. Moreover, TBM machine data (torque, thrust force, chamber pressure, soil discharge, and ground settlement) can be gathered as the TBM advances. The torque is calculated as the summation of the vector product of the contact force vector and the position vector from the rotation axis of all cutter head elements. The sum of the contact forces in the x-direction of all EPB shield components is regarded as the thrust force. The chamber pressure is estimated as the average stress in the x-direction of the bulkhead wall, as the face pressure is generally estimated with several pressure meters installed on the bulkhead wall. Soil discharge is deduced from the number of ball elements in the model domain by deleting balls discharged outside of the model domain. Lastly, ground settlement is measured by reading the displacement of zone gridpoints.

This study considers several analysis cases to examine the applicability of the simulation model for gathering the appropriate machine data. In addition to tunnelling under constant operational conditions, this study focuses on simulating tunnelling while controlling the TBM operational conditions. Before a simulation can control the operational conditions, their appropriate operational ranges have to be defined. Accordingly, those for torque, thrust force, and chamber pressure are determined using the following empirical and theoretical relationships.

Equation (1) gives an empirical formula widely adopted by tunnel engineers for estimating the required cutter head torque [48,49].

$$T=\alpha D^3,$$

(1)

where $T$ is the required torque (kNm), $\alpha$ is the unitless coefficient of proportionality for the torque, and $D$ is the diameter of the TBM (m).

The proportional relationship between required thrust force and TBM area is well established, because the thrust force is generally designed based on the area of the TBM [50]. Therefore, the relationship is as in Equation (2).

$$F_T=\beta \frac{\pi D^2}{4},$$

(2)

where $F_T$ is the required thrust force (kN) and $\beta$ is the unitless coefficient of proportionality for the thrust force.

When tunnelling with EPB shields, the face pressure is achieved by pressurizing the muck within the chamber (i.e., face pressure equals the chamber pressure). Much research has sought to determine the adequate face pressure [51]; conventionally, the chamber pressure should be between Rankine’s active lateral earth pressure and the lateral earth pressure at rest for preventing surface heaving or settlement, as shown in Equations (3) and (4).

$$P_{c,min}=K_a{\sigma }_v{}^{\prime },$$

(3)

$$P_{c,max}=K_0{\sigma }_v{}^{\prime },$$

(4)

where $P_{c,min}$ is the minimum chamber pressure, $P_{c,max}$ is the maximum chamber pressure, $K_a$ is the Rankine’s active lateral earth pressure coefficient, $K_0$ is the lateral earth pressure coefficient at rest, and ${\sigma }_v{}^{\prime }$ is the vertical effective stress.

Considering the above relationships and model conditions, this study adopts the following operational ranges for simulating the TBM operational condition control: required torque, 2000–2500 kNm; required thrust force, 8500–11,000 kN; chamber pressure, 138.5–185 kPa.

Along with the specified operational range, the penetration rate and rotational speed of the screw conveyor can be automatically adjusted within the set operational range by a user-defined subroutine in program codes that control the TBM operational conditions in response to the monitored TBM data. The rotational speed of the cutter head is fixed (here, at ${\omega }_{cutter}$ = 2 rpm), as is common during shield tunnelling. During the simulation, the penetration rate is controlled by monitored torque and thrust force: i.e., if excessive or insufficient torque or thrust force are recorded, the penetration rate automatically increases or decreases by 20%. The rotational speed of the screw conveyor is adjusted by 20% in response to the measured chamber pressure being outside of the set range. To prevent unlimited variation of the operating conditions, lower and upper bounds are set for the penetration rate (0.2 and 2 mm/s) and the rotational speed of the screw conveyor (2 and 20 rpm). To evaluate the numerical model and identify its applicability to simulating TBM operation, five analyses are conducted and compared.

Table 6. Numerical simulation cases with different TBM operational conditions.

Case Number	Rotational Speed of Cutter Head, ${\mathit{\omega }}_{\mathit{c}\mathit{u}\mathit{t}\mathit{t}\mathit{e}\mathit{r}}$ [rpm]	Penetration Rate, PR [mm/s]	Rotational Speed of Screw Conveyor, ${\mathit{\omega }}_{\mathit{s}\mathit{c}\mathit{r}\mathit{e}\mathit{w}}$ [rpm]
1	2	0.5	5
2			10
3		1	5
4			10
5		Automatically controlled (initial $PR$ = 0.5 mm/s and ${\omega }_{screw}$ = 10 rpm)

lists the TBM operational conditions in each simulation.

3. Discussion of Analysis Results

This section discusses the applicability of the proposed model and the numerical analysis results for the measured torque and thrust force, the monitored chamber pressure, the discharged volume and discharge rate of particles, and the surface settlement.

The measured torque and thrust force are considered the criteria for estimating TBM performance. Their being below the specification of the TBM signifies that more rapid tunnelling could be possible than is currently underway. That is, an increase in the penetration rate could be possible by increasing the jacking force within the capacities of the torque and thrust force. Conversely, excessive torque and thrust force represent increased risk of downtime. Therefore, the penetration rate has to be reduced to bring them within acceptable levels.

When examining the torque data here, the raw data are modified to facilitate control of the operational conditions and comparison of varying data. Figure 9 presents the collected raw and modified torque data during 200 mm of TBM advance under the specific operational conditions of ${\omega }_{cutter}$ = 2 rpm, $PR$ = 0.5 mm/s, and ${\omega }_{screw}$ = 15 rpm. The figure shows the fluctuations of the raw torque data to be so large as to hinder regulation of the operational conditions within the defined bounds. Analyzing the overlapping data for comparison is also difficult. Therefore, a moving average scheme ($n$ = 11) is adopted to adjust the raw data.

Meanwhile, the chamber pressure and discharge of particles are assessed in terms of tunnelling safety, because they are directly related to tunnel face stability and ground settlement. Excessive chamber pressure during tunnelling will cause heaving, which is severely damaging to existing structures. Conversely, insufficient chamber pressure will lead to an abrupt soil inflow causing settlement. Therefore, in addition to the chamber pressure and discharging data, the measured surface settlement in each analysis case is compared after completion of the TBM advance.

All analyses consider a 300 mm advance of the EPB shield. The cases with constant operational conditions are carried out for 300 s (cases 1 and 2) or 600 s (cases 3 and 4). In contrast, the analysis with automatic control of the operational conditions (case 5) has four adjustments during the 300 mm advance, which takes about 422.84 s. After an advance of 48.96 mm, the measured thrust force reaches its lower bound, and the penetration rate reaches its upper bound (2 mm/s). Next, at 164.15 mm, the chamber pressure reaches its lower bound, causing the rotational speed of the screw conveyor to reach its lower bound (2 rpm). The third adjustment comes at 255.69 mm, when the torque hits its upper bound, inducing the penetration rate to drop to 0.2 mm/s. The last adjustment at 287.87 mm is in response to the chamber pressure reaching its upper bound, causing the rotational speed of the screw conveyor to rise to 20 rpm. Finally, an advance of 300 mm is completed with $PR$ = 0.2 mm/s and ${\omega }_{screw}$ = 20 rpm at 422.84 s.

Table 7. Variation of TBM operational conditions during EPB shield tunnelling simulation (case 5).

Advance [mm]	Time [s]	Penetration Rate, PR [mm/s]	Rotational Speed of Screw Conveyor, ${\mathit{\omega }}_{\mathit{s}\mathit{c}\mathit{r}\mathit{e}\mathit{w}}$ [rpm]	Note
0	0	0.5	10	Initiation of tunnelling
48.96	97.92	2	10	Lower bound of thrust force
164.15	155.52	2	2	Lower bound of chamber pressure
255.69	201.29	0.2	2	Upper bound of torque
287.87	362.19	0.2	20	Upper bound of chamber pressure
300	422.84	0.2	20	End of tunnelling

lists the variations of operational conditions during simulation case 5.

3.1. Torque and Thrust Force

Figure 10a,b present torque and thrust force data acquired during EPB shield tunnelling simulation. When the rotational speed of the screw conveyor is relatively fast compared with the penetration rate, the torque and thrust force decrease as tunnelling advances (case 2). Conversely, tunnelling with a rapid penetration rate compared with the rotational speed of the screw conveyor leads to steadily increasing torque and thrust force (case 3). However, if the proper penetration rate and rotational speed of the screw conveyor are applied, torque and thrust force are maintained without large variations (cases 1 and 4). During automatically controlled operation, the torque and thrust force settled within the defined operational range in response to variation of the operational conditions (Table 7, case 5). Figure 11a,b give the average torque and thrust force after the first 100 mm of advance to eliminate the effect of initial unstabilized data. With reference to the results in Figure 10, to increase the rate of TBM advancement (case 2) or to reduce unexpected downtime (case 3), the operational conditions should be optimized by moderating torque and thrust force during tunnelling (cases 1, 4, and 5). Accordingly, the simulated torque and thrust force data indicate that determining the optimal penetration rate and rotational speed of the screw conveyor is important for enhancing the performance of a TBM.

3.2. Chamber Pressure, Discharge, and Surface Settlement

Figure 12 illustrates the variations in the monitored chamber pressure during the EPB shield advance. Each case shows a decreasing trend from the large initial chamber pressure, because the chamber is filled with prestressed balls. Figure 13 compares the average chamber pressures after the first 100 mm of advance, when the initial variations are stabilized. Contrary to the torque and thrust force results, the chamber pressure is maintained within the set operational range only under automatic control of the penetration rate and rotational speed of the screw conveyor (case 5). These results show that when tunnelling with insufficient chamber pressure (case 2), tunnel face instability and large surface settlement can be expected. Conversely, when the chamber pressure is large compared with the estimated pressure (cases 1, 3–5), heaving of the surface will occur.

During simulations, particles inside the chamber are successfully discharged by the rotation of the screw conveyor. Figure 14 compares the discharged volume of particles and the discharge rate with time. As the screw conveyor is unfilled at the beginning of the analysis, there is a time delay in discharge data acquisition. The figure shows that the discharged volume and discharge rate are closely related to the rotational speed of the screw conveyor, with the latter being exactly proportional to the conveyor’s rotation. That is, discharge rates of 0.004, 0.011, 0.02, and 0.04 m³/s occur at rotational speeds of 2, 5, 10, and 20 rpm, respectively. When the operational conditions are automatically adjusted, the discharge rate varies twice with changes in rotational speed of the screw conveyor (case 5). Figure 15 compares the total discharged volume after the TBM has advanced 300 mm. Although direct comparison of the discharged volume is impossible because of differences in the analysis time among the cases, it can be inferred that a large discharge will destabilize the tunnel face and cause settlement (case 2) and that a small discharge will induce surface heaving (case 3).

Figure 16a,b give the transverse and longitudinal surface settlement trough after the completion of tunnelling. As expected from the chamber pressure data described above, heaving occurs due to excessive chamber pressure in cases 1, 3, and 4. In contrast, relatively large settlement with insufficient chamber pressure is observed in case 2. However, the settlement in case 5 is successfully suppressed compared with that in case 2. As noted above, these results can also be deduced from the discharged volumes of particles in Figure 15. Even though the settlement results are not exact, the numerical results confirm the significance of maintaining appropriate chamber pressure and discharge in controlling surface settlement. That is, both are key factors to guarantee safe tunnelling with an EPB shield.

4. Conclusions

A coupled DEM–FDM numerical model simulated EPB shield tunnelling. Its main components are an actual-size spoke-type EPB shield and ground formation. Modeling used large-scale triaxial test analyses for calibration and a two-step ground formulation scheme with DEM–FDM coupling. Subsequent tunnelling simulations adopted varying penetration rate and rotational speed of the screw conveyor, including their automatic control during operation. Torque, thrust force, chamber pressure, discharge, and settlement data were obtained, and their implications were discussed. The conclusions of this numerical study are summarized as follows:

• The numerical model successfully simulated EPB shield tunnelling under control of the operational conditions. During an advance of 300 mm, the penetration rate and rotational speed of the screw conveyor were varied four times in response to the monitored TBM data. The adjustments ensured that the measured torque, thrust force, and chamber pressure remained within their specified operational ranges.
• The simulated torque and thrust force data demonstrated the significance of adopting appropriate penetration rate and rotational speed of the screw conveyor to enhancing the TBM performance. When the torque and thrust force were maintained within their specified operational ranges, the TBM advance was optimal within its specifications.
• The chamber pressure and discharge data were directly related to the safety of tunnelling. Insufficient or excessive chamber pressure could induce surface settlement or heaving. However, controlling the chamber pressure and discharge during operation considerably suppressed surface settlement. Therefore, maintaining appropriate chamber pressure and muck discharge is significant for safe tunnelling with an EPB shield.
• The provided numerical model could be a useful tool for gathering simulated real-time TBM data and finding the optimal TBM operational conditions for efficient and safe tunnelling. In future, with some modifications, the model can be used in automatic TBM driving, simulating multiple-layer soil tunnelling, and designing cutter heads.