Averaged strain energy density estimated rapidly from the nodal stresses by FEM for cracks under mixed mode loadings including the T-stress contribution

The present contribution reviews a recently proposed method to rapidly estimate the averaged SED at the tip of short as well as long cracks under in-plane I+II and long cracks under out-of-plane I+III mixed mode loadings. Short cracks are distinguished from long cracks by considering that the stress fields within the control volume of short cracks are no longer governed solely by the stress intensity factors (SIFs), but also the contribution of higher order terms, and primarily the T-stress, becomes significant to estimate the averaged SED. According to the proposed method, the averaged SED is calculated using the linear elastic nodal stresses evaluated by FEM either at the crack tip, to account for the SIFs contribution, and at selected FE nodes of the crack free edges, to include the T-stress contribution. The advantage of the proposed approach is two-fold: coarse FE meshes can be adopted; moreover, geometrical modelling the control volume is no longer necessary. To validate the approach, cracked plates subjected to in-plane I+II mixed mode loading as well as bars weakened by circumferential outer cracks subjected to out-of-plane mixed mode I+III loading have been analysed. A comparison between approximate values of the averaged SED according to the nodal stress approach and those derived directly from the FE strain energy adopting very refined FE meshes has been successfully performed.


INTRODUCTION
ealing with cracks under in-plane mixed mode I+II loading conditions, according to Williams [1], the local stress fields expressed in terms of Cartesian stress components as functions of the polar coordinates (r,θ), with origin at the crack tip (Fig. 1a), can be written in the following form: With reference to cracks under mode III loading conditions, the asymptotic, singular stress distributions have been determined by Qian and Hasebe [2], following Williams' procedure [1].The local stress field in terms of Cartesian stress components as functions of the polar coordinates (r,θ), with origin at the crack tip (Fig. 1b), is the following:   The mode I and mode II Stress Intensity Factors (SIFs) can be defined according to Gross and Mendelson [3] by means of Eqns.( 3) and (4), respectively.
In the context of fracture mechanics it is largely assumed that the stress field in the close neighborhood of the crack tip can be properly characterized by means of the coefficients of the leading order terms, i.e. the SIFs.However, detailed analyses reported in the literature have highlighted the fundamental role of the T-stress in defining the stress state close to the crack tip [4][5][6][7][8].Larsson and Carlson [4] and later Rice [5] argued on the effect of T-stress on the plastic zone ahead of the crack tip in materials characterized by elastic-plastic behaviour.The influence of T-stress on failure mechanisms of brittle materials was investigated by Ayatollahi et al. [6,8] and Fett and Munz [7], who employed a modified maximum tangential stress approach (MTS), taking into account mode I and mode II SIFs, T-stress and a material-dependent length parameter.
The combined effects of SIFs and T-stress on structural strength problems of cracked components under mixed mode I+II+III loadings can be easily evaluated by means of the strain energy density (SED) averaged over a control volume, thought of as a material property and modelled as a circular sector of radius R 0 , as shown in Fig. 2, according to Lazzarin and Zambardi [9].The averaged SED criterion has been widely adopted in the recent literature for static [10] and fatigue [9][10][11] strength assessments.Dealing with a general mixed mode crack problem, the averaged SED can be expressed in closed-form as a function of the SIFs, i.e.KI, KII and KIII, and of the T-stress, i.e.T, according to the following analytical expression [11,12]: In the above equation, e 1 , e 2 and e 3 are three known parameters dependent on the opening angle of a general sharp V-notch and on the Poisson's ratio ν [9], while E is the Young's modulus of the considered material.The use of Eqn.(7) in engineering problems presents a major drawback, since very refined meshes are required to evaluate the SIFs and the T-stress on the basis of definitions (3)-( 5) and ( 6), respectively, applied to a number of stress-distance FE data.This is due to the fact that the entire local stress field must be determined accurately.The practical application is even more time-consuming in the case of 3D FE models.However, the averaged SED can be evaluated directly from the FE results, FEM W , by summing the strain-energies W FEM,i calculated for each i-th finite element belonging to the control volume and by dividing by the volume (or area in 2D problems, A in Fig. 2): Equation ( 8) represents the so-called direct approach to calculate the averaged SED.According to a recent contribution of Lazzarin et al. [13] coarse FE meshes within the control volume A can be used.
Recently, a method to rapidly estimate the averaged SED at the tip of short as well as long cracks under in-plane I+II [14,15] and long cracks under out-of-plane I+III [16] mixed mode loadings has been proposed.It is based on the nodal stresses evaluated from finite element (FE) analyses, according to the nodal stress approach: the averaged SED is calculated using the linear elastic nodal stresses evaluated by FEM either at the crack tip, to account for the SIFs contribution according to the peak stress method (PSM), and at selected FE nodes of the crack free edges, to include the T-stress contribution.The advantage of the proposed approach is two-fold: there is no need of mesh refinements in the close neighbourhood of the points of singularity, so that coarse FE meshes can be adopted; moreover, geometrical modelling the control volume in FE models is no longer necessary.The present contribution reviews the nodal stress approach and its validation, which is based on the FE analyses of cracked plates subjected to in-plane I+II mixed mode loading as well as bars weakened by circumferential outer cracks subjected to out-of-plane mixed mode I+III loading, while varying (i) the crack lengths, (ii) the mode mixity and (iii) the finite element size adopted in the numerical analyses.Table 1: Geometrical and loading parameters taken into consideration in Refs.[14][15][16].All combinations have been analysed, provided that the ratio a/d was greater than the minimum feasible one, i.e. a/d = 1.

