Determination of stress intensity factors for elements with sharp corner located on the interface of a bi-material structure or homogeneous material

This paper presents a new universal analytical and numerical method allowing for determination of stress intensity factors (SIFs) for notches and cracks located in both a homogeneous and a heterogeneous material. An advantage of the proposed method is that it does not require knowledge of an analytical description of singular stress/displacement fields or connection of SIFs to energy parameters such as energy release factor. In this method, a universal analytical function has been used, which in combination with data obtained using finite element method (FEM) allows for a direct determination of stress intensity factors. One parameter that is necessary to know when using the proposed method is the eigenvalue λ\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\lambda $$\end{document}. The characteristic equation allowing for determination of the eigenvalues for any corner has been given herein. What is more, also a criterion clearly defining the nodes is determined, from which FEM numerical results are implemented to the developed analytical function. In order to verify the proposed method, for a selected group of geometrical and material structures with sharp corner, values of stress intensity factors were determined and compared to data available in the literature. Satisfactory compliance of obtained results with literature ones was found.


Introduction
With the determination of durability and strength of structural elements with sharp notches and cracks, some difficulty is involved. According to theoretical solutions based on elasticity theory, around the concentrator tip, values of both stresses and strains go ad infinitum. Use of the regular method involving comparison of these values to permissible stresses or strains is pointless, and thus a completely different approach is required. Therefore, for forecasting the durability and strength of materials and structural elements with material defects, methods based on fracture mechanics are used. Determination of critical service loads most frequently requires knowledge of stress fields around the concentrator tip area where a local stress concentration is formed. In the literature, analytical solutions can be found which are based on elasticity theory, where particular components of stresses (in a polar coordinate system located in the concentrator tip, Fig. 1) are described by the following equation-σ i j (r, ϕ) = N n=1 K In r λ n −1 f I i j (ϕ) + N n=1 K IIn r λ n −1 f II i j (ϕ) where f I i j (ϕ), f II i j (ϕ) are known eigenfunctions including boundary conditions of the issue, λ n -eigenvalues, n-number of terms of the asymptotic solution, K In , K IIn -stress intensity factors (SIFs).
The description of mechanical fields most frequently requires the use of the first power series term (n = 1, λ n = λ 1 = λ, K In = K I1 = K I , K IIn = K II1 = K II ), including singular term-eigenvalue λ < 1. Concrete analytical solutions and methods of their obtaining for the sharp corners located in a homogeneous material can be found, e.g. in papers [1][2][3].
The author of paper [3], using the eigenfunction expansion method, for the notch located in a homogeneous material obtained the so-called characteristic equation. The solution of this equation shows that the eigenvalue λ is not constant and depends on the notch tip angle. What is more, in the range of < 0, 1 > it occurs only once (single singular) and always belongs to a set of real numbers. In case if a sharp corner occurs in an inhomogeneous material, eigenvalues can be complex numbers-λ = λ r + iδ. For the very first time, such fact was reported for a case of an interfacial crack in paper [4]. With the presence of an imaginary part δ of the exponent λ characteristic oscillating terms of type sin[δln(r )] and cos[δln(r )] are introduced to the solution of stress and displacement fields, which increase the amplitude and "frequency" of changes in stresses and displacements, as the radius r goes to zero. Further investigations given, for example, in papers [5,6], have shown that the real part of the eigenvalue is always constant and equals λ r = 1/2, and the imaginary part depends on material constants. An other stress concentrator, frequently occurring in heterogeneous materials, is structural notch (Fig. 2). Such issue was dealt with by, for example, the authors of papers [7][8][9]12]. On the basis of their investigations, it can be claimed that in the case of a notch located on the interface of a bimaterial structure, stress and displacement fields can be described by both the real and the complex eigenvalue λ. The value of the exponent λ depends on the notch geometry and material constants; however, also multiple singularities can occur.
To sum up, the qualitative nature of stress and displacement fields occurring around the sharp edge tip area is determined by eigenvalue λ. Depending on a concentrator type (geometry, material parameters), local stress and displacement fields may have complex and real singularities, both single and multiple.
For forecasting of the fracture process, a quantitative determination of a singular stress field is necessary, and this is obtained by determination of stress intensity factors. The physical interpretation of SIFs was first presented by Irwin [13]. Definitions of SIFs, as well as methods of their calculation, depend on the qualitative nature of a singular stress field. For the sharp corners in a homogeneous material (eigenvalues always have real values-λ = λ r ∈ R) stress intensity factors can be defined as follows [1][2][3]: where K I and K II are factors for mode I and mode II loading, respectively, defined in the plane of symmetry of the corner. In the literature, e.g. [6,[14][15][16][17][18], for the issue of notches and cracks located on the interface of a bi-material structure, many various definitions of SIFs can be found. In this paper, SIFs were defined as in paper [17]. In a paired definition (2), modification allowing for its application for a sharp corner, of any tip angle, located on the interface of a bi-material structure or a homogeneous material (assuming δ = 0) was introduced: where a is a characteristic dimension that can be considered as a crack length or a notch height.
It is worth noting that in the literature for the issue of sharp corners, for which stress fields are described as eigenvalue λ = 1/2, factors K I and K II are often referred to as generalised stress intensity factors GSIF. In this paper, in order to avoid unnecessary confusion, the factors defined by formula (2) regardless of the value λ are always referred to as stress intensity factors SIFs.
When comparing the presented definitions of SIF-formula 1 and 2-it can be noticed that in the case of corners and cracks in a homogeneous material K I depends only on circumferential stresses, while K II depends only on tangential stresses. For the issue of a bi-material structure with a notch, the exponent λ, as already mentioned, can have complex values. Such fact causes that factors K I and K II depend on tangential and circumferential stresses at the same time. In consequence, this makes the analytical description of local stress and displacement fields and calculation procedures for SIFs difficult.
As regards methods of determination of values of SIFs, a closed-form solution exists only for the issue of an element with a crack (in a homogeneous material and a bi-material structure) assuming that the element is subject to tensile or shear load applied in infinity. In the case of finite dimensions of the element or the presence of a complex external load condition, in order to determine values of sought factors it is necessary to use numerical methods.
In the literature, many papers related to the determination of stress intensity factors can be found. Usually, the following analytical and numerical methods are used: I. Determination of sought factors based on change of strain energy related to a virtual expansion of the crack [19]. II. Determination of the value of integral J on the basis of numerical results [20]. III. Application of special finite elements [21]. IV. Comparison of stresses or displacements obtained from numerical solution (e.g. FEM) to a known analytical solution [22,23].
Not all the methods mentioned can be used for elements with a corner of freely defined boundary conditions (e.g. method I is applicable only for elements with cracks), and some have certain limitations. Using, for example, method II, it is possible to determine only the value of integral J which is described by a quadratic function of two variables (sought factors K I and K II ). Thus, as a result of a solution of such quadratic equation (or in the case of a structural notch-system of two equations) more than one set of correct solutions is obtained.
Regarding method III, when using commercial software for FEM modelling it is not always possible to simply apply own finite elements allowing for a direct determination of the sought factors.
Due to a fact that the purpose of this paper was to develop a universal method allowing for calculation of SIFs for sharp corners of freely defined boundary conditions (crack/notch in a homogeneous material or in a bi-material structure), method IV seems to be the justified choice. Use of such approach in most cases requires knowledge of an analytical description of the distribution of singular stress or displacement fields occurring around the corner tip area. This is cumbersome if the method is to be used for corners of various boundary conditions. Apart from that, it is not possible for all stress concentrators to directly compare analytical solutions (stresses/displacements) to a numerical solution. Such situation is encountered in cases of corners generating stress fields of complex singularities. This results from a fact that both factor K I and K II depend on tangential and normal stresses at the same time.
Therefore in this paper universal analytical functions of K I,II(r ) = f σ ϕFEM , τ r ϕFEM , λ type are developed. These functions, using extrapolation methods and in connection with a numerical solution (σ ϕFEM , τ r ϕFEM ), allowed for SIFs' determination of for any corner. Of course the developed functions are correlated with boundary conditions of the issue and the assumed SIFs' definition. In the developed functions K I,II(r ) , a detailed form of that is presented in Sect. 3, eigenvalue λ is responsible for meeting boundary conditions. The eigenvalues are determined on the basis of the so-called characteristic equations. Characteristic equation and methodology of its obtaining (based on eigenfunction technique) is presented in Sect. 2. In the same Section, also analytical formulas combining stress intensity factors are given, defined using formula (2), with singular stress fields. The method of determination of SIFs, along with the developed criterion, selection of area from which the data obtained from FEM are applied to analytical functions, is described in Sect. 3. In order to verify the applied method, for a selected group of geometrical and material structures with sharp corner, stress intensity factors were determined and compared to data available in the literature. The description of analysed specimens with a notch and the method of their modelling using FEM are presented in Sect. 4. Obtained results and their discussion are given in Sect. 5.

Identification of parameters describing the qualitative and quantitative nature of singular stress fields
As already mentioned, the purpose of this paper is to develop a universal method of determination of stress intensity factors. Universality of the method is that it can be used for the issues of sharp corners of freely defined boundary conditions. Thus, this Section analyses the issues of a bi-material structure with a structural notch located on the interface (Fig. 2). Analytical solutions describing qualitative and quantitative features of singular stress fields obtained for such case can be considered as universal solutions of the issue of corner of any geometry and any structural properties. It can be noticed easily that in the obtained solution when assuming some special assumptions, solutions for the following issues are produced: An analytical description was obtained using the approach described by the authors of papers [24,25]. In the applied method by using elasticity theory equations and Airy stress function, general asymptotic solutions are obtained describing individual components of stress and displacement fields: where r, ϕ-coordinates in a polar coordinate system located in corner tip, The method of determining Eq. (3) is discussed in detail, e.g. in [25]. Particular solution involves determination of eigenvalue λ and constants A i , B i , C i , D i . The sought values for the issue of structural notch (Fig. 2) were determined on the basis of both boundary and connection conditions: Using formulas (3), the main matrix of the above boundary conditions can be formulated as: where From zeroing of a determinant of the main matrix of boundary conditions-det M 8×8 = 0-eigenequation (5) was determined, subsequent radicals of which determine eigenvalues λ, where Constants A i , B i , C i , D i were determined by solving the below system of Eqs. (6) (resulting from boundary conditions): where In order to avoid a trivial solution (A i , B i , . . . = 0), it was necessary to modify the main matrix M 8×8 . This modification involved replacement of any two boundary conditions by the following formulas (7) resulting from the assumed definition of SIF (Eq. (2)): where (3), an analytical description of stress fields was obtained. Below is a special form of stress fields for angle ϕ = 0, i.e. in a plane for which SIF was defined, When the exponent is a real value (δ = 0) relations (8) will be simplified to the following form-σ ϕ1,2ϕ=0 = 1 √ 2π K II r λ r −1 . Full description of stress fields for the issue of structural notch can be found, e.g. in paper [9]. Relation (8), combining stress intensity factors and a stress field, has been used to formulate function K I,II(r ) = f σ ϕFEM , τ r ϕFEM , λ , which is discussed in the next Section.

Description of method used for the determination of SIFs
In this part of the paper forms of analytical functions used for the determination of SIFs are given, as well as a criterion for the selection of nodes from which numerical results of FEM are implemented to the analytical solution.
As regards the criterion for the selection of nodes, as is well known, if a diagram for-σ = Cr −b type stresses-in a double logarithmic system is linear (Fig. 3), gradient of the line is b.
For the issue of stress concentrators generating complex singularities, components of stresses are always dependent on factors K I and K II , at the same time, and therefore in order to use a regularity as shown in Fig. 3, it is necessary to determine [on the basis of Eq. (8)] the so-called combined stresses dependent only on a single factor: In Eqs. (9) and (10) by replacing σ, τ with circumferential and tangential stresses, respectively, obtained from FEM, it is possible to calculate gradient b for all nodes. Of course, when determining SIF only pairs of nodes of a gradient corresponding to the theoretical ones are taken into account, i.e. b = (λ r − 1) ± 0.01. The 'combined stresses' σ j(r,0) at distance r m and r m+1 from the tip of the notch can be formulated as: Using formulas (9-11) after simple mathematical transformations, functions (12) are obtained allowing for the determination of factors K I,II(r ) (at a certain distance from the tip of the notch): Factors K I,II(r ) are calculated by replacing-σ (r m+1 ) , σ (r m ) , τ (r m+1 ) τ (r m ) -in Eq. (12) respective values of circumferential and tangential stresses obtained from FEM in nodes with coordinates r m , r m+1 . It is worth noting that when the exponent λ is a real number (δ = 0) dependence, (12) is simplified to a form presented in paper [26].
The factor calculated using formula (12), for all m + 1 selected nodes of gradient b = (λ r − 1) ± 0.01, in theory should be identical. However, due to potential errors in numerical calculations the results can be slightly different. Thus in order to minimise an error it is possible to: -average results (13): -determine, using linear regression, stress intensity factors for r = 0 (14): Differences in the obtained results, according to the applied method of error minimizing-formulas (13) and (14), are discussed in Sect. 5.

Modelling of issues of sharp corners using finite element method
With the development of numerical methods, nowadays an FEM-based numerical approach is very often used. FEM, combined with analytical or experimental solutions, can be used, e.g., to study phenomena such as friction [27][28][29][30], liquid and heat flow [31][32][33] or piezoelectricity [34,35]. This method can also be used to estimate the basic parameters of fracture mechanics, such as the stress intensity factor. So elements with sharp corners were modelled using FEM by means of ANSYS environment (mechanical APDL) application. In numerical simulations, two types of specimens were modelled: -a rectangular plate with an edge sharp corner under uniaxial tension (Fig. 4a); -a rectangular plate with a centre sharp corner subjected to uniform tensile (Fig. 4b) or pure shear loading (Fig. 5).
For specimens with centre sharp corner, due to symmetry, only a right side of plates was modelled. In a plane of symmetry, the following was assumed: symmetry (for tension samples- Fig 4b) and anti-symmetry (for shear samples- Fig. 5) boundary conditions. Specimens were modelled as a mesh of quadrangular eight-node finite elements of an increased concentration around the tip area (Fig. 6).  In all simulations, specimens were subjected to the same loads σ = 1, assuming the plane stress condition. Moreover, it was assumed that Poisson's ratio is always ν 1 = ν 2 = 0.3. Other geometrical and material parameters were variable. By changing the relative stiffness value E 1 /E 2 , the issue of a sharp corner located either in a homogeneous material (E 1 /E 2 = 1) or in a bi-material structure (E 1 /E 2 = 1) was modelled. By selecting values of angular dimensions accordingly, it was possible to develop models of specimens with: (a) crack-α = γ = π, β = 0; (b) notch symmetrical in relation to interface-α = γ, β = 0.
Proportions of linear dimensions (a/w, h/w), according to a modelled issue, were selected so as to make it possible to compare the obtained results with data available in the literature.

Test results
In this part of the paper, stress intensity factors, calculated using the developed method, are presented. The obtained results, as far as it was possible, were compared to data available in the literature and the numerical solution obtained directly using the ANSYS application (this application allows for direct determination of SIF only for the issue of a crack in a homogeneous material).
Calculations were performed for the following cases of elements with sharp corners: -crack located in homogeneous material and bi-material structure; -symmetrical notch located in homogeneous material and bi-material structure.
For all cases, the calculated stress intensity factors were normalised in accordance with formula (15): Values of exponents-λ = λ r +iδ-were calculated on the basis of Eq. (5). Below are the obtained results.

Results for specimens with crack
Normalised values of stress intensity factors F j calculated for elements with central crack (Figs. 4b, 5, α = γ = π, β = 0) are given in Table 1. Values of exponents calculated from Eq. (5) for various relative stiffness values are consistent with data available in the literature [20]. Singular stress fields have single real (E 1 /E 2 = 1) or complex E 1 /E 2 = 1 singularity. In numerical models, the following geometrical dimensions were adopted: a = 1, a/w = 0.1, h = w. Such defined geometry is more or less consistent with the issue of a flat plate of infinite dimensions, with a central crack (the overall dimensions of the plate and the a/w ratio have been determined on the basis of numerical tests). Calculations were performed for two variants of load-tension and shear.
The obtained results were compared to the results produced directly by the ANSYS application and to the closed-form solution (16) [17]: where σ ∞ , τ ∞ is the load applied on the edges of a plate of infinite dimensions. By analysing the data given in Table 1, high consistency of the obtained results with reference results can be concluded. The maximum error was less than 2%.   Table 2 Values of normalised factors F j and eigenvalues λ for tension specimens with single-sided side crack  Table 2 presents results of the investigation of elements with a single-sided side crack (Fig. 4a, α = γ = π, β = 0). In numerical simulations, tension specimens of the following dimensions were modelled: h/w = 3/2, a = 1. Calculations were performed for various values a/w and E 1 /E 2 .
The obtained results were compared to the results produced directly by the ANSYS application and to data available in the literature. Similar to the case of a central crack, satisfactory consistency in solutions was achieved. Figure 7a presents normalised SIFs values (points) calculated using relation (12) for nodes (having gradient consistent with the theoretical one) of various coordinates r . A broken line is used to mark trend lines. The obtained trend line, both for factor K II (r ) and K II(r ) , is actually a horizontal line. Therefore it can be concluded that final SIFs' values can be determined by using both dependence (13), and (14). In this paper, dependence (14) is used. Such determined SIFs were applied to formula (8) and calculated, using dependences (9) and (10), the 'combined stresses'. The resulting analytical solution (solid line) was compared to FEM solution (point), which is shown in Fig. 7b. Nodes with a slope consistent with the theoretical one are marked in black. On the basis of the obtained results, consistency of both solutions in a range of about 10% a can be concluded.

Results for specimens with notch
In this part of the paper, results for elements with both side (one-sided) and central notches (symmetrical α = γ ) are presented. Specimens of tip angles β = 60 • and β = 90 • were analysed. Calculations were performed for various values of relative stiffness E 1 /E 2 . The eigenvalues λ determined using Eq. (5) are given in Table 3 and in a graphic form in Fig. 8. The obtained data are consistent with data available in the literature [37,38].  Table 3 The eigenvalues λ = λ r + iδ determined for various values E 1 /E 2 , (β = 60 • , β = 90 • ) On the basis of results presented in Table 3 and in Fig. 8, it can be concluded that for the analysed specimens, depending on the relative stiffness, the exponent λ within the range of 0 and 1 can take two real values or a single complex value. Thus local stress fields can be characterized by double real singularity or multiple complex singularity.
For the issue of a structural notch with a single complex singularity, stress intensity factors most often are defined in the same way as for the issue of an interfacial crack. When double real singularities occur, to describe local stress fields three approaches are used: I. Use of a single factor related to the first singular term λ r1 , defined as follows [39,40]: II. Use of two factors related to the first singular term λ r1 [41]: III. Use of two factors, whereas the first is related to circumferential stresses and first singular term λ r1 , and the second to tangential stresses and second singular term λ r2 [37]: Analysing formulas (17)- (19), it is easy to see that only K and K I are equal. K II coefficients cannot be compared with each other. Moreover, by using methods I and II, in the case of multiple singularities, stress factor(s) are calculated only for the first singular term [39], or terms of a higher order [40,42] additionally are determined.
In this paper, for the determination of stress intensity factors approach II was used. The obtained results were presented: -for tension specimens with side notch (Fig. 4a) in Table 4; -for tension specimens with central notch (Fig. 4b) in Table 5; -for pure shear specimens with central notch (Fig. 5) in Table 6.
In all calculations, the same proportions of linear and angular dimensions were adopted, as follows: a/w = 0.4, h/w = 2, a = 1, α = γ . For cases in which a double real singularity occurred, SIFs were calculated for both singular terms. Factors were determined in two stages by using the method described in papers [40,42]. At first factor for the first term K (λ r 1 ) I,II (with exclusion of impact of the second term) was calculated. Then circumferential and tangential stresses obtained from FEM were summed with stresses calculated from formula (8). In formula (8), previously determined factors for the first term, K (λ r 1 ) I,II were used. Stresses-σ, τ -summed in such way were applied to Eq. (12), and factors for the second singular term, K (λ r 2 ) I,II were determined. The summed stresses also were used in formulas (9,10). This allowed for selection of nodes of a gradient consistent with the theoretical, i.e. b = (λ r 2 − 1) ± 0.01. When calculating both K (λ r 1 ) I,II , and K (λ r 2 ) I,II it was assumed that in the case of lack of nodes with a slope consistent with the theoretical value-b = (λ ri − 1) ± 0.01-factors take zero-values.
The obtained results, when possible, were compared to data available in the literature. Very good consistency of results for the issue of a notch in a homogeneous material was obtained the maximum error between the obtained results and data available in the literature was 0.65%. As regards a notch in a bi-material structure, due to differences occurring in a method of defining stress intensity factors, in the literature it is hard to find data which the obtained results could be compared to. The issues of structural notches of the same geometrical and material parameters as in the presented article were analysed in paper [37]. The authors defined SIFs in accordance with relation (19). In the presented paper, definitions of SIFs described by a formula (18) were adopted. Therefore it was possible to compare only factors K I (as mentioned before factors K II are defined differently and cannot be compared directly). Differences in the obtained results are about 1.5% (with the exclusion of the case of a notch with tip angle β = 90 • and relative stiffness E 1 /E 2 = 10, 100 where the difference in results was about 4.5%). Table 4 Values of normalised factors F j for tension specimens with symmetrical one-sided side notch of tip angle β = 60 • and β = 90 •   Moreover in the paper variability of stress intensity factors calculated for selected nodes of various coordinate values r was estimated. The obtained results for the tension element with notch of tip angle β = 90 • and relative stiffness E 1 /E 2 = 10, are shown in Fig. 9. Alike the issue of interface crack, the trend line is a horizontal line, and thus SIFs can be determined by using both relations (13) and (14).
A structural notch, of such geometry and relative stiffness, generates stress fields having double real singularity. Therefore the paper analyses the influence of a number of singular terms used in an asymptotic description on the field of application of analytical solutions. Comparison of analytical solutions (8)(9)(10) to FEM solution is presented in Fig. 10.
By analysing the results presented in Fig. 10, it can be concluded that the use of only the first singular term allows for description of stress fields in an area near the notch tip (about 0.1% a). Use of both singular terms increases the area by about 70% a. What is more, by using only the second singular term in the analytical description, a solution deviating significantly from the results produced using FEM is obtained.

Summary and conclusions
The novel method of determination of stress intensity factors was developed. In this method, factors are determined on the basis of a specially developed analytical function to which stresses obtained from FEM solution are implemented. Moreover, a criterion allowing for selection of nodes from which numerical data are Fig. 9 The SIFs versus r/a for the tension specimens with one-sided notch (β = 90 • , E 1 /E 2 = 10): a calculated for the first singular term, and b calculated for the second singular term used, was developed. Use of the criterion increases accuracy in calculations and allows for a use of a mesh for division to finite elements of smaller concentration. The proposed method, in connection with the criterion of selection of nodes, allows for determination of SIFs for elements with stress concentrators generating both real and complex singularities. Moreover, it is possible to determine SIFs for cases of sharp corners with multiple singularities. Undoubtedly an advantage of the proposed method is that it does not require knowledge of an analytical description of singular stress/displacement fields or connection of SIF to energy parameters. The only required parameter determining boundary conditions of the analysed issue is eigenvalue λ. The method does not require the use of special finite elements either and can be used together with any commercial FEM software.
The paper presents stress intensity fields calculated for corners of various tip angles, relative stiffness, and relative notch height. The calculations were made for various alternatives of loading of specimens. The obtained results, when possible, were compared to data available in the literature. Very good consistency was obtained for the following types of stress concentrators: crack in a homogeneous material, interfacial crack, notch in a homogeneous material, notch located on the interface of a bi-material structure. Moreover, a part of the presented results is unique (no data available in the literature), e.g. elements with central notch in a bi-material structure subjected to shear and tension stresses. These results can be used for comparison by other scholars. )·σ −1 +τ rϕ (r,0,λ r2 ) (τ rϕ (r,0,λ r1 ) Fig. 10 The stresses versus r/a for tension specimens with one-sided notch (β = 90 • , E 1 /E 2 = 10): a determined for the first singular term (analytical solution-(9), (10)) and b determined for the first (dotted line), second (thin solid line) singular term, and summed values (thick solid line) for both singular terms (analytical solution (8)) by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.