DETERMINING STRESS INTENSITY FACTOR K I WITH EXTRAPOLATION METHOD

Original scientific paper In this article it is presented how to determine stress intensity factor from stress, using the extrapolation method for crack opening mode I. The simplicity of the extrapolation method is its advantage; basically, no special algorithms for calculating stress intensity factor are required. Stresses are determined with the finite element method; in relation to that the impact that the observance of the position of the integration point has on the preciseness of the calculated stress intensity factor and stress evaluation in the nodes of finite elements is analysed. The impact of using different trend polynomial functions, the size of finite elements and the number of points used for creating trend functions are analysed to find out the suitability of this method for determining the stress intensity factor KI.


Introduction
Fracture mechanics is the field of rigid body mechanics that describes the connection between external load and stress distribution near a crack tip.Griffith [1] was the first author who published the assumption that brittle materials have many tiny cracks, which produce concentrations of stresses that are sufficient to create cohesive strength in localised areas even in case of a much smaller nominal stress.Later, based on this, at certain assumptions, the stress field near the crack tip, which is placed in the infinite plate, was described, using lineaer elastic fracture mechanics (LEML) [2,3,4,5].It is typical of LEML that the crack tip plastic zone is very small.The stress intensity factor K I , which describes the stress-field intensity near the crack tip, is the main parameter of LEML.Based on K, it is possible to establish when the material will break due to critical crack length or when the crack will propagate [6,7,8,9].
There are several methods to determine stress intensity factors K that can be carried out with numerical procedures.Generally, they are divided into [10,11,12,13,14]: • substitution methods (crack opening method, use of stress field and displacement field near the crack tip) and • energy methods (energy release, energy release using J integral, J integral, virtual crack extension, crack opening and closing, weighting functions, etc.).
In relation to the crack opening method, the problem is solved, based on known displacements which are determined with numerical methods; then equations are used to calculate the stress intensity factor at different distances r from the crack tip for certain angles θ.
Based on the results concerning stresses and displacements in the nodes that are closest to the crack tip and using the equations for their theoretical distribution, it is possible to calculate appropriate stress intensity factors [15].In this way, it is possible to obtain equations that make it possible to use simple expressions and, thus, enable a quick calculation.
The energy release method is based on the connection between the stress intensity factor and the value of the released energy, defined on the basis of the Griffith theory.Numerical process of determining the stress intensity factor requires the calculation of deformation energy releases, based on which it is possible to define the corresponding stress intensity factor [12].
The crack closing method is based on the assumption that, when it comes to crack extension, the deformation energy release is proportional to the work required to close the crack to its initial length.The energy release method using J integral is based on the assumption that, in the LEML area, the energy release G can be equalized with J integral [16].Into this method, it is possible to integrate the VCE method [12,17,18], which can be used to calculate the deformation energy release, on the basis of which the stress intensity factor can be calculated.
The J integral method is very often used in numerical procedures of determining the stress intensity factor in relation to two-dimensional problems.This method is very suitable for numerical calculations as it makes it possible to select an arbitrary integration path.Another advantage of J integral as an energy method is that no precise description of crack tip singularity field is required to obtain good results.
The determination of stress intensity factor from stress, using the extrapolation method, is dealt with; in relation to this, the influences of the element size, the position of the integration point, the number of points for constructing a suitable curve and length of the crack are dealt with.The method is simple and, generally, it does not require special algorithms and procedures, which are integrated into various commercial finite elements software packages.

