Discrete-Element Analysis of the Excavation Performance of an EPB Shield TBM under Different Operating Conditions

: This study used a discrete-element analysis to predict the excavation performance of a 7.73 m-diameter earth pressure balance (EPB) shield tunnel boring machine (TBM). The simulation mainly predicted several excavation performance indicators for the machine, under different operating conditions. The number of particles in the chamber and the chamber pressure varied, as the operating conditions changed during the simulated TBM excavation. The results showed that the compressive force, torque, and driving power acting on the TBM cutterhead varied with its rotation speed, increasing as the cutterhead rotation speed rose. The overall compressive force acting on all of the disc cutters and their impact wear increased linearly as the cutterhead rotation accelerated. The position of a disc cutter on the cutterhead had a particularly strong inﬂuence, with higher compressive forces experienced by the cutters closer to the center. In contrast, the gauge disc cutters at the transition zone of the cutterhead showed more wear than those elsewhere. The muck discharge rate and the driving power of the screw conveyor rose with increasing screw conveyor and cutterhead rotation speeds. Finally, this study suggests optimal operation conditions, based on pressure balance and operational management of the TBM.


Introduction
A tunnel boring machine (TBM) can be classified by its use of a shield, method of securing reaction force, face support, cutterhead type, and excavation method [1]. As a closed-type TBM, the earth-pressure-balanced (EPB)-shield TBM is widely used in urban tunnels, where rapid excavation is required with high safety standards.
An earth pressure balanced (EPB) shield is generally used on TBMs excavating soft ground. Ideally, the support pressure of the EPB shield can be controlled by the screwconveyor rotation speed and the TBM advance rate. However, most natural ground conditions are not ideally suited to EPB drives, and soils are often conditioned by the injection of water, polymers, and foams, to allow EPB shields to perform well. In addition, the pressure of the tunnel face is balanced with minimal surface settlement, as the chamber of the machine is filled with conditioned excavated materials (soil pastes) that act as a support medium at the tunnel face.
Successful excavation relies on the machine being optimized for the project conditions, which requires prediction of its excavation performance. Typical prediction methods include laboratory testing [2], analysis of field data [3][4][5][6], theoretical approach [7,8], and numerical analysis [9][10][11][12][13][14][15][16][17][18][19][20][21][22][23]. Among them, numerical analysis using the discrete-element method (DEM) is increasingly used in this field [9][10][11][12][13][14][15][16][17][18][19][20][21][22][23], owing to rapid advances in computing power. Unlike the TBM machine data, ground information for the previous project was not fully collected. Therefore, the tunneling site chosen for ground modeling comprised weathered soil with a cohesion of 26.30 kPa and an internal friction angle of 27.30 • [24]. Based on soil shear-strength design parameters, three-dimensional direct shear tests were simulated to determine the DEM parameters for each particle and its interactions, to verify its geotechnical design parameters. Table 2 lists the DEM parameters of each particle and its interaction parameters. Given this study's simulation time and purpose, 100-mm single spheres were used. Contact forces were calculated using the Hertz-Mindlin contact model developed by various studies [25][26][27][28][29][30] cited elsewhere [31]. Simulated, 3D direct-shear testing estimated the virtual ground model's shear strength. The tests applied each particle's DEM parameters and interactions using two large, frictionless box geometries of size 4 m (width) × 4 m (length) × 0.75 m (height). Figure 1a shows 26,317 randomly generated particles stacked by gravity. Tests were conducted until the shear strain reached 15% under five different normal load conditions (98,147,196,245, and 294 kPa), with a shear displacement rate of 0.01 m/s (Figure 1b). The servo control system controlled the normal load for coupling work with the DEM simulation. During simulation, the reaction force estimated from the lower box in the shear displacement direction was used to calculate the shear stress, considering the changing contact surface area between two boxes. The geotechnical design parameters were compared with the numerical simulation results, in Figure 2, based on the Mohr-Coulomb failure criteria and linear regression. The cohesion and internal friction angle estimated from the numerical simulation were about 25.42 kPa and 26.72°, respectively. Overall, the numerically estimated cohesion and internal friction angle were underestimated by about 3.3% and 2.1%, respectively, relative to those based on geotechnical design parameters. For the calibration work, this study set the shear speed to 0.01 m/s, which was much faster than the standard shear displacement rate. As a faster shear displacement rate means a higher estimated shear stress, the expected shear strength of the ground model would be lower than that in the simulation. The geotechnical design parameters were compared with the numerical simulation results, in Figure 2, based on the Mohr-Coulomb failure criteria and linear regression. The cohesion and internal friction angle estimated from the numerical simulation were about 25.42 kPa and 26.72 • , respectively. Overall, the numerically estimated cohesion and internal friction angle were underestimated by about 3.3% and 2.1%, respectively, relative to those based on geotechnical design parameters. For the calibration work, this study set the shear speed to 0.01 m/s, which was much faster than the standard shear displacement rate. As a faster shear displacement rate means a higher estimated shear stress, the expected shear strength of the ground model would be lower than that in the simulation.