RANGE OF APPLICABILITY OF THE SED EXPRESSION (7)
In-plane I+II mixed mode loading onsidering the in-plane I+II mixed mode crack problem of Fig. 1a, exact values of the averaged SED, FEM W (Eqn. ( 8)), have been evaluated for the geometrical and loading cases reported in Tab. 1, by adopting the direct approach with very refined meshes (patterns with about 500 FE within the reference volume).The 'exact' K I and K II SIFs and T-stress have been evaluated for the same crack problems by using definitions (3), ( 4) and ( 6), respectively, applied to FE results of numerical analyses characterized by very refined meshes, where the element size close to the crack tip was reduced to approximately 10 -5 mm.After that, also the analytical SED, AN W has been calculated from Eqn. (7) and deviated from the exact value, FEM W , by less than 5% for all considered cases.It should be noted that when the T-stress contribution is negligible, only K I and K II contribute to the averaged SED, K III being null, and Eqn.(7) simplifies to: C Fig. 3a plots the ratio AN,I+II FEM W /W for three mode mixity ratios MM, i.e. 0 (pure mode I), 0.5 and 0.8, where MM = K II /(K I + K II ).The figure highlights that for a crack size equal to the radius of the material dependent control volume (a/R 0 = 1) and pure mode I loading (MM = 0) the error in the averaged SED calculation is about 42% if the T-stress contribution is neglected, i.e. if Eqn. ( 9) is used in place of Eqn.(7).However, for the same crack size to control radius ratio, the error is decreased to 17% if MM = 0.5 and further reduced to only 2% if MM = 0.8.Fig. 3a shows that when reduced crack size to control radius ratios (a/R0 < 10) are considered, Eqn.(7), which includes the T-stress, must be adopted.Accordingly, the condition a/R 0 < 10 identifies the short crack case, while the condition a/R 0  10 identifies the long crack case.

Out-of-plane I+III mixed mode loading
Dealing with the out-of-plane I+III mixed mode crack problem of Fig. 1b, exact values of the averaged SED, FEM W (Eqn. ( 8)), have been evaluated for the geometrical and loading cases reported in Tab. 1, by adopting the direct approach with very refined meshes.It is worth noting that when the T-stress contribution is negligible, only KI and KIII contribute to the averaged SED, KII being null, and Eqn.(7) simplifies to: Once evaluated the 'exact' K I and K III SIFs by using definitions (3) and ( 5), respectively, applied to FE results of numerical analyses characterized by very refined meshes, also the analytical SED, AN,I+III W can be calculated from (Eqn.( 10)).Fig. 3b plots the ratio AN,I+III FEM W /W for pure mode I (MM = 0) and pure mode III (MM = 1) loading, where MM = KIII/(KI + KIII).The figure highlights that for a crack size equal to the radius of the material dependent control volume (a/R0 = 1) the error in the averaged SED estimation is about -15% for pure mode I loading (MM = 0), while it increases to about +30% for pure mode III loading (MM = 1).These deviations are due to the contribution not only of the T-stress, but also of further higher order non-singular terms, O(r 1/2 ) in Eqns.( 1) and ( 2), which are needed to account for the free-edge boundary conditions as discussed in detail in [16], but they are disregarded in the analytical expression, Eqn.(7).It should be noted that averaged SED expressions which account for the contribution of further higher-order, non-singular terms, O(r 1/2 ) in Eqns.( 1) and (2), are not currently available in the literature.That said, it can be observed from Fig. 3b that the analytical SED, W , (Eqn.( 10)), deviates from the exact value, W (Eqn. ( 8)), by less than 5% for a ratio a/R 0 ≥ 10, i.e. the long crack case.Under this condition, the contribution of higher-order, non-singular terms are believed to be negligible from an engineering point of view, so that Eqn. ( 7) is applicable.Therefore, cracks under out-of-plane I+III mixed mode loading characterised by a/R 0 > 10 will be analysed in the following, in order to comply with the range of applicability of Eqn.(7), according to results shown in Fig. 3b.

