Crack Growth and Energy Release Rate for an Angled Crack under Mixed Mode Loading

Featured Application

Damage fracture for engineering materials.

Abstract

The evaluation of energy release rate with angle is still a challenging task in metal crack propagation analysis, especially for the mixed Mode I-II-III loading situation. In this paper, the energy release rate associated with stress intensity factors at an arbitrary angle under mixed mode loadings has been investigated using both a numerical method and theoretical derivation. A relatively simple and precise numerical method was established through a series of spatial-inclined ellipses in Mode I-II and ellipsoids in Mode I-II-III, with different propagation angles computed from simulation. Meanwhile, a theoretical expression of the energy release rate with angle for a crack tip under a 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 is confirmed that the theoretical expression deduced could provide results as accurately as the present numerical method. The present results were confirmed to be effective and accurate by comparison with experimental data and other literature.

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 } (\mathrm{plane} \mathrm{stress})$$

(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 Poisson’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; Hussian’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 Hussain’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. Wu, C. H [13] gave the numerical relationship of $K-G$ for non-crack-parallel propagations by using explicit asymptotic analysis. K. Hayashi and S. Nemat-Nasser [14] obtained the 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. Antonin Chambolle [15] revisited energy release rate using the expansion of Amestoy and Leblond [16]. In addition, the validity of Irwin’s formula for the energy release rate for any kink angle, material anisotropy and loading condition was proved by Abbas Azhdari and Sia Nemat-Nasser [17]. G.C. Sih, P. C. Paris and F. Erdogan [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. F. Erdogan and G.C. Sih [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.

2. Determining G for I-II Mixed Mode Crack

2.1. Numerical Analysis of 2D Crack

2.1.1. Crack Simulation with ABAQUS

To find the relationship between the stress intensity factors and energy release rate under the I-II combined loading situation, finite element technology was applied to research the problem.

As shown in Figure 1a, a rectangular plate with crack was under combined Mode I and Mode II loading on the top side, and fixed on the bottom side, $L=W=50 \mathrm{mm}$, $A=10 \mathrm{mm}$. the finite element model was shown in Figure 1b,c shows the meshes around the crack tip, the type of which is a collapsed element side, and duplicate nodes configure the mid-side node parameter to its 1/4 position to simulate the singularity of displacement of the crack tip region. The material is linear elastic homogeneous isotropic aluminum alloy, $E=70000 \mathrm{Mpa}, v=0.33$. The $K$ and $G$ are obtained by changing the value of $\sigma :\tau$.

When $\Delta a=0 \mathrm{mm}, \theta =0^{°}, \tau =0 \mathrm{Mpa}$, the value of $\sigma$ changes from 20 $\mathrm{Mpa}$ to 100 $\mathrm{Mpa}$, with a step of 20. The result is shown in Figure 2. The simulation value output from ABAQUS is basically consistent with the reference value calculated by the literature [20], and the errors between the two values of each loading presented in Figure 2b are about 0.53%. This verified simulation method will be applied in the following sections.

2.1.2. The Theoretical Basis of Energy Release Rate Definition

A cracked plate in plane stress has an infinitesimal kink at an angle $\theta$ from the plane of the crack, as illustrated in Figure 3a–c. The energy release rate $G$, as a crack extension force, is a measure of the energy available for an increment of crack extension, which can be calculated approximately as:

$$G=-\frac{dU}{da}=-{\displaystyle \underset{a_2\to a_1}{lim}}\frac{\left(-U_2\right)-\left(-U_1\right)}{a_2-a_1}\approx \frac{U_2-U_1}{\Delta a}$$

(2)

Irrespective of whether the crack extension is in its initial plane or kinked, the angle is $\theta$, where $U_1$ and $U_2$ are the strain energy before and after kink extension, respectively, and $a$ is the length of the kink, which is measured from the original crack tip.

For an edge crack in the linear elastic body, the J-integral computed along a contour surrounding the crack tip is equal to the energy release rate (Figure 3d), which is revalidated by simulation work. The strain energy release rate and $J$-integral for a kink can be evaluated based on the output data of ABAQUS, as shown in Figure 4. It shows that energy release rate is in good agreement with the $J$-integral evaluated along the contour starting at one kink’s crack surface and ending at the other kink’s crack surface under arbitrary loading and kink angle conditions.

(1)

Obtain the values of $J$ and strain energy at the points $a_0=0 \mathrm{mm}, a_1=0.1 \mathrm{mm}, a_2=0.3 \mathrm{mm}, a_3=0.5 \mathrm{mm}, a_4=1 \mathrm{mm}.$

(2)

The energy differences between $a_0 \mathrm{and} {\mathrm{a}}_1, {\mathrm{a}}_0 \mathrm{and} {\mathrm{a}}_2, {\mathrm{a}}_0 \mathrm{and} {\mathrm{a}}_3, {\mathrm{a}}_0 \mathrm{and} {\mathrm{a}}_4$ are calculated and divided by $a_1-a_0=0.1 \mathrm{mm}, a_2-a_0=0.3 \mathrm{mm}, a_3-a_0=0.5 \mathrm{mm}, a_4-a_0=1 \mathrm{mm}.$ It represents the energy release rate for the crack at 0.05 mm, 0.15 mm, 0.25 mm, 0.5 mm.

(3)

Determine that the minimum value of kink length for calculation is 0.1 mm; since the kink crack becomes infinitesimally small, the values from the energy difference deviate from the $J$ integral.

This validates the assumption made in the connection with strain energy release rate. $J$-integral will be used in the following research as it is reliable and convenient.

2.1.3. Computational Analysis

From the output data of the simulation, the energy release rate and stress intensity factors can be obtained corresponding to the required condition. The method developed in the present study is described as follows:

(1)

Determine $K_I$ and $K_{II}$ under arbitrary combined Mode $\mathrm{I}$ and Mode $\mathrm{I}\mathrm{I}$ loading conditions with the initial crack.

(2)

Determine $J$ for kink cracks repeatedly, varying the kink length from $0.1 \mathrm{mm}$ to $1 \mathrm{mm}$ with four data points $\left(a=0.1 \mathrm{mm}, 0.3 \mathrm{mm}, 0.5 \mathrm{mm}, 1 \mathrm{mm}\right)$. Then, compute energy release rate as the kink propagation vanishes on the curve of $J$ versus kink length by the method of fitting.

(3)

Plot the non-dimensional value ($K_I/\sqrt{EG}, K_{II}/\sqrt{EG}$) for each kink angle ($\theta ={30}^{°}, {60}^{°}, {90}^{°},{120}^{°}, {150}^{°}$)

(4)

Determine the parameters of inclined ellipses that fit the above data points using the curve fitting algorithm by MATLAB, and the corresponding ellipse equation is presented as follows:

$$\frac{{(\frac{K_I}{\sqrt{EG}}cos\alpha +\frac{K_{II}}{\sqrt{EG}}sin\alpha )}^2}{a^2}+\frac{{(\frac{K_{II}}{\sqrt{EG}}cos\alpha -\frac{K_I}{\sqrt{EG}}sin\alpha )}^2}{b^2}=1$$

(3)

where $a$, $b$, $\alpha$ are the semi-major axis, semi-minor axis and inclination of ellipse, respectively.

(5)

Obtain the coefficients of quadratic of energy release rate in terms of stress intensity factors, and defined as

$$G=\frac{1}{E}\left(a_{11}{K}_I^2+a_{12}{K}_I{K}_{II}+a_{22}{K}_{II}^2\right)$$

(4)

As shown in Figure 5, the non-dimensional value ($K_I/\sqrt{EG}, K_{II}/\sqrt{EG}$) for each kink angle can be calculated through the technique proposed above. By fitting the database, a series of inclined ellipses are presented with a certain angle and size. Figure 6 shows the inclined angle $\alpha$ of Iso-energy release rate ellipses with kink angle $\theta$. The variation of inclination for ellipses drops linearly with the increase of the kink angle. The trend of semi-major axis and semi-minor axis of ellipses according to the variation of kink angle has been shown in Figure 7. The coefficients of quadratics for energy release rate in terms of stress intensity factors are calculated for each kink angle based on the ellipse equations (Figure 8).

2.2. Discussion of Numerical Results

2.2.1. Comparison with References

The results obtained by the present method are more evident through the comparison with Anderson [21]. Anderson’s analytic work introduced the local SIFs, which are functions of nominal SIFs—the local SIFs represent the true stress strength in front of a kink crack. Figure 9a shows the coefficients calculated by Anderson’s theory and the present method.

The values of coefficients $a_{11},a_{12}$ and $a_{22}$ calculated by Anderson and the present method are surprisingly consistent. Figure 9b presents the errors between the two results for each kink angle; the average error between each kink angle is less than 1%, which illustrates that the results obtained by the present work agree well with Anderson’s.

Moreover, the validity of extended Irwin’s formula is reiterated by the numerical simulation under two arbitrary loading conditions, as shown in

Table 1. Comparison of G for two different methods.

$\mathsf{\sigma }=-\mathbf{10.4}\times {\mathbf{10}}^{\mathbf{5}} MPa, \mathit{\tau }=\mathbf{2.5}\times {\mathbf{10}}^{\mathbf{5}} MPa$
$\mathsf{\theta }$ (deg)	kI(θ) from ABAQUS ($MPa\sqrt{mm}$)	kII(θ) from ABAQUS $MPa\sqrt{mm}$	G(θ) for Calculation (mJ)	G for Simulation (mJ)
30°	1.31 × 105	9.90 × 105	4.84	4.84
60°	−5.01 × 105	6.26 × 105	3.12	3.17
90°	−7.72 × 105	1.30 × 105	2.98	2.98
120°	−6.88 × 105	−2.42 × 105	2.58	2.58
150°	−4.02 × 105	−3.22 × 105	1.29	1.29
180°	0.00 × 100	0.00 × 100	0.00	0.00
$\mathsf{\sigma }=\mathbf{0} MPa, \mathit{\tau }={\mathbf{10}}^{\mathbf{5}} MPa$
$\mathsf{\theta }$ (deg)	kI(θ) from ABAQUS $MPa\sqrt{mm}$	kII(θ) from ABAQUS $MPa\sqrt{mm}$	G for Calculation (mJ)	G for Simulation (mJ)
30°	2.91 × 106	1.18 × 106	47.82	47.82
60°	1.85 × 106	1.46 × 106	26.84	26.84
90°	8.14 × 105	1.18 × 106	10.02	10.01
120°	1.50 × 105	5.97 × 105	1.84	1.84
150°	−8.86 × 104	1.02 × 105	0.09	0.09
180°	0.00 × 100	0.00 × 100	0.00	0.00

. Stress intensity factors $k_I\left(\theta \right) \mathrm{and} k_{II}\left(\theta \right)$ and energy release rate $G$ with the limiting process as the propagation kink goes to zero are obtained from the output of the ABAQUS. The calculated value of $G\left(\theta \right)$ obtained from Equation (3) is very consistent with the value from ABAQUS.

2.2.2. Experiment Verification

To research the law of crack propagation under combined loading, the uniaxial oblique crack experiment was carried out. For an oblique crack under a tensile loading shown in Figure 10, the stress intensity factors can be calculated by Formula (5):

$$\begin{array}{c}{K}_I=\sigma \sqrt{\pi a}{\sin }^2\beta \\ K_{II}=\sigma \sqrt{\pi a}\sin \beta \cos \beta \end{array}\}$$

(5)

Reference [22] verified the validity of specimen for mixed mode fracture experiment. The composition of the 6061-T6 aluminum alloy which is used to manufacture test specimens is shown in

Table 2. Composition of 6061-T6 aluminum alloy.

Chemical Composition	Cu	Si	Fe	Mn	Mg	Zn	Cr	Ti	other	AL
Ratio	25%	60%	70%	15%	85%	25%	16%	15%	15%	margin

. Three groups were divided by inclination angle (30°, 60°, 90°), and each group has three samples, as shown in Figure 11a. Samples are processed by middle-speed WEDM, and the accuracy is 0.2 mm. Figure 11b shows the Universal Hydraulic Testing Machine, type LD26.105. The samples with oblique crack were of the exerted displacement boundary condition on the top end at a speed of 0.2 mm/s, and fixed on the bottom end until the crack started propagation.

As shown in Figure 12a, it was revealed that when the inclination angle was 90°(Mode I crack), the crack extended in a self-similar manner that propagated along the direction of the initial crack. While the phenomenon was unsuitable to describe 30° and 60° (mixed I-II crack) extension situations, crack propagation under these two angles deviated from the initial direction, which illustrated that the assumption for the mixed-mode crack extension path of Formula (1) was wrong.

From the experimental results, the actual extension angles were measured as shown in Figure 12b; besides, the value $K_I/K_{II}$, calculated by expression (5) for a inclined crack, was plugged into the arbitrary angle energy release rate formula $G\left(\theta \right)$ obtained by numerical fitting. Therefore, the extension angles were deduced by the $MG$-criterion, and

Table 3. Calculated extension angle and experimental extension angle.

Crack Angle (°)	Sample Number	Calculated Extension Angle (°)	Experimental Extension Angle (°)	Experimental Extension Angle (Average) (°)	Error (°)
30	1	−60	−59	−58.70	2.17
	2	−60	−59.5
	3	−60	−57.6
60	1	−43	−41.5	−41.17	4.26
	2	−43	−40.8
	3	−43	−41.2
90	1	0	0	0	0
	2	0	0
	3	0	0

shows the results:

The results between the calculated extension angle and experimental extension angle are extraordinary similar. The error value, considering the specimens are not manufactured by ideal elastic material, was lower than the 5% that can be accepted. Consequently, the validity of the numerical method can be verified by the two comparison results above, which can be applied for the I-II-III mixed mode crack.

3. Fracture Parameters for 3D Crack

3.1. The Theoretical Derivation of Energy Release Rate for I-II-III Mixed Mode Crack

Considering the meaning of SIFs, $K_{I,} K_{II,} K_{III}$ denote the intensity of tensile stress, in-plane stress and out-of-plane stress, as shown in Figure 13. When crack kink occurs, the local SIFs should be recalculated. A local $x^{\prime }-y^{\prime }-z^{\prime }$ coordinate system at the tip of the kink is defined and assumes that the remote stress fields in polar coordinates defined the local stress field. The transformation of the $x^{\prime }-y^{\prime }$ plane is shown in Figure 14.

Stress fields ahead of the crack front for Mode I, Mode II and Mode III at an angle $\theta$ in polar coordinates are shown in Formula (6).

$$\begin{array}{c}{\sigma }_{rr}=\frac{K_I}{\sqrt{2\pi r}}\left[\frac{5}{4}cos\left(\frac{\theta }{2}\right)-\frac{1}{4}cos\left(\frac{3\theta }{2}\right)\right]+\frac{K_{II}}{\sqrt{2\pi r}}\left[-\frac{5}{4}sin\left(\frac{\theta }{2}\right)+\frac{3}{4}sin\left(\frac{3\theta }{2}\right)\right]\\ {\sigma }_{\theta \theta }=\frac{K_I}{\sqrt{2\pi r}}\left[\frac{3}{4}cos\left(\frac{\theta }{2}\right)+\frac{1}{4}cos\left(\frac{3\theta }{2}\right)\right]+\frac{K_{II}}{\sqrt{2\pi r}}\left[-\frac{3}{4}sin\left(\frac{\theta }{2}\right)-\frac{3}{4}sin\left(\frac{3\theta }{2}\right)\right]\\ {\tau }_{r\theta }=\frac{K_I}{\sqrt{2\pi r}}\left[\frac{1}{4}sin\left(\frac{\theta }{2}\right)+\frac{1}{4}sin\left(\frac{3\theta }{2}\right)\right]+\frac{K_{II}}{\sqrt{2\pi r}}\left[\frac{1}{4}cos\left(\frac{\theta }{2}\right)+\frac{3}{4}cos\left(\frac{3\theta }{2}\right)\right]\\ {\sigma }_z=v\left({\sigma }_{rr}+{\sigma }_{\theta \theta }\right)\\ {\tau }_{rz}=\frac{K_{III}}{\sqrt{2\pi r}}sin\left(\frac{\theta }{2}\right)\\ {\tau }_{\theta z}=\frac{K_{III}}{\sqrt{2\pi r}}cos\left(\frac{\theta }{2}\right)\end{array}\}$$

(6)

The directions of the local stress field for the crack tip in a local $x^{\prime }-y^{\prime }-z^{\prime }$ coordinate system are consistent with those of the polar coordinate, as shown in Figure 14. Thus, the stress at the crack tip between the polar coordinate and local $x^{\prime }-y^{\prime }-z^{\prime }$ coordinate has the following relationship:

$$\begin{array}{c}{\sigma }_{y^{\prime }{y}^{\prime }}={\sigma }_{\theta \theta }\\ {\tau }_{x^{\prime }{y}^{\prime }}={\tau }_{r\theta }\\ {\tau }_{y^{\prime }{z}^{\prime }}={\tau }_{\theta z}\end{array}\}$$