TBM Numerical Model
Drawings of the 7.73 m EPB shield TBM from the Japanese company were used for the 3D TBM model. The TBM model was created with the 3D AUTO CAD software, and comprises six parts-cutterhead, disc cutters, cutter bits, chamber, shield, and screw conveyor ( Figure 3). The opening ratio of the cutterhead was ~24.25%. Table 3 lists the specifications for each part. The interaction parameters between the TBM and the particles were not considered in this study. Each disc cutter on the cutterhead face was fixed without rotation during excavation.

TBM Numerical Model
Drawings of the 7.73 m EPB shield TBM from the Japanese company were used for the 3D TBM model. The TBM model was created with the 3D AUTO CAD software, and comprises six parts-cutterhead, disc cutters, cutter bits, chamber, shield, and screw conveyor ( Figure 3). The opening ratio of the cutterhead was~24.25%. Table 3 lists the specifications for each part. The interaction parameters between the TBM and the particles were not considered in this study. Each disc cutter on the cutterhead face was fixed without rotation during excavation. Eight plate geometries are created and placed in the form of a box to be filled with particles for the ground model. One of the plates located on the +x-axis was designed with a hole of the same diameter as the TBM. Once all plates were arranged, the TBM geometry was initially introduced into the hole by nearly 400 mm. About 400,000 particles were then randomly generated and stacked by gravity. A further 52,000 particles were generated in the chamber and screw conveyor to represent TBM excavation. The top and bottom of the ground boundaries were set as parallel boundaries. Once these works were complete, an additional, fixed, overburden load of 265 kPa was applied. It was also applied continuously to the top plate by the servo control during excavation. As a result, a rectangular parallel-piped ground model of size 1.5 m (width) × 15.6 m (length) × 14.6 m (height) was created over the stabilization time ( Figure 4). Based on the EPB shield TBM's technical specifications, its advance (translational motion) was modeled under different rotational motion conditions. Excavation was simulated for 360 s under 12 combinations of conditions-four cutterhead RPMs and three screw conveyor RPMs. The TBM advance was simulated at ~20 mm/min ( Table 4). As there were no particles around the shield skin, the interaction between the shield skin and the particles was not considered in this study.  Eight plate geometries are created and placed in the form of a box to be filled with particles for the ground model. One of the plates located on the +x-axis was designed with a hole of the same diameter as the TBM. Once all plates were arranged, the TBM geometry was initially introduced into the hole by nearly 400 mm. About 400,000 particles were then randomly generated and stacked by gravity. A further 52,000 particles were generated in the chamber and screw conveyor to represent TBM excavation. The top and bottom of the ground boundaries were set as parallel boundaries. Once these works were complete, an additional, fixed, overburden load of 265 kPa was applied. It was also applied continuously to the top plate by the servo control during excavation. As a result, a rectangular parallel-piped ground model of size 1.5 m (width) × 15.6 m (length) × 14.6 m (height) was created over the stabilization time ( Figure 4). randomly generated and stacked by gravity. A further 52,000 particles were generated in the chamber and screw conveyor to represent TBM excavation. The top and bottom of the ground boundaries were set as parallel boundaries. Once these works were complete, an additional, fixed, overburden load of 265 kPa was applied. It was also applied continuously to the top plate by the servo control during excavation. As a result, a rectangular parallel-piped ground model of size 1.5 m (width) × 15.6 m (length) × 14.6 m (height) was created over the stabilization time ( Figure 4). Based on the EPB shield TBM's technical specifications, its advance (translational motion) was modeled under different rotational motion conditions. Excavation was simulated for 360 s under 12 combinations of conditions-four cutterhead RPMs and three screw conveyor RPMs. The TBM advance was simulated at ~20 mm/min ( Table 4). As there were no particles around the shield skin, the interaction between the shield skin and the particles was not considered in this study. Based on the EPB shield TBM's technical specifications, its advance (translational motion) was modeled under different rotational motion conditions. Excavation was simulated for 360 s under 12 combinations of conditions-four cutterhead RPMs and three screw conveyor RPMs. The TBM advance was simulated at~20 mm/min ( Table 4). As there were no particles around the shield skin, the interaction between the shield skin and the particles was not considered in this study.

Analysis of Numerical Simulation Result
The TBM performance indicators of compressive force, pressure, torque, wear, and mass flow rate were evaluated from the simulation.
The compressive force was determined from the sum of the normal contact force magnitudes (Figure 5a). This sum was calculated from the product of the average normal contact force and the contact frequency at data write-out points. The normal contact force was obtained from two bodies in direct contact and acted perpendicular to each body. The pressure was determined from the compressive force and the geometry's surface area, which was composed of structured meshes [31]. Appl. Sci. 2021, 11, 5119 8 of 26

