Characteristic Tensor for Evaluation of Singular Stress Field under Mixed-mode Loadings

A characteristic tensor is defined using stress tensor averaged in a small circular domain at the crack tip and multiplied by the root of domain radius. It possesses the original stress tensor characteristics and has a simple relationship with conventional fracture-mechanics parameters. Therefore, it can be used to estimate stress intensity factors (SIFs) for cracks of arbitrary shape subjected to multiaxial stress loads. A characteristic tensor can also be used to estimate SIFs for kinked cracks. This study examines the relation between a characteristic tensor and SIFs to demonstrate the correlation between the characteristic tensor and fracture-mechanics parameters. Consequently, a single straight crack and a kinked crack of finite length existing in a twodimensional, infinite isotropic elastic body in a plane stress state, were considered to investigate the properties of the characteristic tensor under mixed-mode loadings. To demonstrate the practical utility of the characteristic tensor, the stress distribution obtained through finite element analysis (FEA) was used to estimate mixed-mode SIFs, and the values of estimated SIFs were compared with those obtained using an analytical solution. Results demonstrate that SIFs estimated under mixed-mode loadings exhibit a good agreement with the analytical values. This indicates that the proposed characteristictensor-based approach is effective in extracting features of singular stress fields at crack tips, and can be employed to estimate values of fracture-mechanics parameters, such as SIFs. Owing to its simplicity, the proposed approach can be easily incorporated in commercial FE codes for practical applications to simulate the crack-growth problem under both static and dynamic loading scenarios. The excellent applicability of the characteristic tensor greatly contributes to efficiency of the design process in industries.