THE PEAK STRESS METHOD TO RAPIDLY EVALUATE K I , K II AND K III
he Peak Stress Method (PSM) is an approximate numerical technique to evaluate the SIFs.The PSM takes its origins by a numerical technique proposed by Nisitani and Teranishi [17] to rapidly estimate by FEM the mode I SIF of a crack emanating from an ellipsoidal cavity.A theoretical justification to the PSM has been provided later on and the method has been extended also to sharp and open V-notches in order to rapidly evaluate the mode I Notch Stress Intensity Factor (NSIF) [18].Subsequently, the PSM has been formalised to include also cracked components under mode II loading conditions [19] and open V-notches subjected to pure mode III (anti-plane) stresses [20].In more detail, the PSM enables to rapidly estimate the SIFs K I , K II and K III (Eqns.( 3)-( 5)) from the crack tip singular, linear elastic, opening, sliding and tearing FE peak stresses σyy,peak, τxy,peak and τyz,peak, respectively, which are referred to the bisector line according to Fig. 4a, concerning the in-plane stress components, and Fig. 4c, as to the out-of-plane stress component.More precisely, the following expressions are valid [18][19][20]: * 0.5 , 1.38 ** 0.5 ,

1.93
where d is the so-called 'global element size' parameter to input in Ansys ® FE code, i.e. the average FE size adopted by the free mesh generation algorithm available in the FE code.Eqns.( 11)-( 13) were derived using particular 2D or 3D finite element types and sizes, so that a range of applicability exists, which has been presented in detail in previous contributions [18][19][20], to which the reader is Here it is worth recalling that for mode I loading (Eqn.( 11)) the mesh density ratio a/d that can be adopted in numerical analyses must be a/d  3, a being the minimum between the crack and the ligament lengths; for mode II loading (Eqn.( 12)) more refined meshes are required, the mesh density ratio a/d having to satisfy a/d  14; in case of mode III loading (Eqn.( 13)) the condition a/d  3 must again be satisfied.
As an example, Fig. 4b shows a free mesh where d = 0.15 mm was input in Ansys ® software, while Fig. 4d shows a free mesh where the average FE size d is in constant proportion with the crack length a, i.e. a/d = 3.The mesh patterns shown in Figs.4b,d are as coarse as possible to estimate the averaged SED with a 10% error using next Eqn.(17).It is important to underline that no additional input parameters other than d and no additional special settings are required to generate an FE mesh according to the PSM.When Eqns. ( 11)- (13) were calibrated [18][19][20], the 'exact' K I , K II and K III SIFs were evaluated using definitions (3)-( 5), respectively, applied to FE results of numerical analyses characterized by very refined meshes, where the element size close to the crack tip was reduced to approximately 10 -5 mm.Therefore, the FE size required to estimate K I , K II and K III from σ yy,peak , τ xy,peak and τ yz,peak , respectively, is likely to be some orders of magnitude larger than that needed to directly calculate the local stress fields in order to apply definitions (3)-( 5).Moreover, while Eqns.( 3)-( 5) require to process a number of stress-distance numerical results, the PSM requires a single stress value to evaluate the SIFs.