Analysis of Ground Model
3.1.1. Pressure Acting on Plates of the Ground Model Figure 7 shows the results for the lateral pressure, depending on the simulation time and depth. Pressure was estimated from the frictionless plate type geometry with size 17 m (length) × 1 m (height) composing the ground, located in front of the cutterhead face. The vertical pressure, defined as the overburden load, was maintained almost continuously at 265 kPa during the simulation. The estimation of lateral pressure appeared to vary with the plate's depth. Unlike the overburden load, lateral pressure was not continuously controlled during the simulation. Rather, the lateral pressure caused by gravity and the overburden load was measured in the virtual plate geometry as a reaction force, with the measured lateral pressure fluctuating more than the overburden load. Wear is defined as the gradual loss of mass from the movement of a body's surface, arising from its relative motion [32]. Relative wear of the cutterhead, cutting tools, and screw conveyor was analyzed using a relative wear model based on energy. Although there are various kinds of wear, impact wear (defined as the normal force on each element of the geometry) was selected for wear evaluation. Based on the energetic approach, the total normal contact energy was derived as follows: where E n is the normal cumulative contact energy, F n is the normal contact force, V n is the normal relative velocity, and δ t represents each time interval.
The energy approach shows that the disc cutter's impact damage could accumulate as energy and gradually increase over time [33]. The variation of the normal cumulative contact energy of the geometry (from the beginning of excavation) was used to indicate the impact wear.
The muck discharge rate was defined as the mass flow (kg/s), owing to the rotation of the screw conveyor-it was the total mass of particles per unit time measured and recorded by a mass-flow sensor ( Figure 6). The sensor's domain in the x-direction was set at 770 mm and fixed relative to the TBM shield geometry.
The torque on the geometry was the product of the tangential contact force between two bodies and the distance from the center of mass to the contact point ( Figure 5b). The tangential force was the force from the tangential overlap. The total torque was calculated from the torque and moment around each geometry element. The sum of the torque from each element was converted to the total torque on the entire geometry about the center of mass, and was plotted along the selected axis [31].   Figure 7 shows the results for the lateral pressure, depending on the simulation time and depth. Pressure was estimated from the frictionless plate type geometry with size 17 m (length) × 1 m (height) composing the ground, located in front of the cutterhead face. The vertical pressure, defined as the overburden load, was maintained almost continuously at 265 kPa during the simulation. The estimation of lateral pressure appeared to vary with the plate's depth. Unlike the overburden load, lateral pressure was not continuously controlled during the simulation. Rather, the lateral pressure caused by gravity and the overburden load was measured in the virtual plate geometry as a reaction force, with the measured lateral pressure fluctuating more than the overburden load.  Figure 7 shows the results for the lateral pressure, depending on the simulation time and depth. Pressure was estimated from the frictionless plate type geometry with size 17 m (length) × 1 m (height) composing the ground, located in front of the cutterhead face. The vertical pressure, defined as the overburden load, was maintained almost continuously at 265 kPa during the simulation. The estimation of lateral pressure appeared to vary with the plate's depth. Unlike the overburden load, lateral pressure was not continuously controlled during the simulation. Rather, the lateral pressure caused by gravity and the overburden load was measured in the virtual plate geometry as a reaction force, with the measured lateral pressure fluctuating more than the overburden load.  Figure 8 plots the average lateral pressure acting on the plate with respect to the depth. The average pressure tends to increase as the depth increases, although that acting on plate No. 7 was estimated to be lower than that on plate No. 8. The pressure acting on the lowest plate (No. 8) was 1.5 times that on the uppermost plate (No. 1). Until the depth The average pressure tends to increase as the depth increases, although that acting on plate No. 7 was estimated to be lower than that on plate No. 8. The pressure acting on the lowest plate (No. 8) was 1.5 times that on the uppermost plate (No. 1). Until the depth of 6.0 m, the pressure acting on each plate increased 1.1 times for every 1 m increase in the depth of the plate from the top of the TBM cutterhead. Based on the earth pressure theory, depth and lateral pressure are directly related, but as lateral pressure was not controlled by the servo system, in this study they were not related. However, it is quite clear that they exhibit an almost linear proportional relationship.