Stresses near a crack tip
The presence of a crack in the element that is exposed to external stress causes local redistribution of stresses; with the theory of elasticity it is possible to calculate stresses in any point (r, θ) near a crack tip (Fig. 1).
In case of crack opening mode I, the crack-tip stress field can be described, using the development in series [6,8,9]: 3 Using the extrapolation method to numerically determine stress intensity factor K from stresses When using the extrapolation method, the virtual value of the stress intensity factor is calculated based on stress distribution near a crack tip; then the value is extrapolated to the crack tip.In the crack tip, the relevant (trend) curve (usually a trend line) is placed [19], where points are positioned; in these points, stress values are numerically evaluated.In view of the Eq. ( 1), it is simplest to select that the angle θ = 0° as in such a case the virtual stress intensity factor K I * is calculated as follows [19]: After having calculated non-real stress intensity factors K I * in N E extrapolation points, relevant values K I (P) in the crack tip can be determined (Fig. 2).Generally speaking, it is possible to express the equation of the trend line as follows [11]: which presents the dependence between the virtual stress intensity factor and the distance from the crack tip.If the least squares method is used to determine the coefficients A and B in Eq. ( 3), they are determined as follows [19]: where Y i stands for the stress intensity factor in the extrapolation point at the distance r i from the crack tip.
After having determined the coefficients, the coefficient A gives the value of the stress intensity factor in the crack tip, whereas the coefficient B indicates the direction coefficient of the trend line.
4 Analysing the suitability of using the extrapolation method to determine the stress intensity factor KI from stress The analysis was performed, using the model of semi-infinite plate with a crack at the edge and with tension load.The semi-infinite plate with tension load and with a crack on the surface (Fig. 3) is one of the simplest and most basic examples in fracture mechanics with known analytical value of the stress intensity factor K I .
For such a case, the stress intensity factor K I is calculated as follows [20]: For the selected data: a = 10 mm and s 0 = 6 MPa, it applies that K I = 37,7 MPa•mm 1/2 .As in relation to FEM, stresses are evaluated in integration points (Fig. 4), the results for K I are first determined on the basis of the results for stresses, defined in integration points that are closest to the line along the crack [21].Fig. 4 shows their position in front of the crack tip; the position of an integration point is determined, based on the distance r and angle θ from the crack tip.Eight node quater-point plane strain elements were used around the crack tip in order to accurately model r −1/2 stresses behaviour in FEM.The other part of the model was descretised with regular eight node plane strain elements.
Then the impact upon the results is analysed with the position of the integration point in front of the crack tip being taken into account, in view of the distance from the crack tip, described with the coordinate r and not with the coordinate θ; it is assumed that points are located in the direction θ = 0°.Then the impact of the element size and of the crack length upon K is analysed.
To simulate a semi-infinite plate, a square with the side of 1000 mm is used.In view of the symmetry of the problem, half of the model is taken into consideration and relevant boundary conditions are prescribed, as presented in Fig. 5. Results for determining the stress intensity factor K I from stress, using the extrapolation method, are presented and compared with reference analytical solution determined from Eq. ( 6).When creating a trend curve, the number of points taken into account is based on the following: • Firstly, 24 points are taken into account, which means that 8 elements are considered (Fig. 4).• Then points are gradually taken away, from the crack tip away, until up to 9 points are taken away.
Results for the element size of 0,5 mm are presented in Fig. 6.Technical Gazette 23, 6(2016), 1673-1678 Figure 6 The influence of considering the position of the integration point, the number of points that are taken away when constructing a trend line and the level of the trend curve on the stress intensity factor It is evident from the results that the differences in results for K I are negligible when taking and when not taking the angle of θ into account.The differences in results are larger when four or five points are taken away.Also it is evident that the results are better when using a trend line than when using a 2 nd order trend polynomial.When a trend line is used, the difference is below 2 % with regard to reference solution K I = 37,7 MPa•mm 0,5 if at least four points are taken away.
If the values of virtual stress intensity factors K I in the nodes of elements along a straight line in front of the crack tip are calculated, comparable results are obtained when values of virtual K I for the first two elements are not taken into account.
The influence of the element size (1 mm, 0,5 mm and 0,25 mm) is dealt with by adding points, on the basis of which the trend line is constructed; the values for the first two elements in front of the crack tip are not taken into account as the analysis has shown that the calculated values of virtual crack intensity factors for these two elements are poor.A similar impact of not taking the results for the first element into consideration was established also in [16].It is evident from the results in Fig. 7 that, for all three dealt with element sizes, it is required to take at least 21 points into consideration when constructing a trend line to ensure that the deviation from the reference result is below ±2 %.
It is not always justifiable to use a trend line.This is particularly important when using a large number of points to construct a suitable curve at a larger distance from the crack tip.The values for the stress intensity factor presented in Tab. 1 refer to the interval length of 100 mm; the values for the first two elements are not taken into account.
It is evident from the results in Tab. 1 that they match the reference value best when using a 3 rd order polynomial, while deviations are largest when a line is used for interpolation.Also it is evident from the results that in such a case the number of points used for constructing a function does not influence the value as, at the interval length of 100 mm and at the element size of 0.5 mm, the number of points is two times larger than at the element size of 1 mm and, at the element size of 0,25 mm, the number of points is four times larger than at the element size of 1 mm.Our establishments are checked also for the example of calculating the stress intensity factor K I for the crack lengths of 10 mm, 20 mm and 30 mm, and for the element size of 0,625 mm.The results are presented in Tab. 2 and apply when using a line.It is evident from the results in Tab. 2 that the presented method can be used to reliably calculate the values of stress intensity factors K I also for various crack lengths and element sizes.