(7)

The local stress intensity factors at the arbitrary extend angle $\alpha$ were deduced by summing the normal and shear stresses, respectively:

$$\begin{array}{c}{k}_I^{\prime }\left(\alpha \right)={\sigma }_{y^{\prime }{y}^{\prime }}\sqrt{2\pi r}=C_{11}{K}_I+C_{12}{K}_{II}\\ k_{II}^{\prime }\left(\alpha \right)={\tau }_{x^{\prime }{y}^{\prime }}\sqrt{2\pi r}=C_{21}{K}_I+C_{22}{K}_{II}\\ k_{III}^{\prime }\left(\alpha \right)={\tau }_{y^{\prime }{Z}^{\prime }}\sqrt{2\pi r}=C_{33}{K}_{III}\end{array}\}$$

(8)

Combining Formula (6)–(8), the coefficients can be calculated using Formula (9):

$$\begin{array}{c}{C}_{11}=\frac{3}{4}\cos \left(\frac{\alpha }{2}\right)+\frac{1}{4}\cos \left(\frac{3\alpha }{2}\right)\\ C_{12}=-\frac{3}{4}sin\left(\frac{\alpha }{2}\right)-\frac{3}{4}sin\left(\frac{3\alpha }{2}\right)\\ C_{21}=\frac{1}{4}sin\left(\frac{\alpha }{2}\right)+\frac{1}{4}sin\left(\frac{3\alpha }{2}\right)\\ C_{22}=\frac{1}{4}\cos \left(\frac{\alpha }{2}\right)+\frac{3}{4}\cos \left(\frac{3\alpha }{2}\right)\\ C_{33}=\cos \left(\frac{\alpha }{2}\right)\end{array}\}$$

