## 2 FEM-CZM Method  
### Separation-Traction Criteria  
In the ABAQUS 2021 software, the cohesive zone model (CZM) adheres to the traction-separation law, with its constitutive model illustrated in Figure 2 Specifically, the model exhibits linear elastic hardening prior to the onset of damage evolution and transitions to linear softening after damage evolution initiates. By introducing a modified parameter to characterize the functional relationship between the external force acting on the cohesive element and the initial damage stress, the model becomes more suitable for analyzing elastoplastic problems [25; 26]. Prior to reaching the damage evolution stage, the external force applied to the cohesive element has not yet reached the initial damage stress, and the relationship between stress and strain follows Hooke's Law. To facilitate parameter calibration, the constitutive thickness of cohesive elements was set to 1 while maintaining the default geometric thickness of 0, thereby ensuring that the nodal separation displacement equals the nominal strain. The linear elastic properties of the cohesive element can be described by the following matrix:  
\[t=\begin{cases}t_{n}\\ t_{s}\\ t_{t}\end{cases}=\frac{I}{L}\begin{bmatrix}E_{nn}&E_{ns}&E_{nt}\\ E_{ns}&E_{ss}&E_{st}\\ E_{nt}&E_{st}&E_{tt}\end{bmatrix}\begin{cases}\delta_{n}\\ \delta_{s}\\ \delta_{t}\end{cases}=\begin{bmatrix}K_{nn}&K_{ns}&K_{nt}\\ K_{sn}&K_{ss}&K_{st}\\ K_{nt}&K_{st}&K_{tt}\end{cases}\begin{cases}\delta_{n}\\ \delta_{s}\\ \delta_{t}\end{cases}. \tag{1}\]  
In the equation, \(t\) represents the traction stress vector; \(t_{n}\) is the normal traction force; \(t_{s}\) and \(t_{t}\) are the two tangential traction forces; \(\delta_{n}\), \(\delta_{s}\), \(\delta_{t}\) are the three components of strain; \(L\) is the initial thickness of the cohesive element; and \(K\) is the stiffness matrix of the cohesive element.  
After reaching the stage of damage evolution, the stiffness of the cohesive elements will linearly decrease. In ABAQUS, a scalar \(D\) is used to represent the degree of damage, with an initial value of 0 and gradually increasing to 1. To accurately describe the damage evolution within the cohesive elements, a function related to the effective displacement \(\delta_{m}\) is introduced as follows:  
\[\delta_{m}=\sqrt{\delta_{n}^{2}+\delta_{s}^{2}+\delta_{t}^{2}}. \tag{2}\]Here, \(\langle\ \rangle\) denotes the Macaulay brackets, which ensure that damage in the cohesive elements only occurs under tensile forces. When \(\langle\delta_{n}\geq 0\rangle\), it indicates tension, and \(\langle\sigma_{n}\rangle\) is equal to \(\sigma_{n}\) itself. When \(\langle\delta_{n}\leq 0\rangle\), it indicates compression, and \(\langle\delta_{n}\rangle\) is equal to zero. The damage variable \(D\) evolves according to the following expression:  
\[D=\frac{\delta_{mf}(\delta_{mm}-\delta_{mo})}{\delta_{mm}(\delta_{mf}-\delta_ {mo})}. \tag{3}\]  
In the context of loading history, let \(\delta_{mm}\) represent the maximum effective displacement, \(\delta_{mo}\) denote the effective displacement at the onset of cracking, and \(\delta_{mf}\) represent the effective displacement at failure. The normal stress during the softening phase can be expressed as follows:  
\[t_{n}=\begin{cases}(1-D)t_{no},&t_{no}\geq 0\\ t_{no},&t_{no}<0\\ t_{s}=(1-D)t_{so}\\ t_{t}=(1-D)t_{to}\end{cases}. \tag{4}\]  
The stiffness of the cohesive element can be expressed using the following formula:  
\[\begin{array}{l}K_{n}=(1-D)K_{no}\\ K_{s}=(1-D)K_{so}\\ K_{t}=(1-D)K_{to}\end{array}. \tag{5}\]  
The change in separation displacement is represented as follows:  
\[\begin{array}{l}\delta_{n}=\frac{t_{n}}{(1-D)}K_{n}\\ \delta_{s}=\frac{t_{s}}{(1-D)}K_{s}\\ \delta_{t}=\frac{t_{t}}{(1-D)}K_{t}\end{array} \tag{6}\]