Introduction
Fracture-mechanics parameters are a topic of intense research, as they regard the evaluation of singular stress fields near crack tips and generic crack-propagation analyses. The stress intensity factor (SIF) is a basic elastic fracture-mechanics parameter that is used to express the singular displacement and stress fields near crack tips. SIF values can be directly calculated from either displacement or stress fields obtained via finite element analysis (FEA) [Chan, Tuba and Wilson (1970)]. One of the most popular fracturemechanics parameters is Rice's J-integral [Rice (1968)], which was originally proposed as a contour integral around crack tips. Owing to its characteristics of being pathindependent during evaluation and applicability to both elastic and elastic-plastic mechanics, the J-integral has wide applications in practical crack-propagation analyses. However, in actual investigations dealing with stress fields obtained using the finite element method (FEM), the J-integral demonstrates an inherent inaccuracy [Moran and Shih (1987)], to avoid which Moran et al. developed the domain integral, which demonstrates a more robust and accurate computational performance in comparison with the J-integral when calculating contour integrals around crack tips. Additionally, under mixed-mode loadings, which define practical operating conditions for most structures, fracture-mechanics parameters contribute considerably to the assessment of material failure due to crack propagation. When considering mixed-mode loadings, the use of SIFs is preferred to that of the J-integral, because the J-integral does not provide an individual SIF for each mode separately. To extract mixed-mode SIFs from contour integrals, Chen et al. [Chen and Shield (1977)] developed interaction integrals to determine SIFs from J-integrals in cases involving actual deformations and auxiliary fields. Yau et al. [Yau, Wang and Corten (1980)] applied a J-integral-based interaction to a single straight crack formed in a 2D specimen and determined mixedmode SIFs under several loading conditions. Combining the auxiliary field and domain integral methods, Shih et al. [Shih and Asaro (1988)] developed domain forms of interaction energy integrals for analysis of 3D cracks. The method was employed for both 2D and 3D cracks. Gosz et al. [Gosz, Dolbow and Moran (1998)] employed domain forms of interaction-energy integrals to calculate mixed-mode SIFs of curved 3D cracks by considering additional terms in the domain integral. These additional terms vanish at a straight crack front. Gosz et al. [Gosz and Moran (2002)] evaluated several crackgeometry types using the method, thereby introducing an orthogonal curvilinear coordinate system that eliminates issues introduced by the consideration of auxiliary stress and displacement fields for crack shapes. Walters et al. [Walters, Paulino and Dodds (2005)] investigated the extent to which terms concerning the crack-front curvature contribute to the computation accuracy, and in turn, proposed a simplified formulation of the interaction integral with judicious levels of mesh refinement. Subsequently, interaction integrals based on domain forms, which have wide applications in industries, were implemented in non-linear FEM software [Gadallah, Osawa, Tanaka et al. (2018)]. A virtual crack closure-integral method (VCCM) offers another means to identify fracture-mechanics parameters in elastic bodies subjected to various loading conditions [Rybicki and Kanninen (1977)]. Through implementation of this method, the energy release rate can be derived from nodal forces and displacements calculated by FEA. The said method can be applied to both 2D and 3D bodies [Shivakumar, Tan and Newman (1988)]. SIFs can be estimated by considering the asymptotic stress field closer to the crack tip. In practical crack-propagation analyses, the method has been extended to be applied using a tetrahedral solid-element formulation suitable to 3D remeshing procedures for new surfaces created owing to crack extension ; Kawai, Okada and Araki (2008)]. Because singular-field evaluations (based on identification of fracture-mechanics parameters described above) are essential for efficient product development and manufacturing, extensive research has been conducted in this regard. A characteristic tensor can be derived from the characteristic stress field near crack tips, and is strongly related to fracture-mechanics parameters in principle as mentioned by Murakawa [Murakawa (2018)]. Each component of a characteristic tensor comprises an invariant stress singularity. Because these tensor components contain information regarding the intensity and direction of the said singularity, they can be considered a potential aid in the evaluation of crack initiation as well as its propagation direction under complex loading conditions encountered in practical analyses. A major advantage of characteristic tensor is its ease of for calculation. Although several parameters based on fracture mechanics and methods have been successfully implemented for industrial problems, their applications are limited because they are not as easy as the characteristic tensor. This is because conventional methods, such as domain-based interaction integrals and VCCM, require intensive and complex computation implementation. These limitations are critical, especially for dynamic analyses using explicit FEM. Because the time step, which is determined by element size and stress wave speed, is extremely small in explicit FEM, a large number of calculation steps is generally necessary. Hence, currently crack propagation computation based on reliable fracture-mechanics parameters is seldom performed using dynamic explicit FEM. On the contrary, computing procedures for the characteristic tensor are comparatively simple and easy to incorporate with commercial FEM software. Because the characteristic tensor is only derived from a stress field, it does not require special element formulations or complex computation procedures. If the characteristic tensor is estimated using an appropriate fitting function, the deterioration in prediction accuracy owing to a coarse mesh is not significant in comparison with that computed by the direct method [Chan, Tuba and Wilson (1970)]. Moreover, the direct relationship between the characteristic tensor and fracture-mechanics parameter makes it convenient to utilize conventional methods for crack-propagation analysis. Existing reliable methods based on fracture-mechanics parameters that were validated in extant studies can be used with the characteristic tensor. When used in combination with other element technologies, such as the extended finite element method (XFEM) [Moës, Dolbow and Belytschko (1999)], the tensor proves to be efficient in reducing the computing time. Analysis using the characteristic tensor can be effectively incorporated in industrial crack-propagation-analysis workflows, thereby enhancing the efficiency of the overall design cycle. To validate the effectiveness and accuracy of the characteristic tensor, it is essential to identify its relevance with regard to conventional fracture-mechanics parameters. This study aims to demonstrate the reliability associated with the use of a characteristic tensor by identifying its relationship with SIFs of structures subjected to mixed-mode loadings. In this paper, the definition of a characteristic tensor is first presented, followed by the verification of its properties in cases involving a single-straight crack subjected to uniaxial and mixed-mode loadings. Subsequent sections describe the several FEA-based numerical studies performed to investigate the feasibility of the characteristic tensor for practical analyses. A single straight crack and a kinked crack with finite kink length were considered, and their behavior was investigated under several loading conditions. In all cases, characteristic tensor components were determined from the stress distribution obtained via FEA. Mixed-mode SIFs were estimated using corresponding characteristic tensor components, and their values were compared with those obtained analytically, as reported in extant literature to demonstrate the reliability of the proposed method. Results obtained in this study demonstrates the feasibility of the characteristic tensor in general crack-propagation analyses. In Section 2, the definition of the characteristic tensor is presented, and then, its relationship with fracture-mechanics parameters under a mode-I loading case is presented. In Section 3, properties of the characteristic tensor under a mixed-mode loading case are discussed. The relation between each tensor component and SIFs is then introduced. Finally, Section 4 presents the results of several numerical simulations to evaluate the applicability of characteristic tensor.

