## 2 Fundamental of FDEM
### Contact Model for Triangular Element  
Solid materials come into contact due to the closing of cracks and the sliding or collision of blocks. The generation of contact is dynamic, and it is impossible to know in advance when and where the contact occurs. Therefore, a contact detection algorithm is needed that can dynamically determine which elements are in contact at the current moment during the simulation. Once the elements are in contact, a contact force calculation is performed on these contact elements, and contact force is applied to these contact elements to prevent excessive embedding between them. In the FDEM, the NBS algorithm [42] is used to judge the contact, and the potential contact force algorithm is used to calculate the contact force. This potential contact force algorithm can ensure energy conservation in the contact process, and the contact force is a distributed force rather than a concentrated force.  
The basic idea of the NBS contact detection algorithm is to map the solid element to the background cells, and the size of the cells can just fit the largest solid element. In this way, for a given element, it is only necessary to judge whether the element in the cell is in contact with the element in this cell and adjacent cells. In this way, the number of elements involved in the detection will be greatly reduced. Thus, the efficiency of contact detection is improved.  
The calculation of the contact force between the blocks or the media on both sides of the fracture is transformed into the calculation of the contact force between the triangular elements in the FDEM. When judging the contact of triangular elements, the interaction forces with equal force but opposite directions are applied to them to prevent the two triangular elements from being embedded in each other. The contact force includes the normal contact force and the tangential force. The following will focus on the normal contact force and the tangential contact force.  
(i) Normal Contact Force  
As shown in Figure 3, the two contact triangular elements are called the active triangle \(\beta_{\rm c}\) and the passive triangle \(\beta_{\rm t}\), respectively. A closed area will be formed when they are in  
Figure 2: Fracture constitutive model of joint elements: (**a**) tensile failure; (**b**) shear failure; (**c**) mixed tensile and shear failure [11; 12; 20].  
contact. The normal contact force between these two triangular elements can be expressed as the area integral of potential function gradient in the overlapping area as:  
\[f_{\rm c}=p_{\rm n}\int\limits_{\Gamma_{\beta_{\rm tr}\cap\beta_{\rm c}}}[{\bf grad }\varphi_{\rm c}(P_{\rm c})-{\bf grad}\varphi_{\rm t}(P_{\rm t})]{\rm d}A \tag{10}\]  
Applying [PERSON]'s formula to Equation (10), the area integral can be transformed into a line integral of the boundary of the overlapping area:  
\[f_{\rm c}=p_{\rm n}\oint\limits_{\Gamma_{\beta_{\rm tr}\cap\beta_{\rm c}}}{\bf n }_{\Gamma}(\varphi_{\rm c}-\varphi_{\rm t}){\rm d}\Gamma \tag{11}\]  
where \(\Gamma_{\beta_{\rm tr}\cap\beta_{\rm c}}\) is the boundary of the overlapping area of the two triangular elements, \({\bf n}_{\Gamma}\) is the outer normal direction vector of the boundary of the overlapping area.  
(ii) Tangential contact force  
After obtaining the total normal contact force and the force action point, the tangential friction force can be calculated by:  
\[{\bf F}_{\rm t}^{t+\Delta t}={\bf F}_{\rm t}^{t}-p_{\rm t}\Delta{\bf u}_{s} \tag{12}\]  
where the \(p_{\rm t}\) is the tangential penalty parameter and the \(\Delta{\bf u}_{s}\) is the relative displacement increment at the current time step. When the tangential contact force calculated by the equation satisfies \(\left|{\bf F}_{\rm t}^{t+\Delta t}\right|\geq F_{n}\tan\phi\), the tangential friction force can be calculated by:  
\[{\bf F}_{t}^{t+\Delta t}=\frac{{\bf F}_{t}^{t+\Delta t}}{\left|{\bf F}_{t}^{t+ \Delta t}\right|}u|{\bf F}_{n}| \tag{13}\]  
where \(u\) is the coefficient of friction.