(9)

Reviewing the expression of the classical energy release rate, the new expression for the I-II-III mixed mode crack can be modified as (10):

$$G\left(\alpha \right)=\frac{k_I^{\prime }\left(\alpha \right)}{E}+\frac{k_{II}^{\prime }\left(\alpha \right)}{E}+\frac{k_{III}^{\prime }\left(\alpha \right)}{2\mu }$$

(10)

3.2. Numerical Analysis of 3D Crack

3.2.1. 3D Model of Crack Simulation

A 3D model with a mixed mode I-II-III crack is considered here: $L=50 \mathrm{mm}, W=50 \mathrm{mm}, D=20 \mathrm{mm}, A=10 \mathrm{mm}$; the bottom is fixed, and normal stress and shear stress are applied at the top. The 20-nodes collapsed element used in the crack tip region is used to simulate the singularity of the displacement at the crack tip (Figure 15).

3.2.2. Computational Analysis

The data processing for $K$ and $G$ outputted by the I-II-III mixed mode crack simulation is more complex than for the I-II mixed mode crack—not only for the value $K_{III}$ added, but also for much more data on the crack tip line. To obtain the fracture parameter $K$, the present study is described as follows:

(1)

Determine $K$ and $J$ for kink cracks repeatedly, varying the kink length from $0.1 \mathrm{mm}$ to $1 \mathrm{mm}$ with four data points $\left(a=0.1 \mathrm{mm}, 0.3 \mathrm{mm}, 0.5 \mathrm{mm}, 1 \mathrm{mm}\right)$. Then, compute the energy release rate as the kink propagation vanishes on the curve of $K$ and $J$ versus kink length through fitting each node on the crack tip line.