Pressure Acting on Plates of the Ground Model
Appl. Sci. 2021, 11, 5119 10 of Figure 8. Lateral pressure on the plates of the ground model plotted with respect to location. Figure 9 shows the variation of the soil particles' velocity with respect to cutterhe rotation speed from the rear view of the TBM advance at a specific simulation time. At rev/min, only particles around the TBM have a relatively high velocity as compared those elsewhere in the ground. In contrast, increasing the cutterhead rotation to rev/min greatly increased the velocity of particles around the TBM. The propagation the disturbance zone around the tunnel face appeared to increase throughout the grou model, as the areas colored green and red were larger.  Figure 9 shows the variation of the soil particles' velocity with respect to cutterhead rotation speed from the rear view of the TBM advance at a specific simulation time. At 1.0 rev/min, only particles around the TBM have a relatively high velocity as compared to those elsewhere in the ground. In contrast, increasing the cutterhead rotation to 3.0 rev/min greatly increased the velocity of particles around the TBM. The propagation of the disturbance zone around the tunnel face appeared to increase throughout the ground model, as the areas colored green and red were larger.  Figure 10 shows the soil particles' velocity in the chamber with respect to the cutterhead's rotational speed at a specific simulation time. Their overall velocity increased as the cutterhead rotated more quickly. A particle further from the chamber's center will tend to move more quickly than one centrally located. Therefore, particles at the chamber's center showed lower fluidity than those further from the center. A lack of opening in the central area of this TBM's cutterhead appeared to be the cause of the low flow of particles in the center of the chamber. This tendency is frequently observed in the fieldthe center of a TBM's cutterhead is generally the first area to become blocked or clogged, rather than the face or gauge area [34]. Figure 11 shows the variation of the total number of particles in the chamber during the simulation. As the cutterhead's rotational speed increased, there were overall more particles in the chamber. The number of particles in the chamber greatly increased in the first 90 s. The rate of increase tended to decrease as the screw's rotational speed increased. The simulation started with some particles already in the chamber and screw conveyor. A sudden inflow of particles from the cutterhead to the opening area was expected to occur immediately upon starting the simulation. This inflow affected the chamber system regardless of the chosen screw rotational speed. Table 5 summarizes the standard deviation (SD) of the number of particles in the chamber with respect to simulation time under different operating conditions. Case 8 had   Figure 10 shows the soil particles' velocity in the chamber with respect to the cutterhead's rotational speed at a specific simulation time. Their overall velocity increased as the cutterhead rotated more quickly. A particle further from the chamber's center will tend to move more quickly than one centrally located. Therefore, particles at the chamber's center showed lower fluidity than those further from the center. A lack of opening in the central area of this TBM's cutterhead appeared to be the cause of the low flow of particles in the center of the chamber. This tendency is frequently observed in the field-the center of a TBM's cutterhead is generally the first area to become blocked or clogged, rather than the face or gauge area [34]. Figure 11 shows the variation of the total number of particles in the chamber during the simulation. As the cutterhead's rotational speed increased, there were overall more particles in the chamber. The number of particles in the chamber greatly increased in the first 90 s. The rate of increase tended to decrease as the screw's rotational speed increased. The simulation started with some particles already in the chamber and screw conveyor. A sudden inflow of particles from the cutterhead to the opening area was expected to occur immediately upon starting the simulation. This inflow affected the chamber system regardless of the chosen screw rotational speed. the lowest value (SD = 18) with 2.0 rev/min cutterhead speed and 9.0 rev/min screw-conveyor speed, followed by Cases 7 and 10. Apart from the initial 90 s of simulation, the SD of the number of particles was highest at 6.0 rev/min screw-conveyor speed.  Figure 12 shows the pressure acting on the chamber surface during simulation. As explained in Figure 11, the expected sudden inflow of particles from the cutterhead to the opening area would affect the chamber system regardless of the conditions of screw rotation speed in the first 90 s. Table 6 summarizes the average chamber pressures and the SD of chamber pressures, with respect to operating conditions. As the cutterhead speed increased, the chamber pressure increased at the same screw-conveyor speed. The chamber pressure increased as the screw-conveyor speed decreased to 6.0 rev/min at the same cutterhead speed. Further- Figure 11. Variation of the total number of particles in the chamber during simulation. Table 5 summarizes the standard deviation (SD) of the number of particles in the chamber with respect to simulation time under different operating conditions. Case 8 had the lowest value (SD = 18) with 2.0 rev/min cutterhead speed and 9.0 rev/min screwconveyor speed, followed by Cases 7 and 10. Apart from the initial 90 s of simulation, the SD of the number of particles was highest at 6.0 rev/min screw-conveyor speed. Table 5. Standard deviation of the number of particles in the chamber depending on operating conditions.

Case
Operating  Figure 12 shows the pressure acting on the chamber surface during simulation. As explained in Figure 11, the expected sudden inflow of particles from the cutterhead to the opening area would affect the chamber system regardless of the conditions of screw rotation speed in the first 90 s.  Case 4 shows the lowest value of the SD of chamber pressure (0.72) under the conditions of 1.5 rev/min cutterhead rotation and 12.0 rev/min screw conveyor rotation, followed by Cases 3 and 7. The pressure acting on the chamber surface could be estimated differently, depending on the total mass of particles in the chamber and each particle's mechanical energy contacting with the chamber surface. Therefore, it was expected that the trends of the two plots in Figures 11 and 12 would be slightly different.   Table 6 summarizes the average chamber pressures and the SD of chamber pressures, with respect to operating conditions. As the cutterhead speed increased, the chamber pressure increased at the same screw-conveyor speed. The chamber pressure increased as the screw-conveyor speed decreased to 6.0 rev/min at the same cutterhead speed. Furthermore, the chamber pressure increased by 6.0% to 16.0%, as the screw-conveyor speed decreased to 6.0 rev/min and the cutterhead speed increased to 3.0 rev/min. For Case 4 to Case 6, the pressure increased by~6.0%; for Case 10 to Case 12, it increased bỹ 16.0%. Comparing cases 3 and 4, the average chamber pressures were almost constant with the cutterhead speed increasing from 1.0 to 1.5 rev/min and screw-conveyor speed from 6.0 to 12.0 rev/min. In cases 9 and 10, the average chamber pressure increased bỹ 3.8%, implying that the screw-conveyor rotation speed should be >12.0 rev/min to balance chamber pressure. Case 4 shows the lowest value of the SD of chamber pressure (0.72) under the conditions of 1.5 rev/min cutterhead rotation and 12.0 rev/min screw conveyor rotation, followed by Cases 3 and 7. The pressure acting on the chamber surface could be estimated differently, depending on the total mass of particles in the chamber and each particle's mechanical energy contacting with the chamber surface. Therefore, it was expected that the trends of the two plots in Figures 11 and 12 would be slightly different.

