## 2 Theoretical Model  
### Bilinear CZM  
Based on the bilinear CZM, this study established a cyclic loading CZM. The bilinear CZM is shown in Figure 8.  
When the CZM reached the critical stress \(\sigma_{\text{max}}\) under loading, damage began to occur in the material layer. The crack propagation in the brazed joint manifests as a mixed-mode cracking pattern. The square stress criterion was adopted with the following formula [18]:  
\[\left\{\frac{<\sigma_{n}>}{\sigma_{nc}}\right\}^{2}+\left\{\frac{\sigma_{t}}{ \sigma_{tc}}\right\}^{2}+\left\{\frac{\sigma_{s}}{\sigma_{sc}}\right\}^{2}=1 \tag{1}\]  
Figure 8: Bilinear CZM.  
Figure 7: Fracture morphology of butt brazing joints at different angles.  
(1) Static loading-unloading phase  
The stiffness K in the elastic deformation stage remained unchanged. During the static damage stage, when the loading displacement exceeds the critical damage displacement, damage parameter D is introduced, and the traction force formula is expressed as follows:  
\[\left[\begin{array}{c}\sigma_{n}\\ \tau_{s}\\ \tau_{t}\end{array}\right]=\left[\begin{array}{cc}K_{n}&&\\ &K_{s}&\\ &&K_{t}\end{array}\right](1-\text{D})\left[\begin{array}{c}\delta_{n}\\ \delta_{s}\\ \delta_{t}\end{array}\right] \tag{4}\]  
where \(\sigma\) and \(\tau\) represent the normal and tangential stresses, respectively (MPa); \(\delta_{n}\) represents the normal displacement, and \(\delta_{s}\) and \(\delta_{t}\) represent the two tangential displacements (mm).  
Based on Equations (2)-(3), the static damage \(D_{m}\) under the mixed loading mode in the normal and tangential directions can be transformed into the following formula:  
\[D_{m}=\frac{\delta_{mf}(\delta_{m}-\delta_{mc})}{\delta_{m}(\delta_{mf}- \delta_{mc})} \tag{5}\]  
where \(\delta_{mc}\) represents the critical displacement of damage under mixed loading (mm), and 2 represents the failure displacement for damage under mixed loading (mm).  
(2) Cyclic loading-unloading phase  
After the first static loading-unloading cycle, the fatigue equation must be determined based on the static damage. In this section, the fatigue damage equation proposed by [PERSON] [21] is adopted. When considering mixed loading, the fatigue equation is transformed into  
\[\overset{\bullet}{D}_{mc}=A\frac{\left|\Delta\overset{\bullet}{u}\right|}{B \delta_{mc}}\left[\frac{\sigma_{m}}{\sigma_{mc}}-C_{mf}\right]H(\Delta u- \delta_{mc}) \tag{6}\]  
where \(\overset{\bullet}{D}_{mc}\) represents the fatigue damage corresponding to the strain value during each loading; \(\Delta\overset{\bullet}{u}\) represents the strain increment; \(\Delta u\) represents the current strain value; \(\delta_{mc}\) represents the critical displacement for damage; \(\sigma_{m}\) represents the current stress, and \(\sigma_{mc}\) represents the current critical stress for damage. A, B, and \(C_{mf}\) represent fatigue damage parameters.  
The total damage accumulation in the cyclic CZM is the sum of static and cyclic damages. The total damage \(D_{t}\) under the mixed loading mode of normal and tangential directions is obtained as  
\[D_{t}=\int{(D_{m}+D_{mc}^{\bullet})dt} \tag{7}\]  
Through the accumulation of damage in these two stages, when the cyclic damage of the material reached a maximum value of 1, the stiffness decreased to 0, indicating crack propagation and ultimate failure.  
Based on the above theoretical model, the VUMAT subroutine is compiled to read the basic parameters of the cohesive zone model in the ABAQUS software. Combined with the subroutine, the failure calculation of cohesive elements is performed. The calculation steps are shown in Figure 10.