# Bonded-Particle Model with Nonlinear Elastic Tensile Stiffness for Rock-Like Materials
## 2 Method Description
### Numerical Simulation Model  
The disc-shaped particles are considered as rigid bodies and their translational motions and rotational motions are governed by the standard equations of rigid body dynamics:  
\[\mathbf{F}_{p}+\mathbf{F}_{p}^{ext}+\mathbf{F}_{p}^{d}=m_{p}\ddot{\mathbf{u}} _{p} \tag{10}\]  
\[M_{p}+M_{p}^{ext}+M_{p}^{d}=J_{p}\ddot{\rho}_{p} \tag{11}\]  
where \(\mathbf{F}_{p}^{ext}\) and \(M_{p}^{ext}\) are the external force vector and external moment applied to particle \(p\) respectively. \(\mathbf{\dot{}}_{p}^{d}\) and \(M_{p}^{d}\) are damping force and damping moment respectively, which are introduced to dissipate kinetic energy and gain a quasi-static state of equilibrium of the particles and are expressed as [34]:  
\[\mathbf{F}_{p}^{d}=-\alpha^{t}\|\mathbf{F}_{p}+\mathbf{F}_{p}^{ext}\|\frac{ \dot{\mathbf{u}}_{p}}{\|\dot{\mathbf{u}}_{p}\|} \tag{12}\]  
\[M_{p}^{d}=-\alpha^{r}\Bigl{|}M_{p}+M_{p}^{ext}\Bigl{|}\frac{\dot{\theta}_{p}}{ \Bigl{|}\dot{\theta}_{p}\Bigr{|}} \tag{13}\]  
where \(\alpha^{t}\) and \(\alpha^{r}\) are translational and rotational damping constant respectively. The centroid displacement and rotation angle of particle \(p\) can be obtained by using an explicit central difference method to integrate Equations (10) and (11). The algorithm is described as follows:  
\[\dot{\mathbf{u}}_{p}^{i+0.5}=\dot{\mathbf{u}}_{p}^{i-0.5}+\ddot{\mathbf{u}}_{ p}^{i}\Delta t \tag{14}\]  
\[\mathbf{u}_{p}^{i+1}=\mathbf{u}_{p}^{i}+\ddot{\mathbf{u}}_{p}^{i+0.5}\Delta t \tag{15}\]  
Figure 4: Shear strength criterion for the NET-BPM (nonlinear elastic tensile bonded particle model). \[\dot{\theta}_{p}^{i+0.5}=\dot{\theta}_{p}^{i-0.5}+\ddot{\theta}_{p}^{i}\Delta t \tag{16}\]  
\[\dot{\theta}_{p}^{i+1}=\theta_{p}^{i}+\dot{\theta}_{p}^{i+0.5}\Delta t \tag{17}\]  
where the superscript \(i\) indicates \(i\)th step. \(\Delta t\) should not exceed the critical time step \(\Delta t_{c}\) to remain numerical stability as introduced in [34]. The forces and moments are transmitted through the contacts between particles or the contacts between particles and plates. There are two kinds of particle contacts according to whether an effective bond exists between the contact particles: the bonded particle contacts and the unbonded particle contacts. Two particles \(p\) and \(q\) are bonded together in the original rock specimens if their relative position satisfies the Equation (18):  
\[\|\mathbf{u}_{p}-\mathbf{u}_{q}\|\leq(1+\delta)(r_{p}+r_{q}) \tag{18}\]  
where \(\delta\) is a nonnegative constant defined by user to eliminate the numerical errors caused during the particle generating. The errors may result in separation of two particles which should touch each other. A bond is created between every bonded particles. If \(\delta\) is too large, it would lead to forming bond between two particles far away from each other. In our model, \(\delta\) = 0.001 is suggested. When the bond is broken, the bonded particle contact becomes unbonded particle contact. Furthermore, two particles \(p\) and \(q\) which are not close to each other may approach each other during the simulation process, and form an unbonded particle contact when their relative position satisfies the Equation (19):  
\[M_{p}+M_{p}^{ext}+M_{p}^{d}=J_{p}\ddot{\rho}_{p} \tag{11}\]0  
The forces and moments on the bonded particle contacts are calculated by Equations (1)-(3). While forces on the unbonded particle contacts, where no moment exists are calculated by Equations (3) and (7). The bonded particle contacts are identified at the beginning of the simulation and updated in each time steps if any bond is broken. The unbonded contacts should be identified in every time steps. In order to improve the computational efficiency, background grids and cells are introduced in searching particle contacts. In each step, all the particles are stored in the cells formed by the grids according to their centroid position coordinates.