3D Continuum Modelling of PDC Cutting of Rock with a Simple Contact-Erosion Scheme

Featured Application: The numerical method presented in this paper is designed to estimate the cutting force and the removed material volume during rotary rock drilling by a single PDC cutter. Abstract: This paper presents a relatively simple numerical approach to predict the cutting force during PDC (polycrystalline diamond contact) cutting of rock. The rock failure model is based on a damage-viscoplasticity model, with the Drucker–Prager yield surface and the modiﬁed Rankine surface as the tensile cut-off. The damage part of the model has separate scalar damage variables for tension and compression. The PDC cutter is idealized to a rigid surface and its interaction with the rock is modelled by contact mechanics, while solving the global equations of motion explicitly in time. A damage-based erosion criterion is applied, to remove the contact nodes surrounded by heavily damaged elements. The eroded elements are left in the mesh as ghost elements that do not contribute to the load transfer but preserve the mass conservation. Numerical simulations on granite, demonstrate that the method reliably predicts the cutting force of a single PDC cutter at different cutting depths and rake angles.


Introduction
Drilling a well into a deep reservoir in a hard formation is technologically challenging and incurs major costs in hydrocarbon operations and harvesting geothermal energy. The major problem therein is the tool wear, which results in long idle times, due to the bit replacement. For this reason, the drill bit (see Figure 1) cutters are made of extremely hard materials. Presently, the leading drilling technology in hard formations (e.g., granite) uses the so-called polycrystalline diamond compact (PDC) cutters, made of synthetic polycrystalline diamonds bonded to a tungsten carbide blade. The drilling action by cutting is achieved by applying a certain weight on the bit (WOB) and then rotating the bit with a torque, T. Design of such drill bits for the extreme conditions, involving high pressure, faced in deep well drilling, requires good understanding of the tool-rock interaction dynamics.
A numerical model, implemented in a computer code capable of predicting the cutting force, the initiation of fracture, and the progression of fragmentation leading to rock material removal, is of utmost importance in the research and development of PDC cutters and rotary bits. There are basically two numerical approaches in the literature related to rock cutting, the discrete element method (DEM), or particle methods more generally, and the finite element method (FEM). The DEM was applied, e.g., in [16][17][18]21,22,25], while ref. [27] gives a review on numerical methods for chip formation in metal cutting. Being based on discontinuum assumption, this method naturally allows fragmentation and crushing to take place in the simulations of rock cutting, as the rock medium to be modelled is discretized into deformable or rigid particles with interaction laws between them. The major drawback is the computational labor, especially in 3D models, required to track the configurations of the particles (neighbors of a particle are not fixed) and possible contacts between them. The FEM, in contrast, is a continuum method, where the discretization results in elements with fixed neighbors. Therefore, modelling large deformations and fragmentation requires special techniques, such as remeshing, element erosion, and contact modelling. However, the FEM is a mature numerical method, that has been successfully used to simulate rock cutting problems, as exemplified in [19,20,23,24,26]. Therefore, as the computational effort needed is less intensive, and the calibration of the model parameters is less involved, the FEM is chosen for the present study.
The above mentioned numerical studies use commercial multipurpose FEA software packages (Abaqus, Ansys, etc.) for rock cutting simulation. The present study aims to develop a 3D numerical platform, tailored to simulate rock cutting with PDC cutters. At this preliminary stage, only the single cutter case is considered. The target is to reliably predict the cutting force and the volume of removed rock in a reasonable CPU time, say over a night or two, for a cutter with a variable size, rake angle, and cutting depth. The method is implemented in the Matlab (alternatively GNU Octave) language, with Fortran based mex-files for the most computationally laborious parts. The main novelty of the approach is the unique combination of the Lee-Fenves type of stiffness recovery scheme for damage, the forward increment Lagrange multiplier method to solve the contact forces in explicit time stepping, and a contact node erosion algorithm. A judicious coupling of these components with a calibrated rock material model allows us to meet the set target, i.e., to reliably predict the cutting force in rock cutting with PDC cutters. A numerical model, implemented in a computer code capable of predicting the cutting force, the initiation of fracture, and the progression of fragmentation leading to rock material removal, is of utmost importance in the research and development of PDC cutters and rotary bits. There are basically two numerical approaches in the literature related to rock cutting, the discrete element method (DEM), or particle methods more generally, and the finite element method (FEM).
The DEM was applied, e.g., in [16][17][18]21,22,25], while ref. [27] gives a review on numerical methods for chip formation in metal cutting. Being based on discontinuum assumption, this method naturally allows fragmentation and crushing to take place in the simulations of rock cutting, as the rock medium to be modelled is discretized into deformable or rigid particles with interaction laws between them. The major drawback is the computational labor, especially in 3D models, required to track the configurations of the particles (neighbors of a particle are not fixed) and possible contacts between them. The FEM, in contrast, is a continuum method, where the discretization results in elements with fixed neighbors. Therefore, modelling large deformations and fragmentation requires special techniques, such as remeshing, element erosion, and contact modelling. However, the FEM is a mature numerical method, that has been successfully used to simulate rock cutting problems, as exemplified in [19,20,23,24,26]. Therefore, as the computational effort needed is less intensive, and the calibration of the model parameters is less involved, the FEM is chosen for the present study.
The above mentioned numerical studies use commercial multipurpose FEA software packages (Abaqus, Ansys, etc.) for rock cutting simulation. The present study aims to develop a 3D numerical platform, tailored to simulate rock cutting with PDC cutters. At this preliminary stage, only the single cutter case is considered. The target is to reliably predict the cutting force and the volume of removed rock in a reasonable CPU time, say over a night or two, for a cutter with a variable size, rake angle, and cutting depth. The method is implemented in the Matlab (alternatively GNU Octave) language, with Fortran based mex-files for the most computationally laborious parts. The main novelty of the approach is the unique combination of the Lee-Fenves type of stiffness recovery scheme for damage, the forward increment Lagrange multiplier method to solve the contact forces in explicit time stepping, and a contact node erosion algorithm. A judicious coupling of these components with a calibrated rock material model allows us to meet the set target, i.e., to reliably predict the cutting force in rock cutting with PDC cutters.