A FE-BASED TECHNIQUE TO EVALUATE RAPIDLY THE T-STRESS
he analytical expressions of the stress components σ xx along the crack free edges are obtained substituting the polar coordinate θ = +π and -π in Eqn. ( 1) and are given in Eqns.(14a) and (14b), respectively [12]: Therefore, the T-stress contribution can be derived according to Eqn. (15), as previously highlighted by Ayatollahi et al. [6], Lazzarin et al. [11] and Radaj [21]: T Equation ( 15) suggests to evaluate the T-stress numerically using the nodal stresses σxx along the crack free edges.The obtained results are shown in Fig. 5a for a mesh density ratio a/d = 2 applied to the crack problem of previous Fig.1a.Fig. 5a show that due to numerical errors caused by the crack tip singularity, Eqn. ( 15) based on FE results is satisfied with a reduced error lower than 3% only at a distance from the crack tip r ≥ 2d.
On the basis of the obtained results, a FE-based technique to rapidly evaluate the T-stress can be defined according to the following expression: To verify the applicability of Eqn. ( 16) to the considered small cracks subjected to mixed mode I+II loading (see Tab. 1), the ratio between the Nodal T-stress according to Eqn. ( 16) and the exact T-stress is shown in Fig. 5b for a mode mixity ratio MM = 0.50, the results for other MM values being identical.Fig. 5b shows that in all cases the minimum feasible mesh density ratio a/d = 2 assures the applicability of Eqn. ( 16), because all numerical results fall within a restricted scatter-band of ±3%.

THE NODAL STRESS APPROACH TO RAPIDLY ESTIMATE THE AVERAGED SED WITH INCLUSION OF THE T-STRESS
n this Section, the nodal stress (NS) approach to estimate the averaged SED for short as well long cracks under inplane I+II [14,15] and long cracks under out-of-plane I+III [16] mixed mode loading is recalled.The technique is referred to as the nodal stress approach that can immediately be formalized by substituting Eqns.( 11)-( 13) and ( 16) into the appropriate analytical SED formulation, Eqn.(7):  Equation (17) shows that only few selected nodal stresses calculated from coarse FE mesh patterns can be used to estimate the averaged SED.Fig. 6 shows a sketch of the five nodal stresses involved in Eqn. ( 17): the crack tip, linear elastic opening (σyy,peak), sliding (τxy,peak) and tearing (τyz,peak) peak stresses referred to the crack bisector line and the linear elastic stresses σxx evaluated at the FE nodes located along the crack free edges at r = 2d.Coarse FE mesh patterns with an average FE size equal to d are adopted.

VALIDATION OF THE NODAL STRESS APPROACH TO ESTIMATE THE AVERAGED SED
Short as well long cracks under in-plane I+II mixed mode loading [14,15] o validate the nodal stress approach based on Eqn.(17), an infinite plate weakened by a central crack and subjected to in-plane mixed mode loading was considered according to Fig. 1a.Several crack lengths 2a (from 1 to 160 mm) have been considered, while width W and length L of the plate were set both equal to 10 times the crack length.The mean FE size d to evaluate σ yy,peak , τ xy,peak , σ xx,θ=π,r=2d and σ xx,θ=-π,r=2d (see Fig. 6) in Eqn.(17) was varied from 0.0125 to 10 mm.All different geometrical and loading parameters taken into account are listed in Tab. 1.According to the PSM, FE analyses have been carried out by using Ansys ® software and by adopting free mesh patterns consisting of 4-node quadrilateral elements (PLANE 182 with K-option 1 set to 3).Exact values of the averaged SED, W (Eqn. ( 8)), have been evaluated by adopting the direct approach with very refined meshes (patterns with about 500 FE within the reference volume).FIgs.7a,b,c report the ratio between approximate ( NS W , nodal stress approach Eqn. ( 17)) and exact ( FEM W , direct approach Eqn. ( 8)) averaged SED values for selected mode mixity ratios MM.The ratio W /W is seen to converge to unity inside a ±10% scatter-band for all considered MMs.In particular, Figs.7a,b,c show that convergence occurs for mesh density ratios a/d greater than 3 for MM = 0, 10 for MM = 0.50 and 14 for MM = 1.The obtained results show that the minimum mesh density ratio a/d to apply the nodal stress approach increases with increasing the mode mixity ratio MM.The minimum mesh density ratio to apply Eqn.(17) with a given level of approximation has been determined in [14,15] depending on the mode mixity ratio (MM) and it is not reported here for sake of brevity.[16] To validate the nodal stress approach based on Eqn.(17) under out-of-plane mixed mode loading, numerical analyses have been performed by considering a bar weakened by a circumferential outer crack according to Fig. 1b.Several crack lengths a (from 3 to 50 mm) have been considered, while diameter D and length L of the bar were set both equal to 10 times the crack length.The average FE size d to evaluate the peak stresses σyy,peak and τyz,peak (see Fig.  8)), have been evaluated by adopting the direct approach with very refined meshes.Fig. 7d reports the ratio between approximate (W , nodal stress approach Eqn. ( 17)) and exact ( FEM W , direct approach Eqn. ( 8)) averaged SED values for all mode mixity ratios MM.The ratio W /W is seen to converge to unity inside a ±10% scatter-band.In particular, Fig. 7d shows that convergence occurs for a mesh density ratio a/d greater than 3 for all mode mixity ratios MM taken into account.The obtained results show that the minimum crack size to FE size ratio a/d to apply the nodal stress approach (Eqn.( 17)) remains constant regardless of the mode mixity ratio MM.This differs from what was obtained previously dealing with cracks subjected to in-plane mixed mode I+II loading, for which the minimum mesh density ratio a/d to apply the nodal stress approach increased with increasing the mode mixity ratio MM, since mode II loading is more demanding in terms of mesh density ratio a/d than mode I loading.

