## 1 Introduction  
Composite materials are now widely used in many fields, such as aircraft fuselages and spacecraft propellant tanks, due to their ability to withstand harsh conditions. At present, numerous composite materials and structures with heterogeneous properties have been successfully developed and manufactured. Composite materials usually contain complex interfaces, which serve as nucleation sites for cracks, voids, and other defects. Therefore, in comparison to homogeneous materials, composite materials are more susceptible to fracture and failure [1]. A crucial aspect of analyzing crack propagation issues is the effective determination of crack-tip parameters. Crack tip fracture mechanics' characteristic parameters mainly include stress intensity factors (SIFs), strain energy release rate (G), and crack tip opening displacement (COD). Among them, the stress intensity factor, the most widely used fracture mechanics parameter [2], is of great value in characterizing the stress, strain, and displacement fields near the crack tip. Methods for solving stress intensity factors can be divided into direct and indirect methods. Direct methods include displacement and stress methods, which calculate results via interpolation formulas but suffer from limited accuracy. In contrast, the J-integral method, a path-independent contour integral method, has the advantage of obtaining integral values independent of path selection, making it especially suitable for solving stress intensity factors at crack tips in heterogeneous materials. Taking the crack problem of two-dimensional linear elastic materials as an example, its J-integral can be defined as follows:  
\[J\ =\ \lim_{\Gamma\to 0}\int_{\Gamma}\ (W\delta_{1i}-\sigma_{ij}u_{j,1})n_{i}d\Gamma, \tag{1}\]  
In this equation, \(W\) represents the strain energy density. \[W\ =\ \frac{1}{2}\sigma_{ij}\varepsilon_{ij}\ =\ \frac{1}{2}C_{ijkl}\varepsilon_{ij} \varepsilon_{kl}\ =\ \frac{1}{2}S_{ijkl}\sigma_{ij}\sigma_{kl}, \tag{2}\]  
In this equation, \(\delta_{ij}\) is the Kronecker delta, \(C_{ijkl}\) and \(S_{ijkl}\) are the stiffness and compliance tensors, respectively, and \(n_{i}\) is the unit outward normal vector to the integration path \(\Gamma\), as shown in Figure 1.  
Nevertheless, it is difficult to separate Mode I and Mode II stress intensity factors in mixed-mode fracture problems using the J-integral method. In the development of the J-integral method, the interaction integral method has shown significant advantages in calculating the mixed-mode stress intensity factors. By superimposing the actual and auxiliary fields, the mixed-mode stress intensity factors can be efficiently separated [3]. After years of theoretical refinement and algorithmic optimization by researchers, this method has proven to provide a advantage for the fracture analysis of linear elastic materials with complex interfaces [4; 5]. Figure 2 illustrates a simulation of crack propagation in complex interfacial materials using the interaction integral method. [PERSON] et al. [6] used the interaction integral method to calculate the stress intensity factors of mixed-mode cracks in anisotropic homogeneous materials. Subsequently, [PERSON] [7] established a three-dimensional integral solution method for the stress intensity factors of bimaterial interface cracks, enabling the precise extraction of crack tip fracture parameters through a suitable auxiliary field. [PERSON] et al. [8] proposed an unstructured tetrahedral element interaction integral method, addressing the path independence issue in three-dimensional functionally graded material crack analysis. [PERSON] and [PERSON] [9] successfully applied the interaction integral method to axisymmetric interface crack analysis, while [PERSON] et al. [10] applied it to three-dimensional interface cracks with curved crack fronts.