Numerical Methods
The numerical method for simulating rock cutting includes a description of fracturing rock based on a damage-viscoplasticity model, and an explicit time stepping scheme to solve the global equations of motion with frictional contact boundary conditions. A simple erosion scheme is also presented, to deal with the heavily damaged elements.

Rock Failure Model
The rock is described as an isotropic homogeneous material up to the elastic limit, after which the material softens, and its strength and stiffness degrade. Temperature effects are, as usual in numerical studies on PDC cutting, neglected. However, it should be noted that considerable heat generation may occur due to friction, especially in rock cutting with conical picks [8]. In cutting, on the one hand, the rock material undergoes large deformations. On the other hand, rock, as a brittle material, cannot withstand large deformations, at least in normal conditions. Therefore, it seems to be enough if the constitutive model can describe finite rigid body rotations. For this reason, the St. Venant-Kirchhoff material model is applied to write the elasticity part of the rock model as: where σ is the stress tensor, C e is the elasticity tensor, ε GL is the Green-Lagrange strain tensor, ε vp is the viscoplastic strain tensor, and C is the right Cauchy-Green deformation tensor, defined by the deformation gradient F, and 1 is the second order identity tensor. The colon (:) in (1) denotes the double contraction operator for tensors, i.e., A : A = A ij A ij where repeated indices imply summation. The stress in (1) should be the 2nd Piola-Kirchhoff stress, S, but under the present small deformation framework, all stress measures reduce to the same numerical values. The softening process, in both compression and tension, is governed by a phenomenological damage-viscoplasticity model. For present purposes, the viscoplasticity part is formulated with the linear Drucker-Prager surface, cut by the rounded Rankine criterion, to match the tensile strength. Mathematically, these surfaces read: where the symbol meanings are as follows: I 1 and J 2 are the first and the second invariants of the stress tensor σ, respectively; c 0 and σ t0 are the intact cohesion and tensile strength; σ i is the ith principal stress, from which only the positive part is considered through McAuley brackets in the MR criterion (5); . κ DP , .
κ MR are the rates of the internal variables in compression and tension; s DP , s MR are the constant viscosity moduli; α DP = 2 sin φ/(3 − sin φ), k DP = 6 cos φ/(3 − sin φ) are DP parameters defined with the internal friction angle φ, so as to match the uniaxial compressive strength. Equation (4) is the viscoplastic potential, where β DP is a parameter, controlling the dilation, similar to α DP but defined with the dilation angle ψ. Thermodynamic consistency is guaranteed for ψ ≤ φ [28], while associated viscoplasticity is recovered when ψ = φ. Equation (6) accounts for the strain rate dependency of the rock.
The damage part of the model has separate scalar damage variables for tensile and compressive stress regimes. The damage variables and evolution equations are:  (9) . .
where the symbols are: ω t and ω c are the tensile and compressive damage variables, respectively; A t and A c control the final value of the damage variables; β t and β c control the initial slope and the amount of damage dissipation, and are defined by the fracture energies G Ic and G IIc ; h e is a characteristic length of a finite element; ε vp eqvt is the equivalent viscoplastic strain in tension defined, in the rate form, as the trace of the viscoplastic strain rate tensor,  (8) is Koiter's rule for bi-surface plasticity, while Equation (12) defines the consistency conditions. The present formulation is thus the viscoplastic consistency approach, used by Wang et al. [29].
The final model component is the nominal-effective stress relation, which specifies how the damage variables operate on the stress. The Lee and Fenves [30] relation is chosen for this purpose: where σ is the nominal stress tensor in Equation (1); s t and s c are stiffness recovery functions depending on the principal nominal stresses, σ i ; and parameters w t and w c control the degree of recovery. The damage and viscoplasticity parts of the model are combined in the effective stress space [31]. With this formulation, the viscoplasticity and damage computations can be performed independently, so that the stress return mapping is first performed in the effective principal stress space by standard methods [29]. Then, the damage variables are updated and, finally, the nominal stress is calculated by Equation (13).