Characteristic tensor
Considering a two-dimensional isotropic elastic body having a single straight crack of length 2 shown in Fig. 1, the mean stress near the crack tip can be defined as the average stress over the domain Ω, which corresponds to the zone enclosed within a circle of infinitesimally small radius . That is, Here, Ω denotes the volume of domain Ω, denotes the distance from the crack tip along the radial direction, and denotes the stress component at the crack tip in the Cartesian coordinate system. The characteristic tensor can be calculated by multiplying the square root of radius with the mean stress as The William's general stress function around the crack tip can be described as a superposition of polynomial functions [Williams (1956)] given by Here, is a function of , which represents the angle along the circumferential direction. Using the stress function relation, the corresponding stress distribution can be expressed as where is a function that depends on and its derivatives. The stress distribution described in Eq. (4) includes higher-order terms, which can be ignored if the value of is infinitesimally small (i.e., ⁄ ≪ 1). In this case, the stress function has a singularity of 1 √ ⁄ in the local neighborhood of a crack tip. Therefore, if the value of the singularity under the applied load and crack geometry is calculated, the characteristic tensor described in Eq. (2) becomes a constant near the crack tip. That is, Since the characteristic tensor has no dependency on any variable other than the singularity of 1 √ ⁄ , it is highly suitable to capture the singularity produced by the existence of cracks. Moreover, the characteristic tensor may exclusively be sufficient to assess the singularity of each stress component owing to its tensorial expression. As a primitive investigation of the characteristic tensor, a mode-I loading case, wherein a single straight crack of length 2 forms in a two-dimensional, isotropic elastic body subjected to a uniaxial load acting perpendicular to the crack as shown in Fig. 2, was considered in this study. The local coordinate system was placed at the crack tip, and the same was embedded along the crack direction. The observed stress distribution under the mode-I loading around a crack tip can be expressed as where denotes the stress component, and denotes the SIF under mode-I loading. Additionally, denotes the distance from the crack tip, and represents a function of angle , which determines the angle along the circumferential direction. Assuming ( , ) = �(1,1), (2,2)� and 0 ≤ ≤ , values of functions 11 and 22 can be determined using The second term on the right in Eq. (6) can be eliminated if → 0 . Thus, the characteristic tensor can be expressed as Eqs. (9) and (10) indicate that the characteristic tensor has a proportional relationship with mode-I SIF near the crack tip. However, the stress field around the crack contains a singular term 1 √ ⁄ as well as higher-order terms contained within the second term on the right in Eq. (6). This may be attributed to the averaging radius , which is generally defined to be of finite length. Therefore, the characteristic tensor derived from the stress field can be approximated using the polynomial function where , denotes the polynomial coefficient. Use of the Westergaard function [Westergaard (1939)] also yields the stress field around a single crack formed in a structure subjected to mode-I loading, and the corresponding stress distribution can be obtained analytically by solving the Westergaard function [Sun and Jin (2012)]. Therefore, the mean stress defined in Eq. (1) can be calculated by numerically integrating the analytical stress distribution. Lastly, the characteristic tensor defined in Eqs. (9) and (10) can be calculated using the mean stress and corresponding radius of the averaging region. Using the stress field along with Eqs. (9) and (10), the SIF value can be estimated using the below expressions.
As the averaging radius becomes infinitesimally small, the estimated SIF value approaches the theoretical value. On the other hand, the larger the radius, the greater is the observed deviation between the two values. This difference between theoretical and estimated SIF values stems from the second term-on the right of the equality in Eq. (11)which depends on the averaging radius. Therefore, it is necessary to extract the characteristic tensor comprising only the singular term 1 √ ⁄ when considering a practical integration domain with a finite averaging radius.   To be able to utilize a characteristic tensor effectively during practical analyses, it is important to understand its properties in more general cases. This section describes how the characteristic tensor is related to SIFs when considering structures subjected to mixed-mode loadings. Fig. 4 depicts a two-dimensional infinite isotropic elastic body with a single crack with a Cartesian coordinate system considered aligned with crack direction. The general stress field in the local neighborhood of the single crack can be described as the superposition of two fundamental modes-I and II. The resulting stress distribution can be expressed as where and denote SIFs corresponding to modes I and II, respectively; denotes a constant stress along the 1-direction in the local coordinate system; and and are functions of the variable , which determines the angle in the circumferential direction asymmetric to the crack surface. Values of the said functions in the interval 0 ≤ ≤ can be obtained using the following relations.
Considering → 0, all higher-order terms in Eq. (11) can be ignored. Using Eq. (1), the mean stress can be evaluated using the expression Considering the asymmetric property of the stress distribution around a crack tip, several terms concerning the integration of 12 , 11 , and 22 conveniently vanish. Thus, the characteristic tensor can be evaluated using the expressions 11 = √ 11 = 11 + √ (23) χ 12 = √ 12 = 12 (25) where coefficients 11 , 22 , and 12 are constants derived by integrating with respect to ( 11 =1.6, 22 =2.4, and 12 =1.6). Likewise, coefficient = �8 (9 3 ) ⁄ is also a constant evaluated by integrating the stress-distribution function with respect to radius . These equations demonstrate that each component of the characteristic tensor, which is based on the coordinate system aligned with the crack direction, is directly proportional to SIFs under mixed-mode loading conditions. This implies that the characteristic tensor is directly related to reliable fracture-mechanics parameters in principle, and the same can be used to identify conventional parameters that can be used to tackle the problem of crack initiation and propagation. However, as described in the previous section, the characteristic tensor estimated using a finite averaging radius is only an approximation obtained from a polynomial function, and its value cannot be considered constant owing to the influence of terms that depend on the averaging radius. Therefore, some numerical manipulations are required to extract constant terms from the said polynomial function and identify valid parameters for crack analyses. In this study, the least-square method was employed, as described in the following section.