Conclusion
The suitability of the extrapolation method for determining stress intensity factor K I from stress is analysed.Based on the example of semi-infinite plate with tensile load and with a crack at the edge, for which an analytical result exists, the influence of the position of an integration point on the value of virtual stress intensity factors is analysed: its position defined with the polar angle is taken into account and, in another example, this angle is not taken into consideration.It is established that if the angle is not taken into consideration, this has no significant impact on the stress intensity factor, determined with extrapolation to the crack tip what simplifies using of this method.For determination of stress intensity factor from Eq. ( 1) based on stresses evaluated in integration points their position must be known, what is not always easy task, because many FE codes do not provide their position as an output result.Results from presented investigation showed that the method of extrapolation of stresses for determination of stress intensity factor can be used accurately enough in case when stresses are evaluated from nodes in front of the crack tip.Results obtained with a trend line used for extrapolation and results acquired with higher order trend polynomials used for extrapolation are compared.When using virtual stress intensity factors at a short interval, the stress intensity factor at a crack tip can be described well with extrapolation, using a trend line.Furthermore, the impact of the number of points used for the construction of a trend function (in this case of a line) is analysed.It is established that at least 21 points in a trend line must be considered to minimize the influence of the vicinity of a crack tip.If a very large number of points is used to construct a trend function, a line, according to our results, does not suffice for the process of determining the stress intensity factor from stress, using the extrapolation method, as the deviation from the reference results is too large.In such a case, the use of a 3 rd order trend polynomial has proved to be suitable for a more precise description of points of virtual stress intensity factors, using a function, and for a more precise determination of extrapolated value of the stress intensity factor at a crack tip.

Figure 1
Points on a line, where stress intensity factor is evaluated[19]

Figure 2
Figure 2Typical distribution of non-real stress intensity factors near a crack tip[19]

Figure 3
Figure 3 Semi-infinite plate with a surface crack

Figure 4 Figure 5
Figure 4The positons of integration points in view of the crack tip in a numerical model

Figure 7
Figure 7The influence of the element size and of the number of points when constructing a trend line on the stress intensity factor Results are shown in Fig.7.The name in the legend indicates the initial point that is taken into consideration when constructing a trend line.It is evident from the results in Fig.7that, for all three dealt with element sizes, it is required to take at least 21 points into consideration when constructing a trend line to ensure that the deviation from the reference result is below ±2 %.It is not always justifiable to use a trend line.This is particularly important when using a large number of

Table 1
The influence of the polynomial order on the stress intensity factor KI for various element sizes and for the interval length of 100 mm

Table 2
Stress intensity factor KI for various crack lengths