Simulation Model for Single PDC Cutter
The principle of the rock cutting model is shown in Figure 2. As the aim is to predict the cutting force, the PDC cutter is idealized to a rigid body, so that its behavior can be described by a single node while the cutter geometry is virtual. The node representing the cutter is given a constant velocity boundary condition, v0 in Figure 2. Moreover, a viscous damper (dashpot) is also attached to this node, to absorb the possible wave effects due to the cutting velocity. The damper properties are defined by the computational area, Ab, of the cutter, the density ρb, and the bar velocity, / , with Eb being the Young's modulus of the bit. The forces Fbit and Fc are the force resisting the bit movement and its reaction, the cutting force. In addition, β is the rake angle and d is the cutting depth. Finally, ph and pv are the horizontal (confining) and vertical (overburden, e.g., from the drilling mud) pressures prescribed for simulating deep with E b being the Young's modulus of the bit. The forces F bit and F c are the force resisting the bit movement and its reaction, the cutting force. In addition, β is the rake angle and d is the cutting depth. Finally, p h and p v are the horizontal (confining) and vertical (overburden, e.g., from the drilling mud) pressures prescribed for simulating deep drilling.
The equation of motion for the cutter node is written as where m b is a computational mass for the cutter, .. u b,z and . u b,z are the acceleration and velocity of the cutter in the cutting (z) direction, and F µ is the friction force. This equation is added to the equation of motion of the rock model and solved simultaneously with it, as explained next.

The Golobal Equations of Motion for Cutter-Rock System with Contact Constraints
As implicit time integration is not a feasible option in modelling cutting with severely distorted elements and erosion schemes (element deletion), explicit time stepping based on the modified Euler method [32] is chosen. Coupled with the forward increment Lagrange multiplier method [33], the system of equations governing the present problem is written as [34]: M ..
where the symbols are as follows: A is the standard finite element assembly operator; M is the mass matrix to be diagonalized by the row sum technique; C is the material damping matrix; f int t is the internal nodal force vector, where B T e is the transpose of the gradient of the interpolation function matrix N e ; f ext t is the external nodal force vector including the nodal forces from the confining and overburden pressures; f dp t is the damping force due to the viscous damper in Equation (16); G is the normal contact constraint matrix, having the normal vectors n bit = [0, − sin β, cos β ] T , n rock = −n bit , of the cutter and the rock contact nodes, as its components; b is the initial distance between all the potential contact nodes in the rock and the cutter node; λ t is the Lagrange multiplier vector representing the normal contact force; f µ,t is the friction force vector, calculated as a sum for N conts active contact nodes; µ is the friction coefficient; F e N is the normal contact force for node e (component nd of λ t ); g tan,nd is the relative tangential motion at contact node nd, depending on the unit tangent vectors t i and the relative motion vector u rel .
It should be noted that only the slip condition is accounted for in the present application of rock cutting, with a constant velocity BC. Moreover, as can be noticed in Equation (18), the impenetrability constraint is written at t + ∆t, while Equation (17) is written at t. This is the essential trick of the forward increment Lagrange multiplier method, which results in the following predictor-corrector scheme to solve the previous equations [34]: ..
. u t+∆t = The solution with this equation system is as follows: First, the acceleration and displacement are predicted with (25) and (26) as if there were no contacts. Then, if contacts are detected (interpenetration occurs so that the R.H.S. of (27) has positive entries), the normal contact force, i.e., the Lagrange multiplier vector, is solved with (27). Finally, the acceleration is corrected with (22), while taking into account the friction contact force due to sliding by (21). A new velocity and displacement can then be predicted with (23) and (24).

