## 1 Introduction  
The understanding of mixed-mode fracture is an important subject in fracture mechanics, as material flaws or pre-cracks may inevitably occur in the manufacturing process. To describe crack propagation under mixed-mode loading, the classical formula of energy release rate \(G\) was expanded in this hypothesis that crack extends collinearly with the initial crack:  
\[G=\frac{K_{I}^{2}}{E}+\frac{K_{II}^{2}}{E}+\frac{K_{III}^{2}}{2\mu}\ (\text{plane stress}) \tag{1}\]  
where \(K\) is stress intensity factor, \(\mu\) is shear modulus, and for plane strain, elastic modulus \(E\) must be replaced by \(E/(1-v^{2})\), \(v\) is [PERSON]'s ratio. While much research shows that the actual expanding direction is not collinear with its initial path, Formula (1) cannot be applied to a mixed-mode crack [1; 2; 3; 4; 5; 6; 7; 8; 9; 10; 11]. For a I-II mixed-mode crack, the analytic method is widely used in the early development of fracture mechanics; [PERSON]'s [12] work as a founder is also worth mentioning, which simplified the star-shaped crack to an L-shaped crack and used a mapping function to solve the elastic problem. However, some inaccuracy in [PERSON]'s equation can be found as the derivatives of two stress functions in the mapping function during the \"simultaneous expansion\" procedure are incorrect. Meanwhile, complex numerical calculations were used. [PERSON] [13] gave the numerical relationship of \(K-G\) for non-crack-parallel propagations by using explicit asymptotic analysis. [PERSON] and [PERSON] [14] obtainedthe energy release rate in the form of a quadratic of stress intensity factors with the coefficients tabulated for various kink angles by the method that models a kink as a continuous distribution of edge dislocations. [PERSON] [15] revisited energy release rate using the expansion of [PERSON] and [PERSON] [16]. In addition, the validity of [PERSON]'s formula for the energy release rate for any kink angle, material anisotropy and loading condition was proved by [PERSON] and [PERSON] [17]. [PERSON], [PERSON] and [PERSON] [18] used a complex variable method to evaluate the strength of stress singularities at crack tips in plane problems and plate bending problems, which can extend the Griffith-Irwin fracture theory to arbitrary crack extension. [PERSON] and [PERSON] [19] analyzed the stresses of a plate crack in cylindrical coordinates and discussed the experimental results of plane stress and transverse bending of plates, respectively.  
Although many achievements have been made for energy release rate, the validity of conclusions reached through a complex process cannot be evaluated. Most of the researches are aimed at cracks under two kinds of loading situation, the mixed-mode I-II-III loading situation is rarely mentioned, and the methods they used neither analytical method nor numerical method is complex; therefore, they are difficult to generalize to other crack problems.  
The objective of this paper was to evaluate the energy release rate associated with stress intensity factors at any angle under mixed mode loading. In this paper, a more uncomplicated and precise numerical method was established based on the concept that the energy release rate is equal to the change rate of the energy difference before and after crack kinking. A series of spatial inclined ellipses in Mode I-II and ellipsoids in Mode I-II-III with different propagation angles computed from non-dimensional value (\(K/\sqrt{EG}\)) were fitted by MATLAB. Meanwhile, a theoretical expression of energy release rate with angle for a crack tip under I-II-III mixed mode crack was deduced based on the propagation mechanism of the crack tip under the influence of a stress field. It was confirmed that the theoretical expression deduced could provide results as accurately as the present numerical method. The validity of the present methods will be shown by comparing them with experimental data and other literature.