Numerical study
To verify applicability of the characteristic tensor under mixed-mode loading, several numerical simulations were performed. SIF values were estimated using the characteristic tensor calculated from the stress distribution obtained for mixed-mode loading scenarios, and the mixed-mode stress state was realized in a two-dimensional infinite plane body having a single crack inclined to the uniaxial loading direction. At first, a single straight crack was considered, and mixed-mode SIFs are estimated for different inclination angles followed by their comparison against theoretical values. Next, a single kinked crack having a finite kink length was considered, and corresponding SIF values were determined. All characteristic tensors were evaluated using the stress distribution obtained via FEA.

Numerical model
For numerical analysis performed in this study, a two-dimensional infinite thin plate containing a single crack of length 2 was considered. Fig. 5 depicts a schematic of the numerical model. An in-plane biaxial load was applied to the model as a boundary condition and the crack was inclined at angle with respect to the loading direction to ensure attainment of several mixed-mode stress states around the crack tip. All loads were applied parallel to either the 1-or 2-direction of the global coordinate system. In addition to the global coordinate system, a local coordinate system aligned with the crack direction was defined in the model to facilitate calculation of the characteristic tensor. The finite element model used in this analysis is depicted in Fig. 6. The size of the model was determined to reproduce the infinite plane considered in this numerical study. A 4-noded fully integrated 2D plate element, which is under plane stress condition, was employed for geometry modeling. The material was assumed be isotropic and elastic.