Erosion Scheme
Element deletion or erosion schemes are mainly used to get rid of highly distorted elements in explicit analyses, to reach numerical convergence and subsequently proceed with the simulation. In addition, deletion of elements is also used to visualize fracture propagation. Highly distorted or damaged elements, including the nodes becoming orphans, can either be removed or left in the mesh as ghost elements that cannot bear any loadings. The former approach is more involved, as it requires restructuring and an assembly of the related data structures and computational entities, such as mass and stiffness matrices. In the present case, the latter option is chosen, and the elements marked to be deleted behave like a vacuum, i.e., the R.H.S. of the equation of motion is zero. However, the mass matrix is kept untouched, to secure continuation of the analysis (zeros on the diagonal would result in a halt).
An erosion criterion indicating the elements to be marked as deleted is required. Many options exist in the literature, such as strain, stress, and damage-based criteria [35]. A damage-based criterion is selected for the present purposes in the following form: where ω tresh is a threshold value above which an element is regarded as destroyed. This criterion is a conservative one, requiring that both damage variables reach the threshold.
In the present case of leaving the deleted elements in the mesh, the erosion scheme thus involves removing the contact nodes for which all neighboring elements are marked as deleted. Such a node, being left as an orphan in the space, cannot bear any loading and should thus be removed from the list of contact nodes. This can be achieved with the pseudocode in Algorithm 1. In the algorithm, cn n is a contact node n and N conts is the number of contact nodes. Moreover, a list of all eroded elements, in addition to the one in the algorithm, is generated, so that there is no need to calculate strains and stresses for an element if it is found in the list. In addition to updates in the algorithm, G and b are also updated based on the value of the gap function, i.e., only the rows for which the gap function Gũ t+∆t − b is positive are used in Equations (22) and (27) at each time step. This fulfils the unilateral condition that contacts can carry only compressive loading.

Numerical Simulations
In this section, numerical simulations demonstrating the method are carried out. The validation cases are the single PDC cutter experiments on granite from [4]. First, the material properties and model parameters, along with the finite element mesh, are specified. Then, the validation simulations are carried out. Finally, the effect of confining pressures is tested. As the variation of cutting speed has a negligible effect on the cutting force at low cutting speeds [11], it is not tested here.

Material Properties and Model Parameters
The material properties and model parameters are given in Table 1. The cohesion and friction angle in Table 1 give a uniaxial compressive strength of 165 MPa. Moreover, the viscosity moduli given are chosen, on the one hand, to be small enough so as to not cause any significant strain rate effects, and, on the other hand, to be large enough that they have a stabilizing effect on the stress return mapping at the material point level. Furthermore, the parameters w t and w c , controlling the degree of stiffness recovery upon stress reversals, are chosen so as to match the cutting experiments. Finally, the cutting speed is a compromise between computational labour of the explicit time integration and the requirement of keeping it as close as possible to what was used in the experiments, 5 mm/s. The stable time step for the present simulations is 3.15 × 10 −8 s, so that simulating a 10 mm stroke of cutting at 5 mm/s would require more than 600 million time steps. As performing the analyses of a single time step takes~0.4 s of CPU time, using such a cutting speed is not feasible. This is a common problem in modelling cutting problems with explicit codes, as exemplified by Jaime et al. [20], who used a cutting speed of 4 m/s in their finite element simulations.
The finite element mesh, made of 523,928 linear tetrahedrons, is shown in Figure 3. The mesh has a 10 mm long initial groove, roughly representing a cut groove. The purpose is to avoid excessive edge effects and, thus, to reach stabilized conditions earlier in the cutting process.

