## 2 Numerical Methods
### 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]:  
\[\mathbf{M}\ddot{\mathbf{u}}_{t}+\mathbf{C}\dot{\mathbf{u}}_{t}+\dot{\mathbf{f }}_{t}^{\mathrm{int}}=\mathbf{f}_{t}^{\mathrm{ext}}-\mathbf{f}_{t}^{\mathrm{ dp}}-\mathbf{G}^{T}\boldsymbol{\lambda}_{t}+\mathbf{f}_{\mu,t} \tag{17}\]  
\[\mathbf{G}\mathbf{u}_{t+\Delta t}-\mathbf{b}=0,\;\mathrm{with} \tag{18}\]  
\[\mathbf{M}=\mathbf{A}_{e=1}^{\mathrm{N_{el}}}\int_{\Omega^{T}}\mathbf{N}_{e} \mathrm{d}\Omega,\mathbf{C}=\alpha\mathbf{M} \tag{19}\]  
\[\mathbf{f}_{t}^{\mathrm{int}}=\mathbf{A}_{e=1}^{\mathrm{N_{el}}}\int_{\Omega^ {T}}\mathbf{g}_{e}^{T}\mathbf{\sigma}\mathrm{d}\Omega \tag{20}\]  
\[\mathbf{f}_{\mu,t}=\sum\
olimits_{md=1}^{\mathrm{N_{cons}}}=\mu F_{N}^{ \mathrm{e}}\frac{\mathbf{g}_{\mathrm{tan},nd}}{||\mathbf{g}_{\mathrm{tan},nd }||},\;\mathbf{g}_{\mathrm{tan}}=\left(\sum\
olimits_{i=1}^{2}\mathbf{t}_{i} \otimes\mathbf{t}_{i}\right)\cdot\mathbf{u}_{\mathrm{rel}} \tag{21}\]  
where the symbols are as follows: \(\mathbf{A}\) is the standard finite element assembly operator; \(\mathbf{M}\) is the mass matrix to be diagonalized by the row sum technique; \(\mathbf{C}\) is the material damping matrix; \(\mathbf{f}_{t}^{\mathrm{int}}\) is the internal nodal force vector, where \(\mathbf{B}_{e}^{T}\) is the transpose of the gradient of the interpolation function matrix \(\mathbf{N}_{e}\); \(\mathbf{f}_{t}^{\mathrm{ext}}\) is the external nodal force vector including the nodal forces from the confining and overburden pressures; \(\mathbf{f}_{t}^{\mathrm{dp}}\) is the damping force due to the viscous damper in Equation (16); \(\mathbf{G}\) is the normal contact constraint matrix, having the normal vectors \(\mathbf{n}_{\mathrm{bit}}=\left[0,-\sin\beta,\;\cos\beta\right]^{\mathrm{T}}\), \(\mathbf{n}_{\mathrm{rock}}=-\mathbf{n}_{\mathrm{bit}}\), of the cutter and the rock contact nodes, as its components; \(\mathbf{b}\) is the initial distance between all the potential contact nodes in the rock and the cutter node; \(\boldsymbol{\lambda}_{t}\) is the Lagrange multiplier vector representing the normal contact force; \(\mathbf{f}_{\mu,t}\) is the friction force vector, calculated as a sum for \(N_{\mathrm{cons}}\) active contact nodes; \(\mu\) is the friction coefficient; \(F_{N}^{\mathrm{e}}\) is the normal contact force for node \(e\) (component \(nd\) of \(\boldsymbol{\lambda}_{t}\)); \(\mathbf{g}_{\mathrm{tan},nd}\) is the relative tangential motion at contact node \(nd\), depending on the unit tangent vectors \(\mathbf{t}_{i}\) and the relative motion vector \(\mathbf{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+\Delta t\), while Equation (17) is written at \(t\).