## 2 Theory and Methods
### Cohesive Element Method  
Since [PERSON] et al. (2016) proposed the cohesive element method in 1976, it has gradually become one of the most widely used crack simulation methods. The cohesive element method simulates cracking by inserting cohesive elements between solid elements. Cohesive elements can effectively simulate the entire crack initiation, propagation, and failure process, thereby reflecting the failure process of quasi-brittle materials (K  
The failure mode of rock typically consists of three zones: the intact rock zone unaffected by cracks, the fracture process zone (FPZ), and the fully cracked zone, as shown in Figure 4. In cohesive elements, failure is manifested through the damage constitutive models. The commonly used constitutive forms of cohesive elements are linear, bilinear, and exponential, as shown in Figure 5. When the load \(\sigma_{m}\) on a node is smaller than the tensile or shear strength \(\sigma_{0}\), the stiffness \(K_{p}\) of the cohesive element node remains unchanged. When the load further increases and exceeds the strength \(\sigma_{0}\), damage begins to occur in the cohesive element, and the stiffness decreases in different forms, indicating different constitutive relationships. The damage \(D_{s}\) is defined as:  
\[D_{s}=1-\frac{K_{m}}{K_{p}} \tag{7}\]  
where \(K_{m}\) represents the nodal stiffness at the nodal load of \(\sigma_{m}\).  
For the linear damage constitutive model, \(D_{s}\) can be expressed as the relative displacement of nodes through transformation:  
\[D_{s}=\frac{\delta_{f}\left(\delta_{m}-\delta_{0}\right)}{\delta_{m}\left( \delta_{f}-\delta_{0}\right)} \tag{8}\]  
where \(\delta_{f}\) represents the maximum relative displacement that cohesive element nodes can undergo, and its value is determined by the area of the constitutive relationship curve, which represents the fracture energy (\(G_{f}\)). From Equation (8), it can be inferred that when \(\delta_{m}=\delta_{0}\) and \(D_{s}=0\), the cohesive element shows no damage; when \(\delta_{m}=\delta_{f}\) and \(D_{s}=1\), the cohesive element is completely damaged.  
Figure 4: Fracture process zone (FPZ).  
Figure 5: Constitutive equation for cohesive elements.  
According to the definition of damage, the constitutive equation for cohesive elements can be written as follows:  
\[\sigma=(1-D_{5})K_{p}\delta \tag{9}\]  
No cracks are generated in the intact rock state, and the constitutive model that represents this state corresponds to the \(\delta_{m}<\delta_{0}\) stage. The displacement \(\delta_{m}\) generated in this state is relatively small and can be considered as the absence of relative displacement between the unit nodes. As the crack further expands, the previously unaffected region enters the vicinity of the crack tip and undergoes plastic-yielding, referred to as the FPZ. The cohesive element simulates this process by incorporating damage treatment to the stiffness through the constitutive equation, and it starts to possess thickness in this state. As the load increases, energy accumulates until it reaches the fracture energy, at which point fracture occurs. The cohesive element has reached its maximum thickness and experiences failure at this stage.