Compressive Force Acting at the Cutterhead's Face
Here, the TBM compressive force was calculated from the sum of normal forces between the particle and TBM cutterhead, and it was, therefore, parallel to the advancedirection axis. Figure 13a plots the compressive force acting on the TBM cutterhead as the TBM advanced. The plot of the average compressive force in Figure 13b clearly shows that as the cutterhead RPM increased, the overall compressive force also increased. The DEM simulation showed that the calculated compressive force on the cutterhead and disc cutter could be considered to be a component of the required thrust force. Based on this consideration, the average compressive force on the cutterhead was estimated to be 31-37% of the actual TBM's maximum thrust. direction axis. Figure 13a plots the compressive force acting on the TBM cutterhead as the TBM advanced. The plot of the average compressive force in Figure 13b clearly shows that as the cutterhead RPM increased, the overall compressive force also increased. The DEM simulation showed that the calculated compressive force on the cutterhead and disc cutter could be considered to be a component of the required thrust force. Based on this consideration, the average compressive force on the cutterhead was estimated to be 31%-37% of the actual TBM's maximum thrust.  Figure 14a presents the torque on the TBM cutterhead as the TBM advanced. The torque increased as the cutterhead RPM increased. Figure 14b clearly shows that the overall torque increased to ~7000 kN•m as the cutterhead's rotation increased to 3.0 rev/min. The average torque on the cutterhead was estimated to be 50-70% of its actual maximum, under different conditions of the cutterhead RPM.  Figure 14a presents the torque on the TBM cutterhead as the TBM advanced. The torque increased as the cutterhead RPM increased. Figure 14b clearly shows that the overall torque increased to~7000 kN·m as the cutterhead's rotation increased to 3.0 rev/min. The average torque on the cutterhead was estimated to be 50-70% of its actual maximum, under different conditions of the cutterhead RPM. Power (kW) could be calculated as follows from the torque (kN•m) and RPM:

Torque and Driving Power of the TBM Cutterhead
Power = Torque × RPM 9.5488 (2) Figure 15a presents the driving power required to drive the cutterhead as the TBM advanced. Figure 15b shows that the calculated power increased dramatically as the cutterhead rotation increased to 3.0 rev/min; however, at 3.0 rev/min, the required power exceeded the actual maximum capacity of TBM. Power (kW) could be calculated as follows from the torque (kN·m) and RPM: Figure 15a presents the driving power required to drive the cutterhead as the TBM advanced. Figure 15b shows that the calculated power increased dramatically as the cutterhead rotation increased to 3.0 rev/min; however, at 3.0 rev/min, the required power exceeded the actual maximum capacity of TBM.

Power =
Torque × RPM 9.5488 (2) Figure 15a presents the driving power required to drive the cutterhead as the TBM advanced. Figure 15b shows that the calculated power increased dramatically as the cutterhead rotation increased to 3.0 rev/min; however, at 3.0 rev/min, the required power exceeded the actual maximum capacity of TBM.

Correlation between Compressive Force and Cutterhead Torque
The calculated compressive force on the cutterhead could be evaluated as part of the TBM's total thrust. Figure 16a shows that the compressive force and torque on the cutterhead were related by a regression equation with an exponential function. The cutterhead torque and TBM thrust would thus be closely related as the TBM advanced. Figure 16b shows the ratio of compressive force and torque as the TBM advanced. The ratio gradually decreased as the cutterhead's rotation increased from 1.0 to 3.0 rev/min. Table 7 shows the calculated average and the SD of the ratio for different operating conditions. Case 2 shows the highest average ratio of about 3.72, and Case 11 shows the

Correlation between Compressive Force and Cutterhead Torque
The calculated compressive force on the cutterhead could be evaluated as part of the TBM's total thrust. Figure 16a shows that the compressive force and torque on the cutterhead were related by a regression equation with an exponential function. The cutterhead torque and TBM thrust would thus be closely related as the TBM advanced. Figure 16b shows the ratio of compressive force and torque as the TBM advanced. The ratio gradually decreased as the cutterhead's rotation increased from 1.0 to 3.0 rev/min. lowest SD of the ratio. Overall, the average and SD of the ratio decreased as the cutterhead RPM increased. The ratio of thrust and torque was a critical indicator for evaluating a TBM's excavation performance. If the TBM met adverse geological conditions, the ratio increased abruptly. As successful excavation generally showed low fluctuation in the ratio as the TBM advanced, Case 11 had the most appropriate conditions in terms of operational management. Table 7. Average and standard deviation of the ratio of compressive force to torque for different operating conditions.

