Discrete/Finite Element Modelling of Rock Cutting with a TBM Disc Cutter

This paper presents advanced computer simulation of rock cutting process typical for excavation works in civil engineering. Theoretical formulation of the hybrid discrete/finite element model has been presented. The discrete and finite element methods have been used in different subdomains of a rock sample according to expected material behaviour, the part which is fractured and damaged during cutting is discretized with the discrete elements while the other part is treated as a continuous body and it is modelled using the finite element method. In this way, an optimum model is created, enabling a proper representation of the physical phenomena during cutting and efficient numerical computation. The model has been applied to simulation of the laboratory test of rock cutting with a single TBM (tunnel boring machine) disc cutter. The micromechanical parameters have been determined using the dimensionless relationships between micro- and macroscopic parameters. A number of numerical simulations of the LCM test in the unrelieved and relieved cutting modes have been performed. Numerical results have been compared with available data from in-situ measurements in a real TBM as well as with the theoretical predictions showing quite a good agreement. The numerical model has provided a new insight into the cutting mechanism enabling us to investigate the stress and pressure distribution at the tool–rock interaction. Sensitivity analysis of rock cutting performed for different parameters including disc geometry, cutting velocity, disc penetration and spacing has shown that the presented numerical model is a suitable tool for the design and optimization of rock cutting process.


Your article is published under the Creative
Commons Attribution license which allows users to read, copy, distribute and make derivative works, as long as the author of the original work is cited. You may selfarchive this article on your own website, an institutional repository or funder's repository and make it publicly available immediately.