(2)

Plot the non-dimensional value ($K_I/\sqrt{EG}, K_{II}/\sqrt{EG},K_{III}/\sqrt{EG}$) for each kink angle ($\theta ={30}^{°}, {60}^{°}, {90}^{°},{120}^{°}, {150}^{°}$)

(3)

Determine the parameters of inclined ellipses that fit the above data points using a fitting method for quadric surfaces in space by MATLAB, and the corresponding ellipse equation is presented as follows:

$$\frac{{(\frac{K_I}{\sqrt{EG}}cos\alpha +\frac{K_{II}}{\sqrt{EG}}sin\alpha )}^2}{a^2}+\frac{{(\frac{K_{II}}{\sqrt{EG}}cos\alpha -\frac{K_I}{\sqrt{EG}}sin\alpha )}^2}{b^2}+\frac{{(\frac{K_{III}}{\sqrt{2\mu G}})}^2}{C^2}=1$$

(11)

where $a$, $b$, $c$, $\alpha$ are the semi-major axis, semi-middle axis, semi-minor axis and inclination of ellipse, respectively.

(4)

Obtain the coefficients of quadratic of energy release rate in terms of stress intensity factors, and defined as

$$G=\frac{1}{E}\left[a_{11}{K}_I^2+a_{12}{K}_I{K}_{II}+a_{22}{K}_{II}^2+a_{33}\left(1+v\right)K_{III}^2\right]$$