Compressive Force Acting on the Disc Cutter
The cutterhead face could be divided into three zones-center, face, and gauge area. The compressive force acting on a disc cutter depended on its position on the cutterhead. Figure 17 shows the estimated compressive forces acting on all disc cutters to be <311 kN, the amount allowable on a 19-inch disc cutter. Table 8 shows that the higher compressive forces acted on a few disc cutters (19,30,43,50, and 51) close to the center of the cutterhead, while the lower compressive forces acted on certain other cutters (7, 24, 31, 32, and  Table 7 shows the calculated average and the SD of the ratio for different operating conditions. Case 2 shows the highest average ratio of about 3.72, and Case 11 shows the lowest SD of the ratio. Overall, the average and SD of the ratio decreased as the cutterhead RPM increased. The ratio of thrust and torque was a critical indicator for evaluating a TBM's excavation performance. If the TBM met adverse geological conditions, the ratio increased abruptly. As successful excavation generally showed low fluctuation in the ratio as the TBM advanced, Case 11 had the most appropriate conditions in terms of operational management. Table 7. Average and standard deviation of the ratio of compressive force to torque for different operating conditions.

Case
Operating The cutterhead face could be divided into three zones-center, face, and gauge area. The compressive force acting on a disc cutter depended on its position on the cutterhead. Figure 17 shows the estimated compressive forces acting on all disc cutters to be <311 kN, the amount allowable on a 19-inch disc cutter. Table 8 shows that the higher compressive forces acted on a few disc cutters (19,30,43,50, and 51) close to the center of the cutterhead, while the lower compressive forces acted on certain other cutters (7, 24, 31, 32, and 49). The five central disc cutters, thus, appeared vulnerable to impact damage during excavation, a trend confirmed by analysis of TBM wear during excavation in the field [3]. To protect the disc cutters from impact damage, those at the center were generally twin-or double-types, with a higher maximum load capacity instead of single disc cutter types.
Appl. Sci. 2021, 11, 5119 18 of 26 49). The five central disc cutters, thus, appeared vulnerable to impact damage during excavation, a trend confirmed by analysis of TBM wear during excavation in the field [3]. To protect the disc cutters from impact damage, those at the center were generally twin-or double-types, with a higher maximum load capacity instead of single disc cutter types. Similar to the results for the cutterhead (Figure 13b), the average compressive force on all disc cutters increased as the cutterhead RPM increased, and was estimated to be 16-17% of the compressive force acting on the cutterhead. Furthermore, the total compressive force acting on the disc cutters and cutterhead were estimated to be 37-43% of the actual TBM's maximum thrust ( Figure 18). Figure 19 shows the dependence of the average compressive force acting on each disc cutter on the cutter's distance from the cutterhead center. The force tended to decrease as the distance increased to ~1750 mm. Moving farther outward from 1750 mm to the gauge area, it increased, remained constant, and then decreased. As there was no opening at the center of the cutterhead, a localized high pressure could emerge that decreased the fluidity of the surrounding soil particles [34]. The high pressure could also greatly affect any centrally placed disc cutter. As the radius from the center increased, the radial opening ratio and the fluidity of the soil particle also increased. The pressure on the disc cutter therefore decreased. However, as the radius from the center increased to 2240 mm, close to the connection between the main spoke and two auxiliary spokes, this connection acted as a barrier to decrease the flow rate. Hence, the compressive force acting on the disc cutter increased abruptly, and the fitted trend line of force with respect to the radius from the center also changed. Figure 17. Average compressive force on each numbered disc cutter. Figure 17. Average compressive force on each numbered disc cutter. Similar to the results for the cutterhead (Figure 13b), the average compressive force on all disc cutters increased as the cutterhead RPM increased, and was estimated to be 16-17% of the compressive force acting on the cutterhead. Furthermore, the total compressive force acting on the disc cutters and cutterhead were estimated to be 37-43% of the actual TBM's maximum thrust (Figure 18).
Appl. Sci. 2021, 11, 5119 19 of 26 Figure 18. Average compressive force acting on all disc cutters with respect to cutterhead RPM. Figure 18. Average compressive force acting on all disc cutters with respect to cutterhead RPM. Figure 19 shows the dependence of the average compressive force acting on each disc cutter on the cutter's distance from the cutterhead center. The force tended to decrease as the distance increased to~1750 mm. Moving farther outward from 1750 mm to the gauge area, it increased, remained constant, and then decreased. As there was no opening at the center of the cutterhead, a localized high pressure could emerge that decreased the fluidity of the surrounding soil particles [34]. The high pressure could also greatly affect any centrally placed disc cutter. As the radius from the center increased, the radial opening ratio and the fluidity of the soil particle also increased. The pressure on the disc cutter therefore decreased. However, as the radius from the center increased to 2240 mm, close to the connection between the main spoke and two auxiliary spokes, this connection acted as a barrier to decrease the flow rate. Hence, the compressive force acting on the disc cutter increased abruptly, and the fitted trend line of force with respect to the radius from the center also changed. Figure 18. Average compressive force acting on all disc cutters with respect to cutterhead RPM. Figure 19. Average compressive force on disc cutters depending on their location. Figure 19. Average compressive force on disc cutters depending on their location.