Validation Simulations
The validation simulations concern comparing the cutting force predicted with the present method, against the experiments from [4], at cutting depths of 0.3, 0.9, and 1.5 mm. Moreover, the rake angle is 20 • , while the radius of the cutter is 9.5 mm. Figure 4 shows the results after a 12 mm stroke.
In each case, there are heavily damaged, both in tensile and compressive (shear) mode, regions, representing the grooves of removed rock. Most of the elements in these regions are eroded (ghost) elements, having interaction neither with the cutter nor with the intact rock. The related cutting force graphs exhibit huge initial peaks, reaching 3.6 kN even at the smallest cutting depth, due to the smooth initial contact surface causing massive fracture events at its edge. The finite element mesh, made of 523,928 linear tetrahedrons, is shown in Figure 3. The mesh has a 10 mm long initial groove, roughly representing a cut groove. The purpose is to avoid excessive edge effects and, thus, to reach stabilized conditions earlier in the cutting process.

Validation Simulations
The validation simulations concern comparing the cutting force predicted with the present method, against the experiments from [4], at cutting depths of 0.3, 0.9, and 1.5 mm. Moreover, the rake angle is 20°, while the radius of the cutter is 9.5 mm. Figure 4 shows the results after a 12 mm stroke.

Validation Simulations
The validation simulations concern comparing the cutting force predicted with the present method, against the experiments from [4], at cutting depths of 0.3, 0.9, and 1.5 mm. Moreover, the rake angle is 20°, while the radius of the cutter is 9.5 mm. Figure 4 shows the results after a 12 mm stroke.  After the initial stage, oscillations in the cutting force are much milder, representing more local fracture events. The actual removed or ghost elements are plotted in Figure 4. The volume of removed rock material, represented as the volume of the elements where both damage variables exceed the threshold 0.99, increases linearly as a function of the cutting depth, as can be easily calculated from the values given in Figure 5. Unfortunately, the authors in [4] do not provide data to enable comparison to their experiments. However, the circular cross-sectional shape of the removed element volumes in Figure 5, and the bottom shape of the cut grooves in the rock, agree with the experimental grooves cut into granite with traditional circular PDC cutters in [11].
The volume of removed rock material, represented as the volume of the elements where both damage variables exceed the threshold 0.99, increases linearly as a function of the cutting depth, as can be easily calculated from the values given in Figure 5. Unfortunately, the authors in [4] do not provide data to enable comparison to their experiments. However, the circular cross-sectional shape of the removed element volumes in Figure 5, and the bottom shape of the cut grooves in the rock, agree with the experimental grooves cut into granite with traditional circular PDC cutters in [11].  Figure 6 shows the predicted cutting force vs. cutting length response, compared to the experimental one.   Figure 6 shows the predicted cutting force vs. cutting length response, compared to the experimental one.
The volume of removed rock material, represented as the volume of the elements where both damage variables exceed the threshold 0.99, increases linearly as a function of the cutting depth, as can be easily calculated from the values given in Figure 5. Unfortunately, the authors in [4] do not provide data to enable comparison to their experiments. However, the circular cross-sectional shape of the removed element volumes in Figure 5, and the bottom shape of the cut grooves in the rock, agree with the experimental grooves cut into granite with traditional circular PDC cutters in [11].  Figure 6 shows the predicted cutting force vs. cutting length response, compared to the experimental one.  The predicted and the experimental cutting force vs. length curves, have similar qualitative features, but the oscillations are somewhat larger in the predicted curve. Moreover, the average of the cutting force is 565 N for the predicted and 761 N for the experimental.