(12)

The results of fitting the non-dimensional value are shown in Figure 16; a series of conclusions can be summarized as follows:

(1)

The non-dimensional value ($K_I/\sqrt{EG}, K_{II}/\sqrt{EG},K_{III}/\sqrt{EG}$ ) transformed from the data output by ABAQUS forms a series of ellipsoids.

(2)

At each presupposed extended angle (30°~150°), the scattered points of dimensionless data output by ABAQUS under different loadings are all on the surface of the same tilted ellipsoid.

(3)

The coefficients of quadratics for energy release rate in terms of stress intensity factors are calculated for each kink angle based on the ellipsoid equations.

3.3. Discussion of Theoretical and Numerical Results

This section expanded the I-II mixed mode crack to the I-II-III mixed mode crack, not only added out-of-plane loading, but also considered the thickness factor. Therefore, the nodes along the thickness direction in a finite element model can output multi group values of $K$ and $G$. Based on the result of (12), an arbitrary loading had been imposed on the crack model, and the output data of $K$ and $G$ were obtained. $G\left(\theta \right)$ is calculated with the obtained $K$ from the theoretical expression (10). Figure 17 shows the $G\left(\theta \right)/G$ along the thickness direction, the values at node 4~17 are at about a constant value of 1.1, while the values at node 1, 2, 3 and 18, 19, 20 which are adjacent to the free surface are offset highly from 1.1. It is confirmed that the theoretical expression deduced could provide the results as accurately as the present numerical method. The reason of deviation for the two solutions is not only the influence of the free surface at both ends of the crack front, but also the stress coupling effect.