Impact Wear Estimated from Disc Cutter
Here, the impact wear was calculated from the sum of normal contact forces between particles and each TBM disc cutter. Figure 20a presents the gradual accumulation of average total normal contact energy in all disc cutters, as the TBM advanced during the simulation.

Impact Wear Estimated from Disc Cutter
Here, the impact wear was calculated from the sum of normal contact forces between particles and each TBM disc cutter. Figure 20a presents the gradual accumulation of average total normal contact energy in all disc cutters, as the TBM advanced during the simulation. Figure 20b presents the variation of the average total normal accumulative contact energy in all disc cutters. It clearly increased as the cutterhead RPM rose. Figure 21 and Table 9 show the normal cumulative contact energy in each disc cutter after excavation. The gauge disc cutters (shown in red) at the transition zone of the cutterhead have higher energy, while those centrally placed (shown in blue) have lower energy compared with the rest. Overall, the figure shows that a disc cutter's wear depended significantly on its position on the cutterhead and the layout pattern. Therefore, a gauge disc cutter in the transition zone of the cutterhead was expected to be replaced frequently during the excavation [3][4][5][6][7][8]. Figure 22 shows the results of each disc cutter's normal cumulative contact energy with respect to its radial area on the cutterhead. There was no clear correlation between radial distance and normal cumulative contact energy until the radial distance increased to 2000 mm. Above this distance, the disc cutter's normal cumulative contact energy increased significantly as its radial distance increased.   Figure 20b presents the variation of the average total normal accumulative contact energy in all disc cutters. It clearly increased as the cutterhead RPM rose. Figure 21 and Table 9 show the normal cumulative contact energy in each disc cutter after excavation. The gauge disc cutters (shown in red) at the transition zone of the cutterhead have higher energy, while those centrally placed (shown in blue) have lower energy compared with the rest. Overall, the figure shows that a disc cutter's wear depended significantly on its position on the cutterhead and the layout pattern. Therefore, a gauge disc cutter in the transition zone of the cutterhead was expected to be replaced frequently during the excavation [3][4][5][6][7][8].    Figure 22 shows the results of each disc cutter's normal cumulative contact energy with respect to its radial area on the cutterhead. There was no clear correlation between radial distance and normal cumulative contact energy until the radial distance increased to 2000 mm. Above this distance, the disc cutter's normal cumulative contact energy increased significantly as its radial distance increased.   The muck discharge rate was evaluated from the specific domain's particle mass flow rate, as presented in Figure 6. Figure 23a shows the decreasing trend of the particle mass flow rate with the screw conveyor's decreasing RPM as the TBM advanced. Figure 23b clearly indicates that the average mass flow rate increased linearly as the screw conveyor RPM increased. For a given screw conveyor RPM, the mass flow rate increased 1.2-1.4 times, as the cutterhead RPM changed from 1.0 to 3.0. Overall, the cutterhead's operating condition affected the performance of the screw conveyor.

Driving Power of the Screw Conveyor
The driving power was calculated from the estimated torque and rotation speed. Figure 24a plots the driving power of the screw conveyor as the TBM advanced. During excavation, the calculated driving power was mostly lower than the maximum power of an actual screw conveyor, except for very few points at a rotation speed of 12 rev/min. Figure 24b shows that the average driving power of the screw conveyor increased as its RPM rose. For a given screw conveyor operating condition, the driving power of the screw conveyor increased by 1.5-1.9 times as the cutterhead rotation increased from 1.0 to 3.0 rev/min. This gives further evidence of the cutterhead's operating conditions affecting the performance of the screw conveyor.
The muck discharge rate was evaluated from the specific domain's particle mass flow rate, as presented in Figure 6. Figure 23a shows the decreasing trend of the particle mass flow rate with the screw conveyor's decreasing RPM as the TBM advanced. Figure 23b clearly indicates that the average mass flow rate increased linearly as the screw conveyor RPM increased. For a given screw conveyor RPM, the mass flow rate increased 1.2-1.4 times, as the cutterhead RPM changed from 1.0 to 3.0. Overall, the cutterhead's operating condition affected the performance of the screw conveyor.

Driving Power of the Screw Conveyor
The driving power was calculated from the estimated torque and rotation speed. Figure 24a plots the driving power of the screw conveyor as the TBM advanced. During excavation, the calculated driving power was mostly lower than the maximum power of an actual screw conveyor, except for very few points at a rotation speed of 12 rev/min. Figure  24b shows that the average driving power of the screw conveyor increased as its RPM rose. For a given screw conveyor operating condition, the driving power of the screw conveyor increased by 1.5-1.9 times as the cutterhead rotation increased from 1.0 to 3.0