CONCLUSIONS
he present contribution has reviewed the nodal stress approach recently proposed to estimate the averaged strain energy density (SED) of mixed-mode I+II and I+III crack tip fields including the T-stress contribution.The method takes five FE nodal stresses calculated with coarse FE meshes made of four-node, linear quadrilateral finite elements: three of them are the singular, linear elastic crack tip opening, sliding and tearing peak stresses, respectively, which take into account the stress intensity factor contribution; the remaining two ones are the nodal stresses evaluated along the crack free edges at a selected distance from the crack tip and take into account the T-stress contribution.The conclusions can be summarised as follows:  Taking advantage of the closed-form expression of the averaged SED (Eqn.7), the nodal stress approach has been formalised according to Eqn. (17);  For short as well as long cracks under in-plane I+II mixed mode loading, more refined FE mesh patterns are required, the higher the mode mixity ratio MM is.In particular, the minimum ratio a/d between the semi-crack length a and the average FE size d to apply the nodal stress approach with a level of approximation equal to 10% is found to be equal to 3 in the case of pure mode I (MM = 0), 10 in the case of mixed mode I+II with MM = 0.50 and 14 for pure mode II loading (MM = 1). For long cracks under out-of-plane I+III mixed mode loading, the minimum mesh density ratio a/d to apply the nodal stress approach with a level of approximation equal to 10% is equal to 3, independent of the mode mixity. Even though the averaged SED can be evaluated directly by means of FE analyses using coarse meshes inside the control volume, the T-stress contribution being automatically included, nonetheless some additional advantages of the nodal stress approach to estimate the averaged SED can be singled out: (i) only the linear elastic nodal stresses calculated at selected FE nodes are necessary; (ii) geometrical modelling the control volume in FE models is no longer necessary; (iii) the adopted FE meshes are coarse.

Figure 1 :
Figure 1: Cartesian stress components and polar coordinates with origin at the crack tip for (a) in-plane mixed mode I+II crack problem and (b) out-of-plane mixed mode I+III crack problem.

Figure 2 :
Figure 2: Strain energy density averaged over a control volume (area) of radius R 0 surrounding the crack tip, .W W= A . for (a) in-plane

Figure 5 :
Figure 5: In-plane mixed mode I+II crack problem: T-stress evaluated according to the nodal stress approach (Eqn.(16)).(a) FE stresses σ xx along the crack free edges by adopting a mesh density ratio a/d = 2. Considered case 2a = 3 mm and MM = 0.50.(b) Ratio between approximate and exact T-stress versus the mesh density ratio for the case MM = 0.50.

Figure 6 :
Fig. 1a 2α = 0°2 a = 3 mm MM = 0.50 a/d = 2 6, nodal stresses σxx,θ=±π,r=2d are not required, T-stress contribution being negligible for long cracks) in Eqn.(17) was varied from 0.05 to 10 mm.All different geometrical and loading parameters taken into account are listed in Tab. 1.According to the PSM, FE analyses have been carried out by means of Ansys ® software and by adopting free mesh patterns consisting of two-dimensional, harmonic, 4-elements (PLANE 25 of Ansys ® element library).Exact values of the averaged SED, W (Eqn. (

Figure 7 :
Figure 7: Ratio between approximate (W ) and exact (W ) averaged SED versus the mesh density ratio for (a), (b) and (c) the inplane mixed mode I+II crack problem of Fig. 1a and (d) the out-of-plane mixed mode I+III crack problem of Fig. 1b.W according to the nodal stress approach, Eqn.(17); W according to the direct approach, Eqn.(8).