The case to be considered is the specimen subject to in-plane shear loading. In this case the generated $K_{II}$ values are rather constant in the inner part of the specimen, as shown in Figure 18, but for nodes near the free surface, they increase remarkably, which means that the singularity changed in the vicinity of the free surface. The stress intensity factor $K_{III}$ was generated, even if the external loading of the specimen does not include an out-of-plane shear component. This phenomenon is defined as the stress coupling effect [23,24]. The absolute values of $K_{III}$ were enlarged gradually from the middle node of crack tip to both sides, which illustrated the coupling effect enhanced by the free surface. Due to the absolute values of $K_{II}$ being less than $K_{III}$, especially in the inner part of the specimen, the coupling effect is called weak coupling.

Figure 19 shows the out-of-plane shear loading case. In this case, the pure III mode loading generated $K_{III}$ and $K_{II}$, and the absolute value $K_{II}$ of nodes adjacent to the free surfaces exceeds the value of $K_{III}$, which is called the strong coupling effect. Regarding the I-II-III mixed mode crack simulation, the induced $K_{II}$ and $K_{III}$ have a significant influence on their theoretical value due to the coupling effect, which can result in the deviation between $G\left(\theta \right)$ and $G$.

4. Conclusions

The energy release rate associated with stress intensity factors under mixed mode loading for aluminum alloy material has been investigated using both a numerical method and theoretical derivation. The present results demonstrate the simpler, accurate calculation process and accurate evaluation of the energy release rate. The following conclusions can be drawn:

(1)

A relatively simple and precise numerical method was established to evaluate the energy release rate associated with the stress intensity factors under mixed mode loading based on the concept that the energy release rate is equal to the change rate of the energy difference before and after crack kinking.

(2)

Based on the numerical method, 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, and the expression of the energy release rate with the crack propagation angle was obtained.

(3)

A theoretical expression of energy release rate at any propagation 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 is confirmed that the theoretical expression deduced could provide results as accurately as the present numerical method.

(4)

The present results are consistent with the experimental data. The error, which is lower than 5%, can be accepted considering that the specimens are not manufactured using ideal elastic material. Consequently, the present results demonstrate the simpler, accurate calculation process and the accurate evaluation of energy release rate.