Determination of Optimal Operating Condition
The simulation results could then be comprehensively evaluated to determine the TBM's optimal operating conditions at the fixed advance rate with 20 mm/min. As the EPB shield TBM must balance the pressure conditions at the tunnel face during excavation, the number of particles in the chamber and the surface pressure in the chamber were selected as evaluation indices. The SD of the ratio of compressive force and torque was also chosen as an evaluation index to reflect the basic operational management of TBM. Based on these three indices, Table 10 ranks the top five operating conditions.

Determination of Optimal Operating Condition
The simulation results could then be comprehensively evaluated to determine the TBM's optimal operating conditions at the fixed advance rate with 20 mm/min. As the EPB shield TBM must balance the pressure conditions at the tunnel face during excavation, the number of particles in the chamber and the surface pressure in the chamber were selected as evaluation indices. The SD of the ratio of compressive force and torque was also chosen as an evaluation index to reflect the basic operational management of TBM. Based on these three indices, Table 10 ranks the top five operating conditions.  In terms of pressure balance in the chamber, Case 7 with cutterhead rotation at 2.0 rev/min and screw conveyor rotation at 12.0 rev/min was best (average, 2.5). However, the driving power of the screw conveyor exceeded its maximum capacity of 75 kW at very few points when the screw conveyor rotated at 12.0 rev/min. The screw conveyor was thus likely to break as it exceeded its maximum capacity. Therefore, it was necessary to slow the screw below 12.0 rev/min under the given simulation conditions. Without rotating the screw conveyor at 12.0 rev/min, Case 5 (1.5 rpm cutterhead rotation and 9.0 rev/min screw conveyor rotation) was best in terms of pressure balance in the chamber.
Regarding fluctuation of the ratio of the compressive force to torque, Case 11 was best. However, the driving power of the cutterhead exceeded the maximum capacity (1440 kW) at 3.0 rev/min cutterhead rotation. Considering operational management of the TBM, Case 5 (1.5 rev/min cutterhead rotation and 9.0 rev/min screw conveyor rotation) was also best at the fixed advance rate with 20 mm/min. As a result, it was concluded that the operating condition with 1.5 rev/min cutterhead rotation and 9.0 rev/min screw conveyor rotation was the most optimal for TBM operation at the fixed advance rate with 20 mm/min in this TBM simulation.

Conclusions
This study used discrete-element analysis to predict the excavation performance of a 7.73 m diameter EPB shield TBM under different operating conditions. The model of the machine comprises six parts-cutterhead, disc cutters, cutter bits, chamber, shield, and screw conveyor. It was simulated to advance through a ground model of size 1.5 m (width) × 15.6 m (length) × 14.6 m (height), comprising 400,000 spherical particles of diameter 100 mm, randomly generated and stacked by gravity, after placement of the TBM geometry. Together, the two models simulated excavation for 360 s at four cutterhead RPMs and three screw conveyor RPMs to determine the best operation conditions. The simulation showed that the lateral pressure acting on the side plate of the ground model increased almost linearly with the increasing depth of the plate. The ground also became less stable as the cutterhead RPM increased during excavation. The propagation of ground disturbance on the periphery of the cutterhead steadily increased as the cutterhead RPM increased.
The particles' overall velocity in the chamber increased as the cutterhead RPM increased, with particles further from the chamber's center tending to move more quickly. The low flow of particles was expected in the center of the chamber. The number of particles in the chamber and the surface pressure on the chamber varied with the cutterhead RPM and screw conveyor RPM during the simulation.
The ratio of compressive force and torque and the fluctuation of this ratio decreased as the cutterhead RPM rose. The overall compressive force acting on all disc cutters also increased linearly with the cutterhead RPM. The position of a disc cutter on the TBM cutterhead had a strong influence, with higher compressive forces experienced by the cutters closer to the center. The impact wear on a disc cutter determined from the variation of normal accumulative contact energy depended on the cutterhead RPM and the cutter's position on the TBM. It clearly increased as the cutterhead RPM rose. Particularly, gauge disc cutters at the transition zone of the cutterhead showed more wear than those elsewhere, because of the large turning radius.
The muck discharge rate was measured by a mass flow rate sensor. The average mass flow rate gradually increased as the cutterhead RPM increased. Additionally, the average driving power of the screw conveyor gradually increased as the cutterhead RPM increased. Overall, the cutterhead operating conditions affected the performance of the screw conveyor.
Finally, comprehensive evaluation of TBM excavation under different simulation conditions led to the determination of the optimal operation conditions based on pressure balance and operational management of the TBM.
Further work should focus on two aspects-soil conditioning and load (pressure) control; the excavated material properties generally change with soil conditioning. If the ground is modeled by considering the soil conditioning work and calibration work with relevant soil samples, it might improve the prediction of the excavation before the TBM starts tunneling. This study simulated excavation using only DEM, except for the overburden load control system. Therefore, it is impossible to apply a high-class control system for real-time control of the rate of advance and cutterhead RPM, depending on the estimated torque, compressive force, and the chamber pressure. As additional pressure is determined by how full the chamber is, an effective control system should be considered in future studies. It is also difficult to consider the rotation of disc cutters during excavation. To overcome this limitation, coupling analysis will be attempted based on multi-body dynamics and DEM. Further studies should also compare simulation results with field data.