Effect of the Stiffness Recovery Parameter w c
The effect of the w c parameter, which controls the stiffness recovery effect on the tensile damage variable (Equation (13)), was tested. The model was calibrated by adjusting this parameter so that the experimental cutting force at the cutting depth of 0.9 mm was matched (the predictions are compared to experiments in Section 4).
The simulation results setting w c = 1 (full recovery) are shown in Figure 7. With full recovery, the tensile damage variable has no effect in compression, and hence the cutting force is substantially larger than in the case of w c = 0.15. However, the volume of removed (ghost) elements (Figure 7b) is only 2.4% larger with full recovery. Therefore, the parameter w c seems to have a minor effect on the failure mode. Notwithstanding, the number of elements with low values of tensile damage is much larger here than it is in Figure 4b. was matched (the predictions are compared to experiments in Section 4).
The simulation results setting 1 (full recovery) are shown in Figure 7. With full recovery, the tensile damage variable has no effect in compression, and hence the cutting force is substantially larger than in the case of 0.15. However, the volume of removed (ghost) elements (Figure 7b) is only 2.4% larger with full recovery. Therefore, the parameter seems to have a minor effect on the failure mode. Notwithstanding, the number of elements with low values of tensile damage is much larger here than it is in Figure 4b.

Effect of Rake Angle
The effect of rake angle, β, (see Figure 2) was tested next. Other parameters were kept unaltered. Figure 8 shows the results in the case where d = 0.9 and β = 10°.

Effect of Rake Angle
The effect of rake angle, β, (see Figure 2) was tested next. Other parameters were kept unaltered. Figure 8 shows the results in the case where d = 0.9 and β = 10 • .  The damage patterns in Figure 8 do not exhibit major differences from the case with β = 20°, and the volume of removed elements differs by only 2.3%. However, the cutting forces in Figure 8c show some deviations.

Effect of Confinement
The final simulations concern deep drilling conditions, where high pressures, due to the overburden rock and/or the drilling mud pillar, as well as horizontal tectonic stresses, prevail. Two cases are tested. In the first (Conf1), the overburden pressure (pv in Figure 2) is 50 MPa and the tectonic, horizontal pressure is 25 MPa; while in the second case (Conf2) pv = 20 MPa and ph = 50 MPa. Referring to Figure 3, the overburden pressure is applied to the y-direction and the horizontal pressure in the x-and z-directions.
According to the simulations results, in Figure 9, the removed volumes are 106.9 mm 2 and 94.7 mm 2 in the confined cases Conf1 and Conf2, respectively, while 91.9 mm 3 of rock removal is predicted in the unconfined case. The similarity of the cut grooves and the volume in the unconfined and these confined cases, means that no significant cracking events outside the cutter face area occurred during the simulations (confinement usually suppresses spall cracks). Moreover, in Conf2, the cutting force peaks are generally higher than in Conf1. Therefore, it seems that Conf2 is more detrimental. In confined conditions, a much higher cutting force is required to remove the same volume of material. The damage patterns in Figure 8 do not exhibit major differences from the case with β = 20 • , and the volume of removed elements differs by only 2.3%. However, the cutting forces in Figure 8c show some deviations.

Effect of Confinement
The final simulations concern deep drilling conditions, where high pressures, due to the overburden rock and/or the drilling mud pillar, as well as horizontal tectonic stresses, prevail. Two cases are tested. In the first (Conf1), the overburden pressure (p v in Figure 2) is 50 MPa and the tectonic, horizontal pressure is 25 MPa; while in the second case (Conf2) p v = 20 MPa and p h = 50 MPa. Referring to Figure 3, the overburden pressure is applied to the y-direction and the horizontal pressure in the xand z-directions.
According to the simulations results, in Figure 9, the removed volumes are 106.9 mm 2 and 94.7 mm 2 in the confined cases Conf1 and Conf2, respectively, while 91.9 mm 3 of rock removal is predicted in the unconfined case. The similarity of the cut grooves and the volume in the unconfined and these confined cases, means that no significant cracking events outside the cutter face area occurred during the simulations (confinement usually suppresses spall cracks). Moreover, in Conf2, the cutting force peaks are generally higher than in Conf1. Therefore, it seems that Conf2 is more detrimental. In confined conditions, a much higher cutting force is required to remove the same volume of material.

Effect of Confinement
The final simulations concern deep drilling conditions, where high pressures, due to the overburden rock and/or the drilling mud pillar, as well as horizontal tectonic stresses, prevail. Two cases are tested. In the first (Conf1), the overburden pressure (pv in Figure 2) is 50 MPa and the tectonic, horizontal pressure is 25 MPa; while in the second case (Conf2) pv = 20 MPa and ph = 50 MPa. Referring to Figure 3, the overburden pressure is applied to the y-direction and the horizontal pressure in the x-and z-directions.
According to the simulations results, in Figure 9, the removed volumes are 106.9 mm 2 and 94.7 mm 2 in the confined cases Conf1 and Conf2, respectively, while 91.9 mm 3 of rock removal is predicted in the unconfined case. The similarity of the cut grooves and the volume in the unconfined and these confined cases, means that no significant cracking events outside the cutter face area occurred during the simulations (confinement usually suppresses spall cracks). Moreover, in Conf2, the cutting force peaks are generally higher than in Conf1. Therefore, it seems that Conf2 is more detrimental. In confined conditions, a much higher cutting force is required to remove the same volume of material.