Calculation of characteristic tensor from computed stress field
Numerical integration of a stress field is required for calculating the characteristic tensor. To calculate the mean stress, defined in Eq. (1), stress components at integration points must be integrated over the averaging domain using Gaussian quadrature. The resulting expression for mean stress calculation can be expressed as where denotes the averaging radius, denotes the total number of integration points in the averaging domain, � denotes the stress component, denotes the weight of the Gaussian quadrature formula, and det denotes the Jacobian determinant. In Eq. (26), all quantities with index possess a value at the corresponding integration point. From Eq.
(2), the characteristic tensor can be obtained as As described in Sections 3 and 4, the characteristic tensor derived via practical stress analysis contains additional terms that depend on the averaging radius. Thus, for accurate evaluation of singularities at cracks, extraction of valid terms is necessary. The characteristic tensor derived from Eq. (27) can be approximated as which contains terms up to = 2 in Eq. (11). In this study, component was extracted from Eq. (28) using the least-square method. Eventually, SIF values were estimated using Eqs. (24) and (25)

Estimation of mixed-mode SIFs for single straight crack
In this study, mixed-mode SIFs for a single straight crack were estimated using the characteristic tensor calculated from the stress field obtained via FEA. To obtain several mixed-mode stress states around the crack tip, biaxial stresses were applied along boundaries of the cracked body, thereby changing the stress ratio-defined as = 1 2 ⁄to values of 0, 0.5, and 1.0 (Fig. 5). Additionally, in each case, the crack angle was changed in the 0-75° range in 15° increments. Three cases of element size ( = 0.1 , 0.01 , and 0.001 ) were considered near the crack tip. Fig. 7 depicts the mesh pattern for the case = 0.1 with = 30°, whereas Fig. 8 depicts the corresponding mesh for each case near the crack tip. Values of estimated mixed-mode SIFs were compared to those theoretically obtained in extant studies [Anderson (2005)]. Analytical expressions used for SIF calculation were given by = 0 (cos 2 + sin 2 ) (29) and = 0 (sin cos )(1 − ) (30) where and denote mixed-mode SIFs, and 0 denotes the SIF for mode I when = 0 and = 0. Figs. 9, 10, and 11 compare mixed-mode SIF values obtained for different inclination angles. All SIF values were normalized with respect to the theoretical SIF obtained for = 0, which remains constant irrespective of the applied stress ratio. The graph in Fig. 9 refers to the case with stress ratio = 0, which implies uniaxial loading along the structures boundary. The graph in Fig. 10 corresponds to the case with = 0.5, whereas that in Fig. 11 refers to the = 1.0 case, thereby implying application of pure biaxial tension along the structure boundary. Because its value equals zero under application of pure biaxial tension, the SIF obtained for mode II in the = 1 case is not depicted in Fig.  11. As observed, estimated SIF values demonstrated good agreement with corresponding theoretical values obtained for different stress ratios and inclined angles. Tabs. 2, 3, and 4 list percentage errors incurred during estimation of mixed-mode SIFs in this study. In cases with = 0.001 and 0.01 , the error incurred was observed to be less than 1%. In the case with = 0.1 , which provides acceptable computational costs even in the explicit FEM, in which the element length directly effects the incremental time step size, the error incurred was less than 3%. These results imply that derivation of the characteristic tensor is helpful in evaluating singularities near crack tips along with conventional fracture-mechanics parameters.    Lastly, in this study, estimated SIF values were compared against those obtained using conventional methods-direct method [Chan, Tuba and Wilson (1970)] and domain-integral method [Lei, O'Dowd and Webster (2000)]. In the direct method, SIF can be estimated using either the displacement or stress around the crack tip, where the displacement perpendicular to the crack direction under mode-I loading can be expressed as In the above expression, G denotes shear modulus, and denotes the Poisson ratio. The corresponding SIF was estimated using nodal displacements on the crack surface as follows: Tab. 5 summarizes the percentage error incurred during SIF estimation at = 0 and = 0 using the three methods. In the case of the direct method, the linear approximation was employed for fitting nodal displacements. As can be observed, for each element size, the domain-integral method predicts SIFs with errors measuring less than 1%. The direct method estimates SIF with good accuracy in the case with = 0.001 . However, the error incurred increases to 10% in the = 0.1 case. The proposed characteristictensor-based approach estimates SIF accurately when using a fine mesh. Furthermore, the deterioration in estimation accuracy when using a coarse mesh is much smaller compared to that observed when using the direct method. Through use of an appropriate fitting function, the least-square method works effectively to extract characteristic tensors containing the 1 √ ⁄ singularity. One aspect of characteristic tensor is its simplicity of derivation. The method of deriving SIFs via characteristic tensors is easier than those of other methods. In the direct method, though the derivation of the SIFs is similar to that of the aforementioned tensors and the least square method can be also applied for the fittings, nodal displacements must be identified on the same straight line toward the crack tip. In the domain integral method, the mesh orientation must be designed and searched during computation such that the contour of elements forms a circle that starts from the bottom surface and ends at the upper surface of the crack. On the contrary, the calculation scheme of the proposed method involves adding the stress components of integration points in the corresponding radius. Therefore, the impact to computation costs is not too large. Evaluation with the characteristic tensor can be an effective solution for assessing singularity, and balancing cost and accuracy.

Estimation of mixed-mode SIFs in single kinked crack having finite kink length
In general, the crack direction changes during crack growth. Consequently, a curved crack path forms as the crack deviates from its original direction. Therefore, SIF estimation for kinked cracks is important from the viewpoint of evaluating singular fields for comprehensive crack-propagation analyses. In this study, SIFs for kinked cracks with finite kink length were estimated in a manner similar to that described in the previous section. Figure 12: Single kinked crack with finite kink length subjected to uniaxial tension Fig. 12 depicts a schematic of the numerical model. As can be seen the crack is aligned with the 1-direction of the global coordinate system. One of the tips of the single crack is inclined at angle α with respect to the original crack direction, and its kink length has been considered finite. In this study, the kink angle was increased in 15° increments in the range of 15-60°. For each kink angle, two cases of kink length-0.1 and 0.2 -were investigated. Dimensions of the FE model were the same as those depicted in Fig. 6. The element size of the model was set as 0.001 around the kinked-crack tip. The boundary stress ratio was considered to be = 0, and only the stress along the 2-direction of the global system was applied to the elastic structure's boundary. The characteristic tensor was calculated with respect to the local coordinate system aligned with the kink direction.
Corresponding SIF values were estimated as described in previous sections. These SIF values were compared against analytical results reported in extant literature [Kitagawa, Yuuki and Ohira (1975)] with SIF values being normalized with respect to that obtained for a straight crack of length c. The said normalized SIFs (F 1 and F2) could be expressed as [Kitagawa, Yuuki and Ohira (1975) The length c of the straight crack is given by = 2 + cos Tab. 6 lists reference values of normalized SIF along with the percentage error incurred during their estimation as well as the characteristic tensor. For all kink angles considered, the observed error measured less than 1% for all lengths. These results demonstrate that SIF estimation via use of the characteristic tensor is accurate, and that the characteristic tensor is reliable for use even in cases involving kinked cracks.

Conclusions
A characteristic tensor is a second-order tensor-similar to a stress tensor-and characterizes the singular stress field near crack tips in linear elastic bodies. A characteristic tensor bears a direct relationship with the conventional fracture-mechanics parameters, such as SIF and energy release rate. Several numerical computations were performed to investigate the theoretical validity of the characteristic tensor as well as the expected accuracy associated with its application in the evaluation of fracture-mechanics parameters subjected to multiaxial stress states. Major conclusions drawn from this study are  The characteristic tensor demonstrates a direct relationship with fracture-mechanics parameters under both uniaxial and mixed-mode loadings.  The characteristic tensor can be easily calculated from stress distributions of FEA by employing the least-square method.  Estimated SIF values for a single crack under mixed-mode loading demonstrate good agreement with their analytical counterparts. Additionally, for a single kinked crack, SIFs can be accurately estimated.  The error incurred during SIF estimation is acceptable, albeit the mesh size used for stress analyses is substantially large.  These results suggest that the proposed characteristic tensor characterizes and quantifies the singularity at the crack tip and can be applied to analyze practical crack propagation problems.

Conflicts of Interest:
The authors declare that they have no conflicts of interest to report regarding the present study.