Introduction
Variety of rock cutting technologies is used in civil or mining engineering. Rock cutting consists in fracturing and disintegration of a rock using different methods and different machines. Tunnel boring machines (TBMs), shown in Fig. 1a, are used to perform rock cutting in excavation of tunnels. In excavation with a TBM, a rock is cut by means of cutter discs (Fig. 1b) installed on a rotating cutter head, which is pressed against the tunnel face.
TBMs are widely used in various tunnelling projects in civil engineering (road and railway tunnels), mining industry (tunnels for access to underground excavations, conveyance of ore and waste, drainage, exploration, water supply and diversion) and other geotechnical engineering applications. The use of TBMs is continuously growing mainly due to their efficiency. Nevertheless, there is still a need of improvement of TBM performance depending mainly on rock properties, operational parameters (machine trust, penetration and rate of advance) and design of the cutter head design, including design of the disc cutters (Roby et al. 2008) and design of their layout (Huo et al. 2010).
Historically, design of rock cutting tools for the excavation machinery has been based on a combination of the experience of engineers and real size laboratory tests, resulting many times in an inefficient process, and involving high costs for the excavation companies. Different empirical models have been developed for the estimation of the principal parameters involved (Nilsen and Ozdemir 1993;Rostami and Ozdemir 1996). These models are useful in certain cases; nevertheless, their use is restricted by the availability of historical data and range of rock material properties (Ramezanzadeh et al. 2004).
Application of numerical analysis improves the design methodology and allows to obtain more efficiently optimized designs of rock cutting tools and machines. Numerical methods can be used to optimize the TBM cutter layout (Sun et al. 2015) or simulate interaction between TBM components and rock mass (Zhao et al. 2014). The present work is focused on simulation of rock cutting with a single TBM disc cutter. Such an analysis can help the designer to understand better rock cutting mechanisms, to detect the reasons of the cutting tool wear and failure, and finally to improve the cutter design and determine optimum operational parameters. The aim of the analysis is to substitute or at least reduce number of laboratory tests used in the design process of rock cutting tools (Nilsen and Ozdemir 1993;Rostami and Ozdemir 1993), and finally speed-up the design process and reduce its costs.
Simulation of excavation processes, in general, and of rock cutting process, in particular, is not an easy task. Continuum-based simulation techniques, such as the finite element method (FEM), encounter serious difficulties in modelling of fracture and fragmentation of a rock material occurring in an excavation process (Jonak and Podgórski 2001;Yu and Khair 2007;Shenghua 2004;Loui and Karanam 2012). The discrete element method (DEM) employing a discrete material model offers a more realistic way to simulate discontinuous phenomena. The DEM has been successfully applied to simulation of different rock cutting processes by Stavropoulou (2006), Rojek et al. (2011), Su and Akcin (2011), Huang et al. (2013), Labra et al. (2008), Wyk et al. (2014).
The known disadvantage of the DEM is that it usually requires the use of a large number of elements, which leads to long computation times. This paper put forwards the idea allowing to reduce the computation cost of the DEM simulation of rock cutting by coupling the DEM with the FEM and using them in different subdomains of the cut material according to expected material behaviour. In rock cutting problems, a sufficiently large specimen must be taken in order to avoid artificial boundary effects. In such specimens, a large part of the rock material is not damaged and can be treated as a continuous material. A continuous material deformation is usually modelled efficiently using the FEM. Therefore, an optimum model of rock cutting can be obtained combining the DEM with the FEM in such a way that discrete elements are used only in a portion of the analysed domain where material fracture occurs, while outside the DEM subdomain finite elements are used Labra et al. 2008). The numerical method proposed in this paper is based on the formulation presented by Rojek and Oñate (2007). Preliminary results for a 2D rock cutting problems obtained with the DEM/ FEM coupled method have been presented by Labra et al. (2008). In the present paper, this algorithm is applied to 3D simulation of the rock cutting problem. This work is original with respect to other research. As far as we know, DEM/FEM coupling has not been used for rock cutting simulations by other authors. Other studies of rock cutting, cf. Stavropoulou (2006), Rojek et al. (2011), Su and Akcin (2011), Huang et al. (2013), Wyk et al. (2014, have employed pure DEM models. Numerical simulations presented here have been performed for the linear cutting machine (LCM) test. The LCM test is one of the most representative testing procedures for prediction of TBM disc cutters (Nilsen and Ozdemir 1993;Rostami and Ozdemir 1993), and it is used extensively in the design of rock cutting tools. Both unrelieved and relieved cutting modes have been studied. The simulation results have been compared with available experimental data and the predictions of the empirical model developed by Rostami (1997). The use of numerical models and simulation has allowed us to study some aspects of rock cutting such as the stress field in rock and disc or pressure in the contact zone which are difficult to study using experimental techniques. Direct measurements of the load distribution in the contact between disc cutter and rock performed by Rostami (2013) have proved that the problem of the load and stress concentration in the contact zone is far from being fully understood and requires further studies. The numerical model presented in this work can be a useful tool to verify theoretical models for rock cutting with disc cutters. Rock fragmentation induced by TBM disc cutters has also been studied by other authors using DEM models (Gong et al. 2006;Moon et al. 2006), but in those works the authors employed 2D models so the simulations consisted in analysis of indentation under plane stress conditions. The present work is the first one employing a full 3D model reproducing real cutting conditions for a TBM disc cutter.

Linear Cutting Test
Performance of a TBM depends on the rock breakage mechanism induced by disc cutters (Yagiz et al. 2010). Rock failure mechanism during cutting with a disc cutter is illustrated in Fig. 2. A crushed zone develops beneath the cutter as it is forced into the rock. As stresses continue to build up in the crushed zone, radial cracks begin to form and propagate into the rock. When one or more of these cracks meet those developed from adjacent cut, chips are released (Balci and Bilgin 2007). The most important factors influencing cutting efficiency are the disc geometry (diameter and tip width), spacing between disc cutters, and penetration.
Tool-rock interaction is characterized by the reaction (cutting) forces. The resultant cutting forces over a single disc cutter are decomposed into three components: the normal force, rolling force and side force, as depicted in Fig. 3. Cutting forces have great influence on the efficiency of the cutting process. They are used to estimate global forces (thrust force and torque) over the TBM (Ramezanzadeh et al. 2004;Farrokh et al. 2012). Cutting forces can be estimated experimentally or theoretically.
The linear cutting machine (LCM) test ( Fig. 4) is one of the most common experimental procedures to predict performance of a single cutting tool (EMI 2016;Nilsen and Ozdemir 1993;Rostami and Ozdemir 1993). It is a full-scale test which was originally developed by the Colorado School of Mines (CSM). The LCM test provides a direct measure of cutting forces and rock cuttability under pre-defined process parameters such as the cutter spacing, cutter penetration, cutter thrust and cutting speed.
The LCM test is a suitable test to validate a numerical model of rock cutting. With the simulations of the LCM test for a single disc cutter, it is possible to perform an analysis of the effect of different parameters such as velocity, penetration rate or tool geometry, on the resultant forces, as well as their influence on the rock fracture.

Cutting Forces Prediction Models
Cutting forces estimation is based on the correlation of different parameters, such as the disc cutter geometry, spacing, penetration rate, disc rolling velocity and rock material properties. Early prediction models for single V-shape disc cutters were proposed by Roxborough and Phillips (1975), Sanio (1985) and Sato and Itakura (1991). Colorado School of Mines developed a model for the CCSshape cutters (Rostami and Ozdemir 1993;Rostami 1997). This model will be used here for the comparison with the simulation results for the LCM test. The first version of the Colorado School of Mines model was developed by Ozdemir (1977), and later was updated by Rostami and Ozdemir (1993), Rostami (1997). The CSM model estimates the cutting forces considering a given penetration, rock mass properties, cutter geometry and cutting conditions. The model is based on a large database of full-scale LCM tests and does not consider rock mass conditions such as fractures or joints.
The model proposes a pressure distribution P in the crushed zone as where w is a constant for the pressure distribution function (typically varying between 0.2 for V-shape and very sharp cutters and -0.2 for wider tip cutters), /-the angle of contact between the rock and the cutter, defined as and P o is the base pressure in the crushed zone, established from regression analysis of several tests, and estimated from rock strength and cutting geometry: where C is a dimensionless constant (usually C ¼ 2:12), and s the spacing between cutters.
The total resulting cutting force F T can be obtained by integrating the pressure over the contact area ( Fig. 5), as where T is the cutter tip width and R the cutter radius.
To recover the normal and rolling forces, the cutting coefficient CC (also called rolling coefficient) is used, which is the ratio of both forces defined by the angle b as Assuming a uniform distribution of the pressure in the contact area, the model proposes a geometrical definition of b as the middle point of the contact area, i.e. Scheme of forces acting on a disc cutter and incidence angles so finally, the normal and rolling forces are estimated projecting the total force in each direction as This model has been used for the estimation of the TBM cutterhead performance in many tunnelling projects with a high degree of success (Rostami and Ozdemir 1996;Rostami 2008;Exadaktylos et al. 2008).

Specific Energy
Rock cutting efficiency is commonly assessed using the specific energy as its criterion. The specific energy SE is defined as the amount of energy required to excavate a unit volume of rock. SE for a single disc cutter is defined as where F r is the rolling force, L the cutting distance, and V the cutting volume. Assuming optimum cutting performance, the cutting volume V can be expressed in terms of the penetration depth p and spacing between disc cutters s as V ¼ p L s, which allows us to rewrite Eq. (9) as follows The ratio of the spacing and penetration (s / p) is a suitable variable to investigate cutting efficiency with TBMs (Sato and Itakura 1991;Rostami 1997Rostami , 2008. The specific energy can be used to predict TBM performance. The procedure is described by (Bilgin et al. 1999(Bilgin et al. , 2014. Given the specific energy SE (in kWh/m 3 ), the net production rate NPR (in m 3 /h) of the excavating machine can be calculated from the following equation: where k is the energy transfer ratio from the cutting head to the tunnel face, it is usually taken as equal to 0.8 for TBM, cf. Bilgin et al. (2014), and P (in kW) is the power used to excavate rock expressed in terms of the TBM torque T (in kNm) and the rotational velocity of the cutterhead N (in revolutions/sec) as follows: 4 Formulation of the Model

Basic Assumptions
A numerical model allowing us to simulate the LCM test has been developed. A system consisting of a tool and rock sample is considered in the model. The tool is considered as a rigid body. The rock is modelled using a hybrid discrete-finite element method approach in which a part of the rock near the tool is modelled using the discrete element method (DEM) and the other part is considered using the finite element method (FEM). The DEM and FEM subdomains are coupled using special kinematic constraints. The tool-rock interaction is modelled assuming the Coulomb friction model. Rock fracture during cutting is assumed to be localized in the DEM subdomain, and the FEM subdomain is assumed to be continuous and linearly elastic.

Discrete Element Method Formulation
The discrete element model assumes that material can be represented by an assembly of distinct particles or bodies interacting among themselves. Generally, discrete elements can have arbitrary shape. In this work, the formulation employing spherical rigid particles is used. Basic formulation of the particle-based discrete element method was first proposed by Cundall and Strack (1979). A similar formulation of the DEM has been implemented in the explicit dynamic finite element method programme by Rojek and Oñate (2004), , starting the development of the DEM/FEM code DEMpack (CIMNE 2010). The translational and rotational motion of the rigid spherical elements (particles) is governed by the standard equations of rigid body dynamics. For the i-th element, we have where € u i is the position vector of the element centroid in a fixed (inertial) coordinate frame, x i is the angular velocity, m i is the element mass, J i is the element moment of inertia, F i is the vector of resultant forces, and T i is the vector of resultant moments about the element central axes. Equations of motion (13)-(14) are integrated in time using the central difference scheme.
Vectors F i and T i are sums of all forces and moments applied to the i-th element due to external load, F ext and T ext , contact interactions with neighbouring spheres and other obstacles F cont and T cont , as well as forces and 123 moments resulting from external damping in the system, F damp and T damp , respectively, which can be written as where r c ij is the vector connecting the centre of mass of the i-th element with the contact point with the j-th element.
Similarly as in Inc ICG (2006), the damping terms F damp i and T damp i in Eqs. (15) and (16) in the present work are of non-viscous type and are given by: where a t and a r , are respective damping factors for translational and rotational motion. The damping provides a mechanism to dissipate contact oscillations and represent properly dissipation of the real material. Effect of damping has been investigated by Rojek et al. (2013).

Constitutive Contact Models
The overall behaviour of the system is determined by the contact laws assumed for the particle interaction. Models of contact in the discrete element method can include force and moment interaction between particles. In the present work, however, contact moments are not considered.
Formulation of the constitutive model employs the decomposition of the contact force between two elements F cont (In the next part of this section, indices denoting the elements will be omitted) into normal and tangential components, F n and F t , respectively: where n is the unit vector along the line connecting the centroids of two contacting particles. Modelling of rock or other cohesive materials requires contact models with cohesion allowing tensile interaction force between particles (Potyondy and Cundall 2004;Choi 1992). In the present formulation, rock materials are modelled using the elastic-perfectly brittle model of contact interaction, in which initial bonding between neighbouring particles is assumed. These bonds can be broken under excessive load which allows us to simulate initiation and propagation of material fracture. Contact laws for the normal and tangential direction in the elastic-perfectly brittle model have been investigated by Labra (2012). When two particles are bonded, the contact forces in both normal and tangential directions are calculated from the linear constitutive relationships: where K n is the interface stiffness in the normal direction, k s is the interface stiffness in the tangential direction, u n is the overlap (u n 0) or gap (u n [ 0) at the contact point, u t is the relative displacement at the contact point in tangential direction. Consistently with the sign convention for u n and Eq. (20), the normal force F n is negative in compression and positive in tension. The particle gap/penetration u n is given in terms of the distance between the particle centroids l (also called branch length) and their radii r i and r j u n ¼ l À ðr i þ r j Þ ð 22Þ and the relative tangential displacement u t is updated incrementally where u old t is the vector of the relative tangential displacement from the previous time step rotated to the present contact plane and Du t is the incremental relative tangential displacement with v t being the relative tangential velocity at the contact point. Cohesive bonds are broken instantaneously when the interface strength is exceeded in the normal or tangential direction where R n in the interface strength in the normal direction, and R t the interface strength in the tangential direction. After decohesion, the contact is treated assuming a standard contact model with Coulomb friction. The normal contact force can be compressive only (F n 0), and the tangential contact force is limited by ljF n j kF t k ljF n j ð 27Þ where l is the Coulomb friction coefficient. Two approaches to evaluation of contact model parameters, K t , K n , R t and R n , can be distinguished. In the first approach, the stiffness and strength parameters of the contact model are assumed to depend on the size of contacting particles and are evaluated locally as certain functions of contacting pair radii (Potyondy and Cundall 2004). In the second approach, uniform microscopic properties are assumed in the whole discrete element assembly (Kruyt and Rothenburg 2004). The latter approach is adopted in the present work. Numerical studies performed by Rojek et al. (2012) have shown that uniform global parameters lead to a more brittle behaviour of discrete element models than local size-dependent parameters.

Discrete/Finite Element Method Coupling
It is assumed that the DEM and FEM can be applied in different subdomains of the same body. The DEM and FEM subdomains can overlap each other. The common part of the subdomains is the part where both discretization types are used with gradually varying contribution of each modelling method. This allows us to avoid or minimize unrealistic wave reflections at the interface between the DEM and FEM subdomains. The idea of such coupling follows that used by Xiao and Belytschko (2004) for bridging molecular dynamics with a continuous model.
The FEM subdomain is modelled using the so-called explicit dynamic formulation of the FEM. The explicit FEM is based on the solution of discretized equations of motion written in the current configuration in the following form: where M F is the mass matrix, € u F is the vector of nodal displacements, F ext F and F int F are the vectors of external loads and internal forces, respectively. The form of the FEM equations (28) is similar to the DEM equations (13), which allows us to use the same solution scheme based on the central difference time integration scheme.
The coupling of DEM and FEM subdomains is provided by additional kinematic constraints. The discrete elements in the transitory zone are constrained by the displacement field of overlapping interface finite elements. Making use of the split of the global vector of displacements of discrete elements, u D , into the unconstrained part, u DU , and the constrained one, u DC , u D ¼ fr DU ; r DC g T , additional kinematic relationships can be written jointly in the matrix notation as follows: where N is the matrix containing adequate shape functions. Additional kinematic constraints (29) can be imposed by the Lagrange multiplier or penalty method. The set of equations of motion for the coupled DEM/FEM system with the penalty coupling is as follows where k is the diagonal matrix containing on its diagonal the values of the discrete penalty function, and global matrices M F , M DU , M DC and J D , and global vectors F ext F , F int F , F DU , F DC and T D are obtained by aggregation of adequate elemental matrices and vectors taking into account appropriate contributions from the discrete and finite element parts. Equation (30) is integrated in time using the standard central difference scheme. A detailed description of the coupled DEM/FEM formulation is given by Rojek and Oñate (2007).

Determination of DEM Parameters
The discrete element model can be regarded as a micromechanical material model, with the contact model parameters being micromechanical parameters. Assuming adequate micromechanical parameters, we obtain required macroscopic rock properties. The most important macroscopic rock properties include the Young's modulus E, Poisson's coefficient m, compressive strength r c and tensile strength r t , which will be used for the model calibration in this work. The contact stiffness moduli, K n and K t , and bond strengths, R n and R t , as well as the Coulomb friction coefficient l will be taken as the most significant micromechanical parameters influencing the macroscopic elastic and strength properties.
The parameter identification in the DEM is a typical inverse problem, and it is often solved by a trial-and-error method (Wang and Tonon 2010). Simulations with assumed model parameters are repeated until specific macroscopic properties are reproduced with sufficient accuracy. In the present work, the micromechanical parameters have been determined using the methodology based on the dimensional analysis proposed by Huang (1999) and used later in other works, cf. (Fakhimi and Villegas 2007;Yang et al. 2006;Rojek et al. 2011;Labra 2012). Results of numerical simulations of the standard laboratory tests for rocks, the unconfined compressive strength (UCS) test and the Brazilian tensile strength (BTS), will allow us to establish dimensionless relationships between the contact parameters and the mechanical material properties.
where W represents the dimensionless scale functions, L and A are characteristic lengths and areas of the particle system, and U is a function representing the assembly characterization parameters influence, as the porosity of the particle assembly e, or the average particle radius r. Employing a micromechanical analysis, cf. Rothenburg (2002, 2004), Eqs. (31)-(34), can be redefined as follows (Labra 2012): r c e r 2 R n n c ð1 À eÞ ¼Ŵ c R t R n ; K t K n ; l ð37Þ where the influence of the assembly is considered in the dimensionless numbers by n c as average number of contacts per particle (coordination number) and the solid fraction ð1 À eÞ. The characteristic length e r and area e r 2 are defined based on the size distribution of the particle assembly and the average branch length l over all the bonded contacts of the assembly, as In order to determine the dimensionless scaling functions, a large number of numerical simulations of the UCS and BTS tests have been performed for different specimens. The dimensionless scaling functions obtained from these simulations which will be later used for the estimation of the DEM parameters are presented in Figs. 6 and 7. The influence of the Coulomb friction coefficient in the relationships (37) and (38) has been neglected assuming it affects mainly a post-critical behaviour (Huang and Detournay 2008;Potyondy and Cundall 2004). The relationships plotted in Fig. 7 have been obtained taking l ¼ 0:5.

Numerical Model
Simulation of the LCM unrelieved test has been performed using the model shown in Fig. 8. The disc cutter of diameter 19 00 has been considered. The geometry of the cutting tool has been taken according to the standard constant cross section (CCS) disc profiles used in Herrenknecht AG TBMs. The cutter disc has been discretized with 8880 triangular elements, considering a refinement in the cutter tip in order to reproduce its curvature. It has been treated as a rigid body with all the discretizing nodes slaved to its centre. The disc has a prescribed translational motion and can rotate about its axis of symmetry under the action of the rolling force. A real value of the disc moment of inertia has been used. In this way, a real kinematics and dynamics of the disc cutters during cutting have been reproduced.
The full-scale LCM tests are performed using 1.0 m Â 0.7 m Â 0.7 m block rock samples. In order to obtain results at a reasonable computational cost, a smaller sample with dimensions 0.4 m Â 0.15 m Â 0.4 m has been taken in the numerical model. Numerical tests have shown that it is sufficiently large to avoid boundary effects. The unrelieved rock sample with the adopted boundary conditions is shown in Fig. 8. The displacement of the bottom surface of the sample has been completely restricted, while lateral surfaces have had the out-of-surface displacements restricted.
The DEM discretization has been employed in the zone of interaction with the disc cutter in the subdomain of size 0.15 m Â 0.05 m Â 0.4 m, as shown in Fig. 8. The rest of the sample has been discretized with finite elements taking the 0.01-m-thick overlap zone. 4752 tetrahedral elements have been used in the FE discretization, and the DEM subdomain have been discretized with 35604 spherical particles, with the radius range 1.4-3.9 mm. The characteristic parameters of the particle assembly are shown in Table 1.
Mechanical properties corresponding to granitic gneiss taken from laboratory tests (Labra et al. 2008b) are given in Table 2.
The elastic constants are used in the FEM subdomain. The DEM model parameters, given in Table 3, have been estimated employing the methodology presented in Sect. 4.5, using the properties from Table 2 and considering the particles assembly characterization.
The density of the particles is calculated considering the solid fraction of the assembly, in order to preserve the equivalent mass.
The main cutting process parameters, the velocity and penetration rates, taken from a real TBM drive (Labra et al. 2008b), are given in Table 4. The TBM advance per    cutterhead revolution has been taken as the penetration rate, and the linear velocity of one of the outer disc cutters at considered working conditions has been taken as the cutting velocity for our simulations of the LCM test. The amount of penetration per revolution has been taken as the penetration depth in the LCM model. The penetration depth in all the simulations has been defined as the indentation depth of the cutter tip below the free undamaged surface of the rock sample. It was kept constant under the prescribed cutting tool trajectory. In the LCM test, one of the most important parameters is the spacing between the disc cutters (Rostami and Ozdemir 1993). This parameter cannot be taken into account in the unrelieved cutting. It will be included in simulation of the relieved cutting later on.

Numerical Results
Simulation of the unrelieved cutting is shown in Fig. 9. Material failure induced by the moving disc cutter is illustrated with the damage parameter defined as the ratio of the number of broken bonds to the initial number of cohesive bonds for each particle. The value of the damage parameter ranges from 0 to 1, it is equal to 0 for the undamaged material and 1 for the completely damaged one.
Force-displacement curves and corresponding average forces are plotted in Fig. 10.
The forces have been taken in the range from 0.1 to 0.25 m discarding the results affected by the boundary of the specimen. The average forces obtained in the simulation are given in Table 5 along with available experimental data for comparison.
The experimental normal force has been estimated from in-situ measurements of a real TBM thrust (Labra et al. 2008b). The in situ value of the normal force for one disc has been calculated by dividing the total thrust force by the total number of disc cutters. In principle, the thrust force obtained from a TBM represents the relieved cutting mode Fig. 9 Simulation of the unrelieved LCM test-material damage  Fig. 10 Cutting forces in LCM test in the unrelieved mode because of interaction between cutting grooves in a TBM excavation. However, in this case the cutter spacing was quite large so the interaction between disc cutters was small and the conditions were close to the unrelieved cutting mode.
The average normal force obtained in simulation, 191.08 kN, agrees quite well with the estimated experimental value, 231.8 kN. No comparison with the force predicted by the CSM model has been done because this model considers a spacing between disc cutters, which requires the use of the relieved rock sample.
Using the values of the average normal and rolling forces obtained in simulation, the cutting coefficient CC is evaluated. The value of the cutting coefficient given in Table 5, 0.102, is in the range of the CC values expected for a hard rock (Tarkoy 1983). Obviously, correctness of prediction of the normal forces and cutting coefficients verifies correctness of the predicted values of the rolling forces.
Numerical model allows us to analyse stress distribution in the rock under the forces exerted by the cutting tool. A stress concentration in the rock under the disc can be clearly seen in Fig. 11 showing the distribution of the minimum principal stress (r 3 ) in the longitudinal cross section of the specimen. The stresses in the discrete element model have been evaluated by averaging procedure over representative volume elements taken around each particle (Rojek and Oñate 2007).
The information on stress distribution and evolution during cutting cannot be obtained directly in laboratory tests. Stress analysis is important for investigation of the mechanism of rock cutting and verification of the assumptions made in theoretical models of rock cutting.
The stress concentration under the cutting tool is strictly related to the distribution of the contact pressure. The CSM model developed by Rostami and Ozdemir (1993), Rostami (1997) assumes nearly uniform pressure distribution in the contact zone (see Fig. 5). Stress concentration in a small area close to the attack point observed in Fig. 11 suggests that the contact pressure cannot be uniform in the contact area. This is confirmed in Fig. 12 showing average contact force distribution over the disc cutter. The corresponding pressure distribution is depicted in Fig. 13. The non-uniform distribution obtained in simulations is characterized by a relatively small value (in comparison with a peak value) of the pressure at the attack point, where the disc cutter enters in the contact with the rock and a very small pressure at the rear part of the theoretical contact area. This is contradictory to the assumptions made in the CSM model by Rostami and Ozdemir (1993), Rostami (1997), but it agrees quite well with the results of experimental investigations of the pressure distribution over the cutter discs published by Rostami (2013).
The force distribution plotted in Fig. 13 allows us to estimate the angle b determined by the direction of the resultant interaction force. The b angle and cutting coefficient (CC) obtained in the simulation are compared with the respective values predicted by the CSM model in Table 6. The theoretical values of the angle b according to the CSM model have been calculated as proposed by Gertsch et al. (2007), Rostami (2008) It can be noted that the values of the angle b and cutting coefficient (CC) obtained in the simulation are very close to those estimated by the CSM model.

Sensitivity Analysis of Unrelieved Cutting
Influence of various parameters on the performance of the rock cutting process has been analysed. Sensitivity of the results to changes of process and design parameters is very important for an optimization of the process. Two different categories of the parameters involved in the rock cutting process can be distinguished. The first category is related to the mechanical properties of the rock specimen. Most of the studies reveal that the most relevant mechanical properties of the rock are the compressive and tensile strengths. The elastic constants of the rock are not considered important for the estimation of the cutting forces and are not analysed in this section. The second category involves the geometric settings of the cutting process. Here, we can include the profile of the disc cutter, the spacing between discs and penetration. Sensitivity analysis of unrelieved cutting has been performed for the uniaxial compressive strength, cutting velocity and penetration depth. The effect of disc spacing will be analysed later for relieved cutting.

Effect of the Uniaxial Compressive Strength
Compressive strength of the rock material is one of the most significant mechanical parameters for the performance of the cutting process. Figure 14 shows the relationship between the average normal and rolling forces and the uniaxial compressive strength obtained in numerical simulations.
The trend line shows a linear increase in the forces with the growth of the compressive strength. In the simulation, values of the UCS are estimated using the dimensionless scale functions (Eqs. 37, 38), assuming a constant ratio between the compressive and tensile strength n ¼ r c =r t .
The linear force-strength relationship is also predicted by the CSM model. Introducing the strength ratio n ¼ r c =r t into Eq. (3) we obtain Keeping the ratio n and other parameters in Eq. (15) constant, we have the linear dependence of the pressure and resulting forces on the compressive strength. Figure 15 shows the effect of the compressive strength on the value of the cutting coefficient. In contradiction to the CSM and other theoretical models in the literature, which relate the cutting coefficient with geometric parameters, such as the penetration and disc cutter radius, Fig. 15 shows a linear relationship between the compressive strength and the cutting coefficient.

Effect of the Penetration Depth
The penetration depth, or penetration rate, is one of the most important geometric parameters, together with the disc spacing, that affect the performance of the TBM, as it is directly related with the advance rate of the TBM  cutterhead. All the theoretical models employ this parameter in evaluation of cutting forces as well as in estimation of the volume of rock material excavated and amount of energy required for the excavation process.
The results plotted in Fig. 16 show the influence of the penetration depth on the normal and rolling forces for the disc cutter of diameter 19 00 in both materials.
When the penetration depth increases, the normal and rolling forces also increase. The influence of the penetration depth can be clearly seen in the value of the rolling force, that is directly related to the energy required in the cutting process. Figure 17 shows the comparison of the cutting coefficient obtained in the simulations and estimated by the CSM model as functions of the penetration depth. Estimation of the cutting coefficient in the CSM model is based on purely geometric criteria (see Eq. (5)), while the numerical results indicate a certain influence of material properties.
6 Simulation of Relieved Rock Cutting

Numerical Model
All the results presented until now have been obtained for unrelieved cutting conditions. Rock cutting with disc cutters mounted in a TBM cutterhead is performed mainly in relieved cutting conditions influenced by the previous passes of disc cutters. Interaction of adjacent cutting paths affects directly the chip formation (Fig. 18) and ensures an efficient cutting performance; therefore an optimum spacing between disc cutters is one of the most important parameters in the design of TBMs for hard rock conditions. The optimum spacing between adjacent cutting trajectories for a given penetration of cutters can be determined by minimization of the energy required for the chip formation.
The geometry used to analyse relieved cutting is shown in Fig. 19. The size of the specimen is similar to that of the unrelieved specimen shown in Fig. 8. A larger DEM subdomain has been used in order to reproduce the spacing up to 0.13 m.
The parameter l in Fig. 19 changes depending on the spacing between the disc cutter and the previous pass in order to maintain a sufficient distance between the disc cutter and the DEM/FEM coupling interface. Two different specimens, for penetration of 4 and 8 mm, have been generated. The DEM subdomains have been discretized with 51668 and 50714 particles, respectively. The main parameters characterizing the particles assemblies are summarized in Table 7. The DEM parameters for the new rock specimens are summarized in Table 8.   Results of the numerical simulations have been compared with the predictions of the CSM model. The normal forces obtained in the simulation are slightly higher than those evaluated using the CSM model for both cases of penetration. So are the numerical rolling forces higher than the theoretical predictions although in this case the difference is smaller.

Numerical Results
Numerical and theoretical values of the cutting coefficient are presented in Fig. 21 as functions of disc spacing for penetration depth of 4 and 8 mm. Numerical results confirm the cutting coefficient is independent of the spacing as it is postulated by the CSM model. The cutting coefficients estimated with the CSM model are higher than those computed in the simulations. A higher difference between the numerical and theoretical values can be seen for the lower penetration. Similarly, the difference between the numerical and theoretical values of the b angle is higher for the lower penetration. The theoretical and numerical values of the cutting coefficient and angle b spacing of 80 mm are given in Table 9.
The present model and the CSM models are different in nature, i.e. one is micromechanical and the other empirical, so the differences in their predictions have different sources. Some differences between the CSM model and the simulations can be explained by simplified assumption for the pressure distribution made in the CSM model. Similarly as in case of the unrelieved rock cutting, the distribution of pressure over the disc cutter obtained in numerical simulation differs from the uniform pressure distribution assumed in the CSM model. Non-uniform distributions of the interaction forces over the disc cutters for both penetration cases are presented in Fig. 22, for the spacing of 80 mm.
Having determined the rolling force F r for the given penetration depth p and disc cutters spacing s, the specific energy SE can be calculated from Eq. (10). The specific energy as a function of the spacing/penetration ratio is plotted in Fig. 23.
The specific energy computed in simulation of relieved cutting is compared in this figure with the specific energy estimated by the CSM model for different values of penetration and spacing. A good correlation between the CSM model and the simulation results has been obtained.
The specific energy estimated in simulations agrees quite well with field data reported in the literature. Bilgin et al. (2014) report the field specific energy for the Kadikoy-Kartal Metro where the mean compressive strength of excavated rocks was 50 MPa. The field specific energy varied in the range 7-18 kWh/m 3 depending on the penetration rate. These values are slightly below the values given in Fig. 23, but this is understandable since the compressive strength in our case was higher (147.3 MPa). The field specific energy reported also by Bilgin et al. (2014) for the Beykoz tunnel excavated in rocks with higher compressive strength (100 MPa) varies from 5 to 27 kWh/m 3 . The field data show that the numerical simulations give reasonable predictions of the specific energy. This is an important result of our simulations since confirms a good performance of the model and its potential utility in real applications.

Conclusions
Numerical tests have demonstrated a good performance of the coupled discrete/finite element model of rock cutting. A good agreement of numerical results with experimental measurements and theoretical predictions has been found.  Main parameters characterizing rock cutting with a TBM disc cutters such as cutting and rolling forces, cutting coefficient and specific energy have been estimated correctly in numerical simulations. The numerical model is capable to represent properly complexity of rock cutting with TBM disc cutters. The numerical simulation can provide valuable information about the cutting phenomenon such as stress distribution in the rock and contact pressure distribution at the tool-rock interaction area. Numerical results have confirmed a non-uniform contact pressure distribution revealed in experimental investigations and shown that a uniform pressure distribution in theoretical models is a simplified assumption. The numerical model has provided a new insight into the cutting process enabling us to understand better rock cutting mechanism.
Sensitivity analysis of rock cutting performed for different parameters including disc geometry, cutting velocity, disc penetration and spacing has shown that the presented numerical model is a suitable tool for the design and optimization of rock cutting process.