Discussion
The results are discussed and compared to the experiments. First, the results are quantified in Table 2, with some calculated drilling technical indexes. More specifically, the work done by the cutter, w bit , is the area under the cutting force vs. cutting length curve. In addition, the MSE is the mechanical specific energy, i.e., E MS = w bit /V crater [11]. The cutting force, F c , is defined, not as the average of the force vs. length curve, but as the average of the peaks (as in [4]), marked by red dots in Figure 4a. It should be noted that the peaks included to calculate the cutting forces start at the cutting length of 2 mm, while for E se , the whole response is used, i.e., starting from zero cutting length. The model predicts the cutting force at higher cutting depths with an acceptable accuracy. The best prediction is naturally obtained at the cutting depth of 0.9 mm, because the model was calibrated to that case (by adjusting w c ). At the smallest cutting depth, the error, 48%, is unacceptably large, but at 1.5 mm, the error is acceptable at both tested rake angles. No clear data in [4] was provided on cutting forces under confined conditions. Anyway, confinement results in a substantial increase,~50% here, in the cutting force. Moreover, according to the present predictions, the case with higher horizontal than vertical pressure is more detrimental, i.e., requires a higher cutting force to remove the same volume of rock.
As for the other rock cutting parameters, the work done on the rock by the cutter naturally increases as the cutting depth increases (higher force moves the same cutting length). Interestingly, the MSE does not deviate much at the three tested cutting depths. The reason for this is that, in the present simulations, the volume of removed rock increases in proportion to the work done into the rock, so that the ratio w bit /V crater is almost the same in each case. Confinement increases the MSE, which means that drilling under confinement is inefficient.
The weaknesses of the present approach include the inability to properly describe cracking, and, especially, fragmentation. The standard continuum model, based on damage and viscoplasticity theories, describes a crack in the smeared (integral) sense, smoothing the crack tip effects. Moreover, the PDC cutter itself is idealized to a single node with rigid virtual geometry. Therefore, the only information the present approach yields for the cutter is the cutting force. However, cutter geometry, cutting speed, rake angle, and cutting depth can be modelled. In particular, the present code can easily be extended to simulate hybrid drilling methods, such as PDC cutting combined with a high-pressure CO 2 jet [21].

Conclusions
A numerical method for predicting the cutting force and the volume of removed material, in rock cutting with a single PDC cutter, was developed and validated in this paper. There follow some concluding remarks:

•
The numerical method includes a damage-viscoplasticity model to describe rock fracture in PDC cutting, and an explicit time stepping approach to solve the global equations of motion with frictional contact boundary conditions. The code is implemented in the Matlab (alternatively Octave) language with a mex-file (compiled from Fortran native code) to speed up the loop over element level computations. As such, the code can perform a 12 mm stroke of PDC cutting with a single cutter (9.5 mm diameter), on a mesh with~500,000 linear finite elements, in 2 days on a typical powerful laptop (time steps being~3 × 10 −8 s). If parallelized, a longer stroke or multiple cutters within the same CPU time are feasible. • An erosion algorithm based on a conservative damage criterion was presented. Accordingly, an element is considered eroded if both the tensile and compressive damage variables exceed a threshold value. However, the elements are not removed from the mesh but retained as ghost elements, with no contribution other than masses in the lumped mass matrix, to preserve the balance of mass. Moreover, a node in the potential or actual contact nodes is removed from the list of contact nodes if all elements connected to it are deleted (i.e., are ghost elements).

•
The method predicts the cutting force in PDC cutting of granite for a single cutter with a good accuracy at variable cutting depths and rake angles. The cut grooves have realistic features. As such, the method serves as a possible platform to be extended to simulate hybrid PDC cutting, where the mechanical cutting action is facilitated by some nonmechanical agent, e.g., thermal jets or microwaves. This is